Index: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/lib_mpp.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/lib_mpp.F90	(revision 11646)
+++ NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/lib_mpp.F90	(revision 11646)
@@ -0,0 +1,2094 @@
+MODULE lib_mpp
+   !!======================================================================
+   !!                       ***  MODULE  lib_mpp  ***
+   !! Ocean numerics:  massively parallel processing library
+   !!=====================================================================
+   !! History :  OPA  !  1994  (M. Guyon, J. Escobar, M. Imbard)  Original code
+   !!            7.0  !  1997  (A.M. Treguier)  SHMEM additions
+   !!            8.0  !  1998  (M. Imbard, J. Escobar, L. Colombet ) SHMEM and MPI
+   !!                 !  1998  (J.M. Molines) Open boundary conditions
+   !!   NEMO     1.0  !  2003  (J.M. Molines, G. Madec)  F90, free form
+   !!                 !  2003  (J.M. Molines) add mpp_ini_north(_3d,_2d)
+   !!             -   !  2004  (R. Bourdalle Badie)  isend option in mpi
+   !!                 !  2004  (J.M. Molines) minloc, maxloc
+   !!             -   !  2005  (G. Madec, S. Masson)  npolj=5,6 F-point & ice cases
+   !!             -   !  2005  (R. Redler) Replacement of MPI_COMM_WORLD except for MPI_Abort
+   !!             -   !  2005  (R. Benshila, G. Madec)  add extra halo case
+   !!             -   !  2008  (R. Benshila) add mpp_ini_ice
+   !!            3.2  !  2009  (R. Benshila) SHMEM suppression, north fold in lbc_nfd
+   !!            3.2  !  2009  (O. Marti)    add mpp_ini_znl
+   !!            4.0  !  2011  (G. Madec)  move ctl_ routines from in_out_manager
+   !!            3.5  !  2012  (S.Mocavero, I. Epicoco) Add mpp_lnk_bdy_3d/2d routines to optimize the BDY comm.
+   !!            3.5  !  2013  (C. Ethe, G. Madec)  message passing arrays as local variables 
+   !!            3.5  !  2013  (S.Mocavero, I.Epicoco - CMCC) north fold optimizations
+   !!            3.6  !  2015  (O. Tintó and M. Castrillo - BSC) Added '_multiple' case for 2D lbc and max
+   !!            4.0  !  2017  (G. Madec) automatique allocation of array argument (use any 3rd dimension)
+   !!             -   !  2017  (G. Madec) create generic.h90 files to generate all lbc and north fold routines
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   ctl_stop      : update momentum and tracer Kz from a tke scheme
+   !!   ctl_warn      : initialization, namelist read, and parameters control
+   !!   ctl_opn       : Open file and check if required file is available.
+   !!   ctl_nam       : Prints informations when an error occurs while reading a namelist
+   !!   get_unit      : give the index of an unused logical unit
+   !!----------------------------------------------------------------------
+#if   defined key_mpp_mpi
+   !!----------------------------------------------------------------------
+   !!   'key_mpp_mpi'             MPI massively parallel processing library
+   !!----------------------------------------------------------------------
+   !!   lib_mpp_alloc : allocate mpp arrays
+   !!   mynode        : indentify the processor unit
+   !!   mpp_lnk       : interface (defined in lbclnk) for message passing of 2d or 3d arrays (mpp_lnk_2d, mpp_lnk_3d)
+   !!   mpp_lnk_icb   : interface for message passing of 2d arrays with extra halo for icebergs (mpp_lnk_2d_icb)
+   !!   mpprecv       :
+   !!   mppsend       :
+   !!   mppscatter    :
+   !!   mppgather     :
+   !!   mpp_min       : generic interface for mppmin_int , mppmin_a_int , mppmin_real, mppmin_a_real
+   !!   mpp_max       : generic interface for mppmax_int , mppmax_a_int , mppmax_real, mppmax_a_real
+   !!   mpp_sum       : generic interface for mppsum_int , mppsum_a_int , mppsum_real, mppsum_a_real
+   !!   mpp_minloc    :
+   !!   mpp_maxloc    :
+   !!   mppsync       :
+   !!   mppstop       :
+   !!   mpp_ini_north : initialisation of north fold
+   !!   mpp_lbc_north_icb : alternative to mpp_nfd for extra outer halo with icebergs
+   !!----------------------------------------------------------------------
+   USE dom_oce        ! ocean space and time domain
+   USE lbcnfd         ! north fold treatment
+   USE in_out_manager ! I/O manager
+
+   IMPLICIT NONE
+   PRIVATE
+
+   INTERFACE mpp_nfd
+      MODULE PROCEDURE   mpp_nfd_2d    , mpp_nfd_3d    , mpp_nfd_4d
+      MODULE PROCEDURE   mpp_nfd_2d_ptr, mpp_nfd_3d_ptr, mpp_nfd_4d_ptr
+   END INTERFACE
+
+   ! Interface associated to the mpp_lnk_... routines is defined in lbclnk
+   PUBLIC   mpp_lnk_2d    , mpp_lnk_3d    , mpp_lnk_4d
+   PUBLIC   mpp_lnk_2d_ptr, mpp_lnk_3d_ptr, mpp_lnk_4d_ptr
+   !
+!!gm  this should be useless
+   PUBLIC   mpp_nfd_2d    , mpp_nfd_3d    , mpp_nfd_4d
+   PUBLIC   mpp_nfd_2d_ptr, mpp_nfd_3d_ptr, mpp_nfd_4d_ptr
+!!gm end
+   !
+   PUBLIC   ctl_stop, ctl_warn, get_unit, ctl_opn, ctl_nam
+   PUBLIC   mynode, mppstop, mppsync, mpp_comm_free
+   PUBLIC   mpp_ini_north
+   PUBLIC   mpp_lnk_2d_icb
+   PUBLIC   mpp_lbc_north_icb
+   PUBLIC   mpp_min, mpp_max, mpp_sum, mpp_minloc, mpp_maxloc
+   PUBLIC   mpp_delay_max, mpp_delay_sum, mpp_delay_rcv
+   PUBLIC   mppscatter, mppgather
+   PUBLIC   mpp_ini_znl
+   PUBLIC   mppsend, mpprecv                          ! needed by TAM and ICB routines
+   PUBLIC   mpp_lnk_bdy_2d, mpp_lnk_bdy_3d, mpp_lnk_bdy_4d
+   
+   !! * Interfaces
+   !! define generic interface for these routine as they are called sometimes
+   !! with scalar arguments instead of array arguments, which causes problems
+   !! for the compilation on AIX system as well as NEC and SGI. Ok on COMPACQ
+   INTERFACE mpp_min
+      MODULE PROCEDURE mppmin_a_int, mppmin_int, mppmin_a_real, mppmin_real
+   END INTERFACE
+   INTERFACE mpp_max
+      MODULE PROCEDURE mppmax_a_int, mppmax_int, mppmax_a_real, mppmax_real
+   END INTERFACE
+   INTERFACE mpp_sum
+      MODULE PROCEDURE mppsum_a_int, mppsum_int, mppsum_a_real, mppsum_real,   &
+         &             mppsum_realdd, mppsum_a_realdd
+   END INTERFACE
+   INTERFACE mpp_minloc
+      MODULE PROCEDURE mpp_minloc2d ,mpp_minloc3d
+   END INTERFACE
+   INTERFACE mpp_maxloc
+      MODULE PROCEDURE mpp_maxloc2d ,mpp_maxloc3d
+   END INTERFACE
+
+   !! ========================= !!
+   !!  MPI  variable definition !!
+   !! ========================= !!
+!$AGRIF_DO_NOT_TREAT
+   INCLUDE 'mpif.h'
+!$AGRIF_END_DO_NOT_TREAT
+
+   LOGICAL, PUBLIC, PARAMETER ::   lk_mpp = .TRUE.    !: mpp flag
+
+   INTEGER, PARAMETER         ::   nprocmax = 2**10   ! maximun dimension (required to be a power of 2)
+
+   INTEGER, PUBLIC ::   mppsize        ! number of process
+   INTEGER, PUBLIC ::   mpprank        ! process number  [ 0 - size-1 ]
+!$AGRIF_DO_NOT_TREAT
+   INTEGER, PUBLIC ::   mpi_comm_oce   ! opa local communicator
+!$AGRIF_END_DO_NOT_TREAT
+
+   INTEGER :: MPI_SUMDD
+
+   ! variables used for zonal integration
+   INTEGER, PUBLIC ::   ncomm_znl       !: communicator made by the processors on the same zonal average
+   LOGICAL, PUBLIC ::   l_znl_root      !: True on the 'left'most processor on the same row
+   INTEGER         ::   ngrp_znl        !  group ID for the znl processors
+   INTEGER         ::   ndim_rank_znl   !  number of processors on the same zonal average
+   INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::   nrank_znl  ! dimension ndim_rank_znl, number of the procs into the same znl domain
+
+   ! North fold condition in mpp_mpi with jpni > 1 (PUBLIC for TAM)
+   INTEGER, PUBLIC ::   ngrp_world        !: group ID for the world processors
+   INTEGER, PUBLIC ::   ngrp_opa          !: group ID for the opa processors
+   INTEGER, PUBLIC ::   ngrp_north        !: group ID for the northern processors (to be fold)
+   INTEGER, PUBLIC ::   ncomm_north       !: communicator made by the processors belonging to ngrp_north
+   INTEGER, PUBLIC ::   ndim_rank_north   !: number of 'sea' processor in the northern line (can be /= jpni !)
+   INTEGER, PUBLIC ::   njmppmax          !: value of njmpp for the processors of the northern line
+   INTEGER, PUBLIC ::   north_root        !: number (in the comm_opa) of proc 0 in the northern comm
+   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE, SAVE ::   nrank_north   !: dimension ndim_rank_north
+
+   ! Type of send : standard, buffered, immediate
+   CHARACTER(len=1), PUBLIC ::   cn_mpi_send        !: type od mpi send/recieve (S=standard, B=bsend, I=isend)
+   LOGICAL         , PUBLIC ::   l_isend = .FALSE.  !: isend use indicator (T if cn_mpi_send='I')
+   INTEGER         , PUBLIC ::   nn_buffer          !: size of the buffer in case of mpi_bsend
+
+   ! Communications summary report
+   CHARACTER(len=128), DIMENSION(:), ALLOCATABLE ::   crname_lbc                   !: names of lbc_lnk calling routines
+   CHARACTER(len=128), DIMENSION(:), ALLOCATABLE ::   crname_glb                   !: names of global comm calling routines
+   CHARACTER(len=128), DIMENSION(:), ALLOCATABLE ::   crname_dlg                   !: names of delayed global comm calling routines
+   INTEGER, PUBLIC                               ::   ncom_stp = 0                 !: copy of time step # istp
+   INTEGER, PUBLIC                               ::   ncom_fsbc = 1                !: copy of sbc time step # nn_fsbc
+   INTEGER, PUBLIC                               ::   ncom_dttrc = 1               !: copy of top time step # nn_dttrc
+   INTEGER, PUBLIC                               ::   ncom_freq                    !: frequency of comm diagnostic
+   INTEGER, PUBLIC , DIMENSION(:,:), ALLOCATABLE ::   ncomm_sequence               !: size of communicated arrays (halos)
+   INTEGER, PARAMETER, PUBLIC                    ::   ncom_rec_max = 5000          !: max number of communication record
+   INTEGER, PUBLIC                               ::   n_sequence_lbc = 0           !: # of communicated arraysvia lbc
+   INTEGER, PUBLIC                               ::   n_sequence_glb = 0           !: # of global communications
+   INTEGER, PUBLIC                               ::   n_sequence_dlg = 0           !: # of delayed global communications
+   INTEGER, PUBLIC                               ::   numcom = -1                  !: logical unit for communicaton report
+   LOGICAL, PUBLIC                               ::   l_full_nf_update = .TRUE.    !: logical for a full (2lines) update of bc at North fold report
+   INTEGER,                    PARAMETER, PUBLIC ::   nbdelay = 2       !: number of delayed operations
+   !: name (used as id) of allreduce-delayed operations
+   ! Warning: we must use the same character length in an array constructor (at least for gcc compiler)
+   CHARACTER(len=32), DIMENSION(nbdelay), PUBLIC ::   c_delaylist = (/ 'cflice', 'fwb   ' /)
+   !: component name where the allreduce-delayed operation is performed
+   CHARACTER(len=3),  DIMENSION(nbdelay), PUBLIC ::   c_delaycpnt = (/ 'ICE'   , 'OCE' /)
+   TYPE, PUBLIC ::   DELAYARR
+      REAL(   wp), POINTER, DIMENSION(:) ::  z1d => NULL()
+      COMPLEX(wp), POINTER, DIMENSION(:) ::  y1d => NULL()
+   END TYPE DELAYARR
+   TYPE( DELAYARR ), DIMENSION(nbdelay), PUBLIC, SAVE  ::   todelay         !: must have SAVE for default initialization of DELAYARR
+   INTEGER,          DIMENSION(nbdelay), PUBLIC        ::   ndelayid = -1   !: mpi request id of the delayed operations
+
+   ! timing summary report
+   REAL(wp), DIMENSION(2), PUBLIC ::  waiting_time = 0._wp
+   REAL(wp)              , PUBLIC ::  compute_time = 0._wp, elapsed_time = 0._wp
+   
+   REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::   tampon   ! buffer in case of bsend
+
+   LOGICAL, PUBLIC ::   ln_nnogather                !: namelist control of northfold comms
+   LOGICAL, PUBLIC ::   l_north_nogather = .FALSE.  !: internal control of northfold comms
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   FUNCTION mynode( ldtxt, ldname, kumnam_ref, kumnam_cfg, kumond, kstop, localComm )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mynode  ***
+      !!
+      !! ** Purpose :   Find processor unit
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*),DIMENSION(:), INTENT(  out) ::   ldtxt        !
+      CHARACTER(len=*)             , INTENT(in   ) ::   ldname       !
+      INTEGER                      , INTENT(in   ) ::   kumnam_ref   ! logical unit for reference namelist
+      INTEGER                      , INTENT(in   ) ::   kumnam_cfg   ! logical unit for configuration namelist
+      INTEGER                      , INTENT(inout) ::   kumond       ! logical unit for namelist output
+      INTEGER                      , INTENT(inout) ::   kstop        ! stop indicator
+      INTEGER         , OPTIONAL   , INTENT(in   ) ::   localComm    !
+      !
+      INTEGER ::   mynode, ierr, code, ji, ii, ios
+      LOGICAL ::   mpi_was_called
+      !
+      NAMELIST/nammpp/ cn_mpi_send, nn_buffer, jpni, jpnj, ln_nnogather
+      !!----------------------------------------------------------------------
+      !
+      ii = 1
+      WRITE(ldtxt(ii),*)                                                                  ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) 'mynode : mpi initialisation'                                    ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '~~~~~~ '                                                        ;   ii = ii + 1
+      !
+      REWIND( kumnam_ref )              ! Namelist nammpp in reference namelist: mpi variables
+      READ  ( kumnam_ref, nammpp, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nammpp in reference namelist', lwp )
+      !
+      REWIND( kumnam_cfg )              ! Namelist nammpp in configuration namelist: mpi variables
+      READ  ( kumnam_cfg, nammpp, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nammpp in configuration namelist', lwp )
+      !
+      !                              ! control print
+      WRITE(ldtxt(ii),*) '   Namelist nammpp'                                             ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      mpi send type          cn_mpi_send = ', cn_mpi_send       ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      size exported buffer   nn_buffer   = ', nn_buffer,' bytes';   ii = ii + 1
+      !
+      IF( jpni < 1 .OR. jpnj < 1  ) THEN
+         WRITE(ldtxt(ii),*) '      jpni and jpnj will be calculated automatically' ;   ii = ii + 1
+      ELSE
+         WRITE(ldtxt(ii),*) '      processor grid extent in i         jpni = ',jpni       ;   ii = ii + 1
+         WRITE(ldtxt(ii),*) '      processor grid extent in j         jpnj = ',jpnj       ;   ii = ii + 1
+      ENDIF
+
+      WRITE(ldtxt(ii),*) '      avoid use of mpi_allgather at the north fold  ln_nnogather = ', ln_nnogather  ; ii = ii + 1
+
+      CALL mpi_initialized ( mpi_was_called, code )
+      IF( code /= MPI_SUCCESS ) THEN
+         DO ji = 1, SIZE(ldtxt)
+            IF( TRIM(ldtxt(ji)) /= '' )   WRITE(*,*) ldtxt(ji)      ! control print of mynode
+         END DO
+         WRITE(*, cform_err)
+         WRITE(*, *) 'lib_mpp: Error in routine mpi_initialized'
+         CALL mpi_abort( mpi_comm_world, code, ierr )
+      ENDIF
+
+      IF( mpi_was_called ) THEN
+         !
+         SELECT CASE ( cn_mpi_send )
+         CASE ( 'S' )                ! Standard mpi send (blocking)
+            WRITE(ldtxt(ii),*) '           Standard blocking mpi send (send)'             ;   ii = ii + 1
+         CASE ( 'B' )                ! Buffer mpi send (blocking)
+            WRITE(ldtxt(ii),*) '           Buffer blocking mpi send (bsend)'              ;   ii = ii + 1
+            IF( Agrif_Root() )   CALL mpi_init_oce( ldtxt, ii, ierr )
+         CASE ( 'I' )                ! Immediate mpi send (non-blocking send)
+            WRITE(ldtxt(ii),*) '           Immediate non-blocking send (isend)'           ;   ii = ii + 1
+            l_isend = .TRUE.
+         CASE DEFAULT
+            WRITE(ldtxt(ii),cform_err)                                                    ;   ii = ii + 1
+            WRITE(ldtxt(ii),*) '           bad value for cn_mpi_send = ', cn_mpi_send     ;   ii = ii + 1
+            kstop = kstop + 1
+         END SELECT
+         !
+      ELSEIF ( PRESENT(localComm) .AND. .NOT. mpi_was_called ) THEN
+         WRITE(ldtxt(ii),cform_err)                                                    ;   ii = ii + 1
+         WRITE(ldtxt(ii),*) ' lib_mpp: You cannot provide a local communicator '          ;   ii = ii + 1
+         WRITE(ldtxt(ii),*) '          without calling MPI_Init before ! '                ;   ii = ii + 1
+         kstop = kstop + 1
+      ELSE
+         SELECT CASE ( cn_mpi_send )
+         CASE ( 'S' )                ! Standard mpi send (blocking)
+            WRITE(ldtxt(ii),*) '           Standard blocking mpi send (send)'             ;   ii = ii + 1
+            CALL mpi_init( ierr )
+         CASE ( 'B' )                ! Buffer mpi send (blocking)
+            WRITE(ldtxt(ii),*) '           Buffer blocking mpi send (bsend)'              ;   ii = ii + 1
+            IF( Agrif_Root() )   CALL mpi_init_oce( ldtxt, ii, ierr )
+         CASE ( 'I' )                ! Immediate mpi send (non-blocking send)
+            WRITE(ldtxt(ii),*) '           Immediate non-blocking send (isend)'           ;   ii = ii + 1
+            l_isend = .TRUE.
+            CALL mpi_init( ierr )
+         CASE DEFAULT
+            WRITE(ldtxt(ii),cform_err)                                                    ;   ii = ii + 1
+            WRITE(ldtxt(ii),*) '           bad value for cn_mpi_send = ', cn_mpi_send     ;   ii = ii + 1
+            kstop = kstop + 1
+         END SELECT
+         !
+      ENDIF
+
+      IF( PRESENT(localComm) ) THEN
+         IF( Agrif_Root() ) THEN
+            mpi_comm_oce = localComm
+         ENDIF
+      ELSE
+         CALL mpi_comm_dup( mpi_comm_world, mpi_comm_oce, code)
+         IF( code /= MPI_SUCCESS ) THEN
+            DO ji = 1, SIZE(ldtxt)
+               IF( TRIM(ldtxt(ji)) /= '' )   WRITE(*,*) ldtxt(ji)      ! control print of mynode
+            END DO
+            WRITE(*, cform_err)
+            WRITE(*, *) ' lib_mpp: Error in routine mpi_comm_dup'
+            CALL mpi_abort( mpi_comm_world, code, ierr )
+         ENDIF
+      ENDIF
+
+#if defined key_agrif
+      IF( Agrif_Root() ) THEN
+         CALL Agrif_MPI_Init(mpi_comm_oce)
+      ELSE
+         CALL Agrif_MPI_set_grid_comm(mpi_comm_oce)
+      ENDIF
+#endif
+
+      CALL mpi_comm_rank( mpi_comm_oce, mpprank, ierr )
+      CALL mpi_comm_size( mpi_comm_oce, mppsize, ierr )
+      mynode = mpprank
+
+      IF( mynode == 0 ) THEN
+         CALL ctl_opn( kumond, TRIM(ldname), 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. , 1 )
+         WRITE(kumond, nammpp)      
+      ENDIF
+      !
+      CALL MPI_OP_CREATE(DDPDD_MPI, .TRUE., MPI_SUMDD, ierr)
+      !
+   END FUNCTION mynode
+
+   !!----------------------------------------------------------------------
+   !!                   ***  routine mpp_lnk_(2,3,4)d  ***
+   !!
+   !!   * Argument : dummy argument use in mpp_lnk_... routines
+   !!                ptab   :   array or pointer of arrays on which the boundary condition is applied
+   !!                cd_nat :   nature of array grid-points
+   !!                psgn   :   sign used across the north fold boundary
+   !!                kfld   :   optional, number of pt3d arrays
+   !!                cd_mpp :   optional, fill the overlap area only
+   !!                pval   :   optional, background value (used at closed boundaries)
+   !!----------------------------------------------------------------------
+   !
+   !                       !==  2D array and array of 2D pointer  ==!
+   !
+#  define DIM_2d
+#     define ROUTINE_LNK           mpp_lnk_2d
+#     include "mpp_lnk_generic.h90"
+#     undef ROUTINE_LNK
+#     define MULTI
+#     define ROUTINE_LNK           mpp_lnk_2d_ptr
+#     include "mpp_lnk_generic.h90"
+#     undef ROUTINE_LNK
+#     undef MULTI
+#  undef DIM_2d
+   !
+   !                       !==  3D array and array of 3D pointer  ==!
+   !
+#  define DIM_3d
+#     define ROUTINE_LNK           mpp_lnk_3d
+#     include "mpp_lnk_generic.h90"
+#     undef ROUTINE_LNK
+#     define MULTI
+#     define ROUTINE_LNK           mpp_lnk_3d_ptr
+#     include "mpp_lnk_generic.h90"
+#     undef ROUTINE_LNK
+#     undef MULTI
+#  undef DIM_3d
+   !
+   !                       !==  4D array and array of 4D pointer  ==!
+   !
+#  define DIM_4d
+#     define ROUTINE_LNK           mpp_lnk_4d
+#     include "mpp_lnk_generic.h90"
+#     undef ROUTINE_LNK
+#     define MULTI
+#     define ROUTINE_LNK           mpp_lnk_4d_ptr
+#     include "mpp_lnk_generic.h90"
+#     undef ROUTINE_LNK
+#     undef MULTI
+#  undef DIM_4d
+
+   !!----------------------------------------------------------------------
+   !!                   ***  routine mpp_nfd_(2,3,4)d  ***
+   !!
+   !!   * Argument : dummy argument use in mpp_nfd_... routines
+   !!                ptab   :   array or pointer of arrays on which the boundary condition is applied
+   !!                cd_nat :   nature of array grid-points
+   !!                psgn   :   sign used across the north fold boundary
+   !!                kfld   :   optional, number of pt3d arrays
+   !!                cd_mpp :   optional, fill the overlap area only
+   !!                pval   :   optional, background value (used at closed boundaries)
+   !!----------------------------------------------------------------------
+   !
+   !                       !==  2D array and array of 2D pointer  ==!
+   !
+#  define DIM_2d
+#     define ROUTINE_NFD           mpp_nfd_2d
+#     include "mpp_nfd_generic.h90"
+#     undef ROUTINE_NFD
+#     define MULTI
+#     define ROUTINE_NFD           mpp_nfd_2d_ptr
+#     include "mpp_nfd_generic.h90"
+#     undef ROUTINE_NFD
+#     undef MULTI
+#  undef DIM_2d
+   !
+   !                       !==  3D array and array of 3D pointer  ==!
+   !
+#  define DIM_3d
+#     define ROUTINE_NFD           mpp_nfd_3d
+#     include "mpp_nfd_generic.h90"
+#     undef ROUTINE_NFD
+#     define MULTI
+#     define ROUTINE_NFD           mpp_nfd_3d_ptr
+#     include "mpp_nfd_generic.h90"
+#     undef ROUTINE_NFD
+#     undef MULTI
+#  undef DIM_3d
+   !
+   !                       !==  4D array and array of 4D pointer  ==!
+   !
+#  define DIM_4d
+#     define ROUTINE_NFD           mpp_nfd_4d
+#     include "mpp_nfd_generic.h90"
+#     undef ROUTINE_NFD
+#     define MULTI
+#     define ROUTINE_NFD           mpp_nfd_4d_ptr
+#     include "mpp_nfd_generic.h90"
+#     undef ROUTINE_NFD
+#     undef MULTI
+#  undef DIM_4d
+
+
+   !!----------------------------------------------------------------------
+   !!                   ***  routine mpp_lnk_bdy_(2,3,4)d  ***
+   !!
+   !!   * Argument : dummy argument use in mpp_lnk_... routines
+   !!                ptab   :   array or pointer of arrays on which the boundary condition is applied
+   !!                cd_nat :   nature of array grid-points
+   !!                psgn   :   sign used across the north fold boundary
+   !!                kb_bdy :   BDY boundary set
+   !!                kfld   :   optional, number of pt3d arrays
+   !!----------------------------------------------------------------------
+   !
+   !                       !==  2D array and array of 2D pointer  ==!
+   !
+#  define DIM_2d
+#     define ROUTINE_BDY           mpp_lnk_bdy_2d
+#     include "mpp_bdy_generic.h90"
+#     undef ROUTINE_BDY
+#  undef DIM_2d
+   !
+   !                       !==  3D array and array of 3D pointer  ==!
+   !
+#  define DIM_3d
+#     define ROUTINE_BDY           mpp_lnk_bdy_3d
+#     include "mpp_bdy_generic.h90"
+#     undef ROUTINE_BDY
+#  undef DIM_3d
+   !
+   !                       !==  4D array and array of 4D pointer  ==!
+   !
+#  define DIM_4d
+#     define ROUTINE_BDY           mpp_lnk_bdy_4d
+#     include "mpp_bdy_generic.h90"
+#     undef ROUTINE_BDY
+#  undef DIM_4d
+
+   !!----------------------------------------------------------------------
+   !!
+   !!   load_array  &   mpp_lnk_2d_9    à generaliser a 3D et 4D
+   
+   
+   !!    mpp_lnk_sum_2d et 3D   ====>>>>>>   à virer du code !!!!
+   
+   
+   !!----------------------------------------------------------------------
+
+
+
+   SUBROUTINE mppsend( ktyp, pmess, kbytes, kdest, md_req )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mppsend  ***
+      !!
+      !! ** Purpose :   Send messag passing array
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(inout) ::   pmess(*)   ! array of real
+      INTEGER , INTENT(in   ) ::   kbytes     ! size of the array pmess
+      INTEGER , INTENT(in   ) ::   kdest      ! receive process number
+      INTEGER , INTENT(in   ) ::   ktyp       ! tag of the message
+      INTEGER , INTENT(in   ) ::   md_req     ! argument for isend
+      !!
+      INTEGER ::   iflag
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( cn_mpi_send )
+      CASE ( 'S' )                ! Standard mpi send (blocking)
+         CALL mpi_send ( pmess, kbytes, mpi_double_precision, kdest , ktyp, mpi_comm_oce        , iflag )
+      CASE ( 'B' )                ! Buffer mpi send (blocking)
+         CALL mpi_bsend( pmess, kbytes, mpi_double_precision, kdest , ktyp, mpi_comm_oce        , iflag )
+      CASE ( 'I' )                ! Immediate mpi send (non-blocking send)
+         ! be carefull, one more argument here : the mpi request identifier..
+         CALL mpi_isend( pmess, kbytes, mpi_double_precision, kdest , ktyp, mpi_comm_oce, md_req, iflag )
+      END SELECT
+      !
+   END SUBROUTINE mppsend
+
+
+   SUBROUTINE mpprecv( ktyp, pmess, kbytes, ksource )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpprecv  ***
+      !!
+      !! ** Purpose :   Receive messag passing array
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(inout) ::   pmess(*)   ! array of real
+      INTEGER , INTENT(in   ) ::   kbytes     ! suze of the array pmess
+      INTEGER , INTENT(in   ) ::   ktyp       ! Tag of the recevied message
+      INTEGER, OPTIONAL, INTENT(in) :: ksource    ! source process number
+      !!
+      INTEGER :: istatus(mpi_status_size)
+      INTEGER :: iflag
+      INTEGER :: use_source
+      !!----------------------------------------------------------------------
+      !
+      ! If a specific process number has been passed to the receive call,
+      ! use that one. Default is to use mpi_any_source
+      use_source = mpi_any_source
+      IF( PRESENT(ksource) )   use_source = ksource
+      !
+      CALL mpi_recv( pmess, kbytes, mpi_double_precision, use_source, ktyp, mpi_comm_oce, istatus, iflag )
+      !
+   END SUBROUTINE mpprecv
+
+
+   SUBROUTINE mppgather( ptab, kp, pio )
+      !!----------------------------------------------------------------------
+      !!                   ***  routine mppgather  ***
+      !!
+      !! ** Purpose :   Transfert between a local subdomain array and a work
+      !!     array which is distributed following the vertical level.
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)      , INTENT(in   ) ::   ptab   ! subdomain input array
+      INTEGER                           , INTENT(in   ) ::   kp     ! record length
+      REAL(wp), DIMENSION(jpi,jpj,jpnij), INTENT(  out) ::   pio    ! subdomain input array
+      !!
+      INTEGER :: itaille, ierror   ! temporary integer
+      !!---------------------------------------------------------------------
+      !
+      itaille = jpi * jpj
+      CALL mpi_gather( ptab, itaille, mpi_double_precision, pio, itaille     ,   &
+         &                            mpi_double_precision, kp , mpi_comm_oce, ierror )
+      !
+   END SUBROUTINE mppgather
+
+
+   SUBROUTINE mppscatter( pio, kp, ptab )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mppscatter  ***
+      !!
+      !! ** Purpose :   Transfert between awork array which is distributed
+      !!      following the vertical level and the local subdomain array.
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpnij)  ::   pio    ! output array
+      INTEGER                             ::   kp     ! Tag (not used with MPI
+      REAL(wp), DIMENSION(jpi,jpj)        ::   ptab   ! subdomain array input
+      !!
+      INTEGER :: itaille, ierror   ! temporary integer
+      !!---------------------------------------------------------------------
+      !
+      itaille = jpi * jpj
+      !
+      CALL mpi_scatter( pio, itaille, mpi_double_precision, ptab, itaille     ,   &
+         &                            mpi_double_precision, kp  , mpi_comm_oce, ierror )
+      !
+   END SUBROUTINE mppscatter
+
+   
+   SUBROUTINE mpp_delay_sum( cdname, cdelay, y_in, pout, ldlast, kcom )
+     !!----------------------------------------------------------------------
+      !!                   ***  routine mpp_delay_sum  ***
+      !!
+      !! ** Purpose :   performed delayed mpp_sum, the result is received on next call
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in   )               ::   cdname  ! name of the calling subroutine
+      CHARACTER(len=*), INTENT(in   )               ::   cdelay  ! name (used as id) of the delayed operation
+      COMPLEX(wp),      INTENT(in   ), DIMENSION(:) ::   y_in
+      REAL(wp),         INTENT(  out), DIMENSION(:) ::   pout
+      LOGICAL,          INTENT(in   )               ::   ldlast  ! true if this is the last time we call this routine
+      INTEGER,          INTENT(in   ), OPTIONAL     ::   kcom
+      !!
+      INTEGER ::   ji, isz
+      INTEGER ::   idvar
+      INTEGER ::   ierr, ilocalcomm
+      COMPLEX(wp), ALLOCATABLE, DIMENSION(:) ::   ytmp
+      !!----------------------------------------------------------------------
+      ilocalcomm = mpi_comm_oce
+      IF( PRESENT(kcom) )   ilocalcomm = kcom
+
+      isz = SIZE(y_in)
+      
+      IF( narea == 1 .AND. numcom == -1 ) CALL mpp_report( cdname, ld_dlg = .TRUE. )
+
+      idvar = -1
+      DO ji = 1, nbdelay
+         IF( TRIM(cdelay) == TRIM(c_delaylist(ji)) ) idvar = ji
+      END DO
+      IF ( idvar == -1 )   CALL ctl_stop( 'STOP',' mpp_delay_sum : please add a new delayed exchange for '//TRIM(cdname) )
+
+      IF ( ndelayid(idvar) == 0 ) THEN         ! first call    with restart: %z1d defined in iom_delay_rst
+         !                                       --------------------------
+         IF ( SIZE(todelay(idvar)%z1d) /= isz ) THEN                  ! Check dimension coherence
+            IF(lwp) WRITE(numout,*) ' WARNING: the nb of delayed variables in restart file is not the model one'
+            DEALLOCATE(todelay(idvar)%z1d)
+            ndelayid(idvar) = -1                                      ! do as if we had no restart
+         ELSE
+            ALLOCATE(todelay(idvar)%y1d(isz))
+            todelay(idvar)%y1d(:) = CMPLX(todelay(idvar)%z1d(:), 0., wp)   ! create %y1d, complex variable needed by mpi_sumdd
+         END IF
+      ENDIF
+      
+      IF( ndelayid(idvar) == -1 ) THEN         ! first call without restart: define %y1d and %z1d from y_in with blocking allreduce
+         !                                       --------------------------
+         ALLOCATE(todelay(idvar)%z1d(isz), todelay(idvar)%y1d(isz))   ! allocate also %z1d as used for the restart
+         CALL mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr )   ! get %y1d
+         todelay(idvar)%z1d(:) = REAL(todelay(idvar)%y1d(:), wp)      ! define %z1d from %y1d
+      ENDIF
+
+      IF( ndelayid(idvar) > 0 )   CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received
+
+      ! send back pout from todelay(idvar)%z1d defined at previous call
+      pout(:) = todelay(idvar)%z1d(:)
+
+      ! send y_in into todelay(idvar)%y1d with a non-blocking communication
+#if defined key_mpi2
+      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.)
+      CALL  mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ndelayid(idvar), ierr )
+      IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.)
+#else
+      CALL mpi_iallreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ndelayid(idvar), ierr )
+#endif
+
+   END SUBROUTINE mpp_delay_sum
+
+   
+   SUBROUTINE mpp_delay_max( cdname, cdelay, p_in, pout, ldlast, kcom )
+      !!----------------------------------------------------------------------
+      !!                   ***  routine mpp_delay_max  ***
+      !!
+      !! ** Purpose :   performed delayed mpp_max, the result is received on next call
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in   )                 ::   cdname  ! name of the calling subroutine
+      CHARACTER(len=*), INTENT(in   )                 ::   cdelay  ! name (used as id) of the delayed operation
+      REAL(wp),         INTENT(in   ), DIMENSION(:)   ::   p_in    ! 
+      REAL(wp),         INTENT(  out), DIMENSION(:)   ::   pout    ! 
+      LOGICAL,          INTENT(in   )                 ::   ldlast  ! true if this is the last time we call this routine
+      INTEGER,          INTENT(in   ), OPTIONAL       ::   kcom
+      !!
+      INTEGER ::   ji, isz
+      INTEGER ::   idvar
+      INTEGER ::   ierr, ilocalcomm
+      !!----------------------------------------------------------------------
+      ilocalcomm = mpi_comm_oce
+      IF( PRESENT(kcom) )   ilocalcomm = kcom
+
+      isz = SIZE(p_in)
+
+      IF( narea == 1 .AND. numcom == -1 ) CALL mpp_report( cdname, ld_dlg = .TRUE. )
+
+      idvar = -1
+      DO ji = 1, nbdelay
+         IF( TRIM(cdelay) == TRIM(c_delaylist(ji)) ) idvar = ji
+      END DO
+      IF ( idvar == -1 )   CALL ctl_stop( 'STOP',' mpp_delay_max : please add a new delayed exchange for '//TRIM(cdname) )
+
+      IF ( ndelayid(idvar) == 0 ) THEN         ! first call    with restart: %z1d defined in iom_delay_rst
+         !                                       --------------------------
+         IF ( SIZE(todelay(idvar)%z1d) /= isz ) THEN                  ! Check dimension coherence
+            IF(lwp) WRITE(numout,*) ' WARNING: the nb of delayed variables in restart file is not the model one'
+            DEALLOCATE(todelay(idvar)%z1d)
+            ndelayid(idvar) = -1                                      ! do as if we had no restart
+         END IF
+      ENDIF
+
+      IF( ndelayid(idvar) == -1 ) THEN         ! first call without restart: define %z1d from p_in with a blocking allreduce
+         !                                       --------------------------
+         ALLOCATE(todelay(idvar)%z1d(isz))
+         CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_DOUBLE_PRECISION, mpi_max, ilocalcomm, ierr )   ! get %z1d
+      ENDIF
+
+      IF( ndelayid(idvar) > 0 )   CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received
+
+      ! send back pout from todelay(idvar)%z1d defined at previous call
+      pout(:) = todelay(idvar)%z1d(:)
+
+      ! send p_in into todelay(idvar)%z1d with a non-blocking communication
+#if defined key_mpi2
+      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.)
+      CALL  mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_DOUBLE_PRECISION, mpi_max, ilocalcomm, ndelayid(idvar), ierr )
+      IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.)
+#else
+      CALL mpi_iallreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_DOUBLE_PRECISION, mpi_max, ilocalcomm, ndelayid(idvar), ierr )
+#endif
+
+   END SUBROUTINE mpp_delay_max
+
+   
+   SUBROUTINE mpp_delay_rcv( kid )
+      !!----------------------------------------------------------------------
+      !!                   ***  routine mpp_delay_rcv  ***
+      !!
+      !! ** Purpose :  force barrier for delayed mpp (needed for restart) 
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,INTENT(in   )      ::  kid 
+      INTEGER ::   ierr
+      !!----------------------------------------------------------------------
+      IF( ndelayid(kid) /= -2 ) THEN  
+#if ! defined key_mpi2
+         IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.)
+         CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr )                        ! make sure todelay(kid) is received
+         IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.)
+#endif
+         IF( ASSOCIATED(todelay(kid)%y1d) )   todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp)  ! define %z1d from %y1d
+         ndelayid(kid) = -2   ! add flag to know that mpi_wait was already called on kid
+      ENDIF
+   END SUBROUTINE mpp_delay_rcv
+
+   
+   !!----------------------------------------------------------------------
+   !!    ***  mppmax_a_int, mppmax_int, mppmax_a_real, mppmax_real  ***
+   !!   
+   !!----------------------------------------------------------------------
+   !!
+#  define OPERATION_MAX
+#  define INTEGER_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppmax_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppmax_a_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef INTEGER_TYPE
+!
+#  define REAL_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppmax_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppmax_a_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef REAL_TYPE
+#  undef OPERATION_MAX
+   !!----------------------------------------------------------------------
+   !!    ***  mppmin_a_int, mppmin_int, mppmin_a_real, mppmin_real  ***
+   !!   
+   !!----------------------------------------------------------------------
+   !!
+#  define OPERATION_MIN
+#  define INTEGER_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppmin_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppmin_a_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef INTEGER_TYPE
+!
+#  define REAL_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppmin_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppmin_a_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef REAL_TYPE
+#  undef OPERATION_MIN
+
+   !!----------------------------------------------------------------------
+   !!    ***  mppsum_a_int, mppsum_int, mppsum_a_real, mppsum_real  ***
+   !!   
+   !!   Global sum of 1D array or a variable (integer, real or complex)
+   !!----------------------------------------------------------------------
+   !!
+#  define OPERATION_SUM
+#  define INTEGER_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppsum_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppsum_a_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef INTEGER_TYPE
+!
+#  define REAL_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppsum_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppsum_a_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef REAL_TYPE
+#  undef OPERATION_SUM
+
+#  define OPERATION_SUM_DD
+#  define COMPLEX_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppsum_realdd
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppsum_a_realdd
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef COMPLEX_TYPE
+#  undef OPERATION_SUM_DD
+
+   !!----------------------------------------------------------------------
+   !!    ***  mpp_minloc2d, mpp_minloc3d, mpp_maxloc2d, mpp_maxloc3d
+   !!   
+   !!----------------------------------------------------------------------
+   !!
+#  define OPERATION_MINLOC
+#  define DIM_2d
+#     define ROUTINE_LOC           mpp_minloc2d
+#     include "mpp_loc_generic.h90"
+#     undef ROUTINE_LOC
+#  undef DIM_2d
+#  define DIM_3d
+#     define ROUTINE_LOC           mpp_minloc3d
+#     include "mpp_loc_generic.h90"
+#     undef ROUTINE_LOC
+#  undef DIM_3d
+#  undef OPERATION_MINLOC
+
+#  define OPERATION_MAXLOC
+#  define DIM_2d
+#     define ROUTINE_LOC           mpp_maxloc2d
+#     include "mpp_loc_generic.h90"
+#     undef ROUTINE_LOC
+#  undef DIM_2d
+#  define DIM_3d
+#     define ROUTINE_LOC           mpp_maxloc3d
+#     include "mpp_loc_generic.h90"
+#     undef ROUTINE_LOC
+#  undef DIM_3d
+#  undef OPERATION_MAXLOC
+
+   SUBROUTINE mppsync()
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mppsync  ***
+      !!
+      !! ** Purpose :   Massively parallel processors, synchroneous
+      !!
+      !!-----------------------------------------------------------------------
+      INTEGER :: ierror
+      !!-----------------------------------------------------------------------
+      !
+      CALL mpi_barrier( mpi_comm_oce, ierror )
+      !
+   END SUBROUTINE mppsync
+
+
+   SUBROUTINE mppstop( ldfinal, ld_force_abort ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mppstop  ***
+      !!
+      !! ** purpose :   Stop massively parallel processors method
+      !!
+      !!----------------------------------------------------------------------
+      LOGICAL, OPTIONAL, INTENT(in) :: ldfinal    ! source process number
+      LOGICAL, OPTIONAL, INTENT(in) :: ld_force_abort    ! source process number
+      LOGICAL ::   llfinal, ll_force_abort
+      INTEGER ::   info
+      !!----------------------------------------------------------------------
+      llfinal = .FALSE.
+      IF( PRESENT(ldfinal) ) llfinal = ldfinal
+      ll_force_abort = .FALSE.
+      IF( PRESENT(ld_force_abort) ) ll_force_abort = ld_force_abort
+      !
+      IF(ll_force_abort) THEN
+         CALL mpi_abort( MPI_COMM_WORLD )
+      ELSE
+         CALL mppsync
+         CALL mpi_finalize( info )
+      ENDIF
+      IF( .NOT. llfinal ) STOP 123
+      !
+   END SUBROUTINE mppstop
+
+
+   SUBROUTINE mpp_comm_free( kcom )
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kcom
+      !!
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+      !
+      CALL MPI_COMM_FREE(kcom, ierr)
+      !
+   END SUBROUTINE mpp_comm_free
+
+
+   SUBROUTINE mpp_ini_znl( kumout )
+      !!----------------------------------------------------------------------
+      !!               ***  routine mpp_ini_znl  ***
+      !!
+      !! ** Purpose :   Initialize special communicator for computing zonal sum
+      !!
+      !! ** Method  : - Look for processors in the same row
+      !!              - Put their number in nrank_znl
+      !!              - Create group for the znl processors
+      !!              - Create a communicator for znl processors
+      !!              - Determine if processor should write znl files
+      !!
+      !! ** output
+      !!      ndim_rank_znl = number of processors on the same row
+      !!      ngrp_znl = group ID for the znl processors
+      !!      ncomm_znl = communicator for the ice procs.
+      !!      n_znl_root = number (in the world) of proc 0 in the ice comm.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kumout   ! ocean.output logical units
+      !
+      INTEGER :: jproc      ! dummy loop integer
+      INTEGER :: ierr, ii   ! local integer
+      INTEGER, ALLOCATABLE, DIMENSION(:) ::   kwork
+      !!----------------------------------------------------------------------
+      !-$$     WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ngrp_world     : ', ngrp_world
+      !-$$     WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - mpi_comm_world : ', mpi_comm_world
+      !-$$     WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - mpi_comm_oce   : ', mpi_comm_oce
+      !
+      ALLOCATE( kwork(jpnij), STAT=ierr )
+      IF( ierr /= 0 ) THEN
+         WRITE(kumout, cform_err)
+         WRITE(kumout,*) 'mpp_ini_znl : failed to allocate 1D array of length jpnij'
+         CALL mppstop
+      ENDIF
+
+      IF( jpnj == 1 ) THEN
+         ngrp_znl  = ngrp_world
+         ncomm_znl = mpi_comm_oce
+      ELSE
+         !
+         CALL MPI_ALLGATHER ( njmpp, 1, mpi_integer, kwork, 1, mpi_integer, mpi_comm_oce, ierr )
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - kwork pour njmpp : ', kwork
+         !-$$        CALL flush(numout)
+         !
+         ! Count number of processors on the same row
+         ndim_rank_znl = 0
+         DO jproc=1,jpnij
+            IF ( kwork(jproc) == njmpp ) THEN
+               ndim_rank_znl = ndim_rank_znl + 1
+            ENDIF
+         END DO
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ndim_rank_znl : ', ndim_rank_znl
+         !-$$        CALL flush(numout)
+         ! Allocate the right size to nrank_znl
+         IF (ALLOCATED (nrank_znl)) DEALLOCATE(nrank_znl)
+         ALLOCATE(nrank_znl(ndim_rank_znl))
+         ii = 0
+         nrank_znl (:) = 0
+         DO jproc=1,jpnij
+            IF ( kwork(jproc) == njmpp) THEN
+               ii = ii + 1
+               nrank_znl(ii) = jproc -1
+            ENDIF
+         END DO
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - nrank_znl : ', nrank_znl
+         !-$$        CALL flush(numout)
+
+         ! Create the opa group
+         CALL MPI_COMM_GROUP(mpi_comm_oce,ngrp_opa,ierr)
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ngrp_opa : ', ngrp_opa
+         !-$$        CALL flush(numout)
+
+         ! Create the znl group from the opa group
+         CALL MPI_GROUP_INCL  ( ngrp_opa, ndim_rank_znl, nrank_znl, ngrp_znl, ierr )
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ngrp_znl ', ngrp_znl
+         !-$$        CALL flush(numout)
+
+         ! Create the znl communicator from the opa communicator, ie the pool of procs in the same row
+         CALL MPI_COMM_CREATE ( mpi_comm_oce, ngrp_znl, ncomm_znl, ierr )
+         !-$$        WRITE (numout,*) 'mpp_ini_znl ', nproc, ' - ncomm_znl ', ncomm_znl
+         !-$$        CALL flush(numout)
+         !
+      END IF
+
+      ! Determines if processor if the first (starting from i=1) on the row
+      IF ( jpni == 1 ) THEN
+         l_znl_root = .TRUE.
+      ELSE
+         l_znl_root = .FALSE.
+         kwork (1) = nimpp
+         CALL mpp_min ( 'lib_mpp', kwork(1), kcom = ncomm_znl)
+         IF ( nimpp == kwork(1)) l_znl_root = .TRUE.
+      END IF
+
+      DEALLOCATE(kwork)
+
+   END SUBROUTINE mpp_ini_znl
+
+
+   SUBROUTINE mpp_ini_north
+      !!----------------------------------------------------------------------
+      !!               ***  routine mpp_ini_north  ***
+      !!
+      !! ** Purpose :   Initialize special communicator for north folding
+      !!      condition together with global variables needed in the mpp folding
+      !!
+      !! ** Method  : - Look for northern processors
+      !!              - Put their number in nrank_north
+      !!              - Create groups for the world processors and the north processors
+      !!              - Create a communicator for northern processors
+      !!
+      !! ** output
+      !!      njmppmax = njmpp for northern procs
+      !!      ndim_rank_north = number of processors in the northern line
+      !!      nrank_north (ndim_rank_north) = number  of the northern procs.
+      !!      ngrp_world = group ID for the world processors
+      !!      ngrp_north = group ID for the northern processors
+      !!      ncomm_north = communicator for the northern procs.
+      !!      north_root = number (in the world) of proc 0 in the northern comm.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER ::   ierr
+      INTEGER ::   jjproc
+      INTEGER ::   ii, ji
+      !!----------------------------------------------------------------------
+      !
+      njmppmax = MAXVAL( njmppt )
+      !
+      ! Look for how many procs on the northern boundary
+      ndim_rank_north = 0
+      DO jjproc = 1, jpnij
+         IF( njmppt(jjproc) == njmppmax )   ndim_rank_north = ndim_rank_north + 1
+      END DO
+      !
+      ! Allocate the right size to nrank_north
+      IF (ALLOCATED (nrank_north)) DEALLOCATE(nrank_north)
+      ALLOCATE( nrank_north(ndim_rank_north) )
+
+      ! Fill the nrank_north array with proc. number of northern procs.
+      ! Note : the rank start at 0 in MPI
+      ii = 0
+      DO ji = 1, jpnij
+         IF ( njmppt(ji) == njmppmax   ) THEN
+            ii=ii+1
+            nrank_north(ii)=ji-1
+         END IF
+      END DO
+      !
+      ! create the world group
+      CALL MPI_COMM_GROUP( mpi_comm_oce, ngrp_world, ierr )
+      !
+      ! Create the North group from the world group
+      CALL MPI_GROUP_INCL( ngrp_world, ndim_rank_north, nrank_north, ngrp_north, ierr )
+      !
+      ! Create the North communicator , ie the pool of procs in the north group
+      CALL MPI_COMM_CREATE( mpi_comm_oce, ngrp_north, ncomm_north, ierr )
+      !
+   END SUBROUTINE mpp_ini_north
+
+
+   SUBROUTINE mpi_init_oce( ldtxt, ksft, code )
+      !!---------------------------------------------------------------------
+      !!                   ***  routine mpp_init.opa  ***
+      !!
+      !! ** Purpose :: export and attach a MPI buffer for bsend
+      !!
+      !! ** Method  :: define buffer size in namelist, if 0 no buffer attachment
+      !!            but classical mpi_init
+      !!
+      !! History :: 01/11 :: IDRIS initial version for IBM only
+      !!            08/04 :: R. Benshila, generalisation
+      !!---------------------------------------------------------------------
+      CHARACTER(len=*),DIMENSION(:), INTENT(  out) ::   ldtxt
+      INTEGER                      , INTENT(inout) ::   ksft
+      INTEGER                      , INTENT(  out) ::   code
+      INTEGER                                      ::   ierr, ji
+      LOGICAL                                      ::   mpi_was_called
+      !!---------------------------------------------------------------------
+      !
+      CALL mpi_initialized( mpi_was_called, code )      ! MPI initialization
+      IF ( code /= MPI_SUCCESS ) THEN
+         DO ji = 1, SIZE(ldtxt)
+            IF( TRIM(ldtxt(ji)) /= '' )   WRITE(*,*) ldtxt(ji)      ! control print of mynode
+         END DO
+         WRITE(*, cform_err)
+         WRITE(*, *) ' lib_mpp: Error in routine mpi_initialized'
+         CALL mpi_abort( mpi_comm_world, code, ierr )
+      ENDIF
+      !
+      IF( .NOT. mpi_was_called ) THEN
+         CALL mpi_init( code )
+         CALL mpi_comm_dup( mpi_comm_world, mpi_comm_oce, code )
+         IF ( code /= MPI_SUCCESS ) THEN
+            DO ji = 1, SIZE(ldtxt)
+               IF( TRIM(ldtxt(ji)) /= '' )   WRITE(*,*) ldtxt(ji)      ! control print of mynode
+            END DO
+            WRITE(*, cform_err)
+            WRITE(*, *) ' lib_mpp: Error in routine mpi_comm_dup'
+            CALL mpi_abort( mpi_comm_world, code, ierr )
+         ENDIF
+      ENDIF
+      !
+      IF( nn_buffer > 0 ) THEN
+         WRITE(ldtxt(ksft),*) 'mpi_bsend, buffer allocation of  : ', nn_buffer   ;   ksft = ksft + 1
+         ! Buffer allocation and attachment
+         ALLOCATE( tampon(nn_buffer), stat = ierr )
+         IF( ierr /= 0 ) THEN
+            DO ji = 1, SIZE(ldtxt)
+               IF( TRIM(ldtxt(ji)) /= '' )   WRITE(*,*) ldtxt(ji)      ! control print of mynode
+            END DO
+            WRITE(*, cform_err)
+            WRITE(*, *) ' lib_mpp: Error in ALLOCATE', ierr
+            CALL mpi_abort( mpi_comm_world, code, ierr )
+         END IF
+         CALL mpi_buffer_attach( tampon, nn_buffer, code )
+      ENDIF
+      !
+   END SUBROUTINE mpi_init_oce
+
+
+   SUBROUTINE DDPDD_MPI( ydda, yddb, ilen, itype )
+      !!---------------------------------------------------------------------
+      !!   Routine DDPDD_MPI: used by reduction operator MPI_SUMDD
+      !!
+      !!   Modification of original codes written by David H. Bailey
+      !!   This subroutine computes yddb(i) = ydda(i)+yddb(i)
+      !!---------------------------------------------------------------------
+      INTEGER                     , INTENT(in)    ::   ilen, itype
+      COMPLEX(wp), DIMENSION(ilen), INTENT(in)    ::   ydda
+      COMPLEX(wp), DIMENSION(ilen), INTENT(inout) ::   yddb
+      !
+      REAL(wp) :: zerr, zt1, zt2    ! local work variables
+      INTEGER  :: ji, ztmp           ! local scalar
+      !!---------------------------------------------------------------------
+      !
+      ztmp = itype   ! avoid compilation warning
+      !
+      DO ji=1,ilen
+      ! Compute ydda + yddb using Knuth's trick.
+         zt1  = real(ydda(ji)) + real(yddb(ji))
+         zerr = zt1 - real(ydda(ji))
+         zt2  = ((real(yddb(ji)) - zerr) + (real(ydda(ji)) - (zt1 - zerr))) &
+                + aimag(ydda(ji)) + aimag(yddb(ji))
+
+         ! The result is zt1 + zt2, after normalization.
+         yddb(ji) = cmplx ( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1),wp )
+      END DO
+      !
+   END SUBROUTINE DDPDD_MPI
+
+
+   SUBROUTINE mpp_lbc_north_icb( pt2d, cd_type, psgn, kextj)
+      !!---------------------------------------------------------------------
+      !!                   ***  routine mpp_lbc_north_icb  ***
+      !!
+      !! ** Purpose :   Ensure proper north fold horizontal bondary condition
+      !!              in mpp configuration in case of jpn1 > 1 and for 2d
+      !!              array with outer extra halo
+      !!
+      !! ** Method  :   North fold condition and mpp with more than one proc
+      !!              in i-direction require a specific treatment. We gather
+      !!              the 4+kextj northern lines of the global domain on 1
+      !!              processor and apply lbc north-fold on this sub array.
+      !!              Then we scatter the north fold array back to the processors.
+      !!              This routine accounts for an extra halo with icebergs
+      !!              and assumes ghost rows and columns have been suppressed.
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt2d     ! 2D array with extra halo
+      CHARACTER(len=1)        , INTENT(in   ) ::   cd_type  ! nature of pt3d grid-points
+      !                                                     !   = T ,  U , V , F or W -points
+      REAL(wp)                , INTENT(in   ) ::   psgn     ! = -1. the sign change across the
+      !!                                                    ! north fold, =  1. otherwise
+      INTEGER                 , INTENT(in   ) ::   kextj    ! Extra halo width at north fold
+      !
+      INTEGER ::   ji, jj, jr
+      INTEGER ::   ierr, itaille, ildi, ilei, iilb
+      INTEGER ::   ipj, ij, iproc
+      !
+      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE  ::  ztab_e, znorthloc_e
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  znorthgloio_e
+      !!----------------------------------------------------------------------
+      !
+      ipj=4
+      ALLOCATE(        ztab_e(jpiglo, 1-kextj:ipj+kextj)       ,       &
+     &            znorthloc_e(jpimax, 1-kextj:ipj+kextj)       ,       &
+     &          znorthgloio_e(jpimax, 1-kextj:ipj+kextj,jpni)    )
+      !
+      ztab_e(:,:)      = 0._wp
+      znorthloc_e(:,:) = 0._wp
+      !
+      ij = 1 - kextj
+      ! put the last ipj+2*kextj lines of pt2d into znorthloc_e 
+      DO jj = jpj - ipj + 1 - kextj , jpj + kextj
+         znorthloc_e(1:jpi,ij)=pt2d(1:jpi,jj)
+         ij = ij + 1
+      END DO
+      !
+      itaille = jpimax * ( ipj + 2*kextj )
+      !
+      IF( ln_timing ) CALL tic_tac(.TRUE.)
+      CALL MPI_ALLGATHER( znorthloc_e(1,1-kextj)    , itaille, MPI_DOUBLE_PRECISION,    &
+         &                znorthgloio_e(1,1-kextj,1), itaille, MPI_DOUBLE_PRECISION,    &
+         &                ncomm_north, ierr )
+      !
+      IF( ln_timing ) CALL tic_tac(.FALSE.)
+      !
+      DO jr = 1, ndim_rank_north            ! recover the global north array
+         iproc = nrank_north(jr) + 1
+         ildi = nldit (iproc)
+         ilei = nleit (iproc)
+         iilb = nimppt(iproc)
+         DO jj = 1-kextj, ipj+kextj
+            DO ji = ildi, ilei
+               ztab_e(ji+iilb-1,jj) = znorthgloio_e(ji,jj,jr)
+            END DO
+         END DO
+      END DO
+
+      ! 2. North-Fold boundary conditions
+      ! ----------------------------------
+      CALL lbc_nfd( ztab_e(:,1-kextj:ipj+kextj), cd_type, psgn, kextj )
+
+      ij = 1 - kextj
+      !! Scatter back to pt2d
+      DO jj = jpj - ipj + 1 - kextj , jpj + kextj
+         DO ji= 1, jpi
+            pt2d(ji,jj) = ztab_e(ji+nimpp-1,ij)
+         END DO
+         ij  = ij +1
+      END DO
+      !
+      DEALLOCATE( ztab_e, znorthloc_e, znorthgloio_e )
+      !
+   END SUBROUTINE mpp_lbc_north_icb
+
+
+   SUBROUTINE mpp_lnk_2d_icb( cdname, pt2d, cd_type, psgn, kexti, kextj )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpp_lnk_2d_icb  ***
+      !!
+      !! ** Purpose :   Message passing management for 2d array (with extra halo for icebergs)
+      !!                This routine receives a (1-kexti:jpi+kexti,1-kexti:jpj+kextj)
+      !!                array (usually (0:jpi+1, 0:jpj+1)) from lbc_lnk_icb calls.
+      !!
+      !! ** Method  :   Use mppsend and mpprecv function for passing mask
+      !!      between processors following neighboring subdomains.
+      !!            domain parameters
+      !!                    jpi    : first dimension of the local subdomain
+      !!                    jpj    : second dimension of the local subdomain
+      !!                    kexti  : number of columns for extra outer halo
+      !!                    kextj  : number of rows for extra outer halo
+      !!                    nbondi : mark for "east-west local boundary"
+      !!                    nbondj : mark for "north-south local boundary"
+      !!                    noea   : number for local neighboring processors
+      !!                    nowe   : number for local neighboring processors
+      !!                    noso   : number for local neighboring processors
+      !!                    nono   : number for local neighboring processors
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*)                                        , INTENT(in   ) ::   cdname      ! name of the calling subroutine
+      REAL(wp), DIMENSION(1-kexti:jpi+kexti,1-kextj:jpj+kextj), INTENT(inout) ::   pt2d     ! 2D array with extra halo
+      CHARACTER(len=1)                                        , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
+      REAL(wp)                                                , INTENT(in   ) ::   psgn     ! sign used across the north fold
+      INTEGER                                                 , INTENT(in   ) ::   kexti    ! extra i-halo width
+      INTEGER                                                 , INTENT(in   ) ::   kextj    ! extra j-halo width
+      !
+      INTEGER  ::   jl   ! dummy loop indices
+      INTEGER  ::   imigr, iihom, ijhom        ! local integers
+      INTEGER  ::   ipreci, iprecj             !   -       -
+      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
+      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
+      !!
+      REAL(wp), DIMENSION(1-kexti:jpi+kexti,nn_hls+kextj,2) ::   r2dns, r2dsn
+      REAL(wp), DIMENSION(1-kextj:jpj+kextj,nn_hls+kexti,2) ::   r2dwe, r2dew
+      !!----------------------------------------------------------------------
+
+      ipreci = nn_hls + kexti      ! take into account outer extra 2D overlap area
+      iprecj = nn_hls + kextj
+
+      IF( narea == 1 .AND. numcom == -1 ) CALL mpp_report( cdname, 1, 1, 1, ld_lbc = .TRUE. )
+
+      ! 1. standard boundary treatment
+      ! ------------------------------
+      ! Order matters Here !!!!
+      !
+      !                                      ! East-West boundaries
+      !                                           !* Cyclic east-west
+      IF( l_Iperio ) THEN
+         pt2d(1-kexti:     1   ,:) = pt2d(jpim1-kexti: jpim1 ,:)       ! east
+         pt2d(  jpi  :jpi+kexti,:) = pt2d(     2     :2+kexti,:)       ! west
+         !
+      ELSE                                        !* closed
+         IF( .NOT. cd_type == 'F' )   pt2d(  1-kexti   :nn_hls   ,:) = 0._wp    ! east except at F-point
+                                      pt2d(jpi-nn_hls+1:jpi+kexti,:) = 0._wp    ! west
+      ENDIF
+      !                                      ! North-South boundaries
+      IF( l_Jperio ) THEN                         !* cyclic (only with no mpp j-split)
+         pt2d(:,1-kextj:     1   ) = pt2d(:,jpjm1-kextj:  jpjm1)       ! north
+         pt2d(:,  jpj  :jpj+kextj) = pt2d(:,     2     :2+kextj)       ! south
+      ELSE                                        !* closed
+         IF( .NOT. cd_type == 'F' )   pt2d(:,  1-kextj   :nn_hls   ) = 0._wp    ! north except at F-point
+                                      pt2d(:,jpj-nn_hls+1:jpj+kextj) = 0._wp    ! south
+      ENDIF
+      !
+
+      ! north fold treatment
+      ! -----------------------
+      IF( npolj /= 0 ) THEN
+         !
+         SELECT CASE ( jpni )
+                   CASE ( 1 )     ;   CALL lbc_nfd          ( pt2d(1:jpi,1:jpj+kextj), cd_type, psgn, kextj )
+                   CASE DEFAULT   ;   CALL mpp_lbc_north_icb( pt2d(1:jpi,1:jpj+kextj), cd_type, psgn, kextj )
+         END SELECT
+         !
+      ENDIF
+
+      ! 2. East and west directions exchange
+      ! ------------------------------------
+      ! we play with the neigbours AND the row number because of the periodicity
+      !
+      SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
+      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
+         iihom = jpi-nreci-kexti
+         DO jl = 1, ipreci
+            r2dew(:,jl,1) = pt2d(nn_hls+jl,:)
+            r2dwe(:,jl,1) = pt2d(iihom +jl,:)
+         END DO
+      END SELECT
+      !
+      !                           ! Migrations
+      imigr = ipreci * ( jpj + 2*kextj )
+      !
+      IF( ln_timing ) CALL tic_tac(.TRUE.)
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         CALL mppsend( 2, r2dwe(1-kextj,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, r2dew(1-kextj,1,2), imigr, noea )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      CASE ( 0 )
+         CALL mppsend( 1, r2dew(1-kextj,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, r2dwe(1-kextj,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, r2dew(1-kextj,1,2), imigr, noea )
+         CALL mpprecv( 2, r2dwe(1-kextj,1,2), imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+      CASE ( 1 )
+         CALL mppsend( 1, r2dew(1-kextj,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, r2dwe(1-kextj,1,2), imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      END SELECT
+      !
+      IF( ln_timing ) CALL tic_tac(.FALSE.)
+      !
+      !                           ! Write Dirichlet lateral conditions
+      iihom = jpi - nn_hls
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         DO jl = 1, ipreci
+            pt2d(iihom+jl,:) = r2dew(:,jl,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, ipreci
+            pt2d(jl-kexti,:) = r2dwe(:,jl,2)
+            pt2d(iihom+jl,:) = r2dew(:,jl,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, ipreci
+            pt2d(jl-kexti,:) = r2dwe(:,jl,2)
+         END DO
+      END SELECT
+
+
+      ! 3. North and south directions
+      ! -----------------------------
+      ! always closed : we play only with the neigbours
+      !
+      IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
+         ijhom = jpj-nrecj-kextj
+         DO jl = 1, iprecj
+            r2dsn(:,jl,1) = pt2d(:,ijhom +jl)
+            r2dns(:,jl,1) = pt2d(:,nn_hls+jl)
+         END DO
+      ENDIF
+      !
+      !                           ! Migrations
+      imigr = iprecj * ( jpi + 2*kexti )
+      !
+      IF( ln_timing ) CALL tic_tac(.TRUE.)
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         CALL mppsend( 4, r2dsn(1-kexti,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, r2dns(1-kexti,1,2), imigr, nono )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      CASE ( 0 )
+         CALL mppsend( 3, r2dns(1-kexti,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, r2dsn(1-kexti,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, r2dns(1-kexti,1,2), imigr, nono )
+         CALL mpprecv( 4, r2dsn(1-kexti,1,2), imigr, noso )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+      CASE ( 1 )
+         CALL mppsend( 3, r2dns(1-kexti,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, r2dsn(1-kexti,1,2), imigr, noso )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      END SELECT
+      !
+      IF( ln_timing ) CALL tic_tac(.FALSE.)
+      !
+      !                           ! Write Dirichlet lateral conditions
+      ijhom = jpj - nn_hls
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         DO jl = 1, iprecj
+            pt2d(:,ijhom+jl) = r2dns(:,jl,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, iprecj
+            pt2d(:,jl-kextj) = r2dsn(:,jl,2)
+            pt2d(:,ijhom+jl) = r2dns(:,jl,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, iprecj
+            pt2d(:,jl-kextj) = r2dsn(:,jl,2)
+         END DO
+      END SELECT
+      !
+   END SUBROUTINE mpp_lnk_2d_icb
+
+
+   SUBROUTINE mpp_report( cdname, kpk, kpl, kpf, ld_lbc, ld_glb, ld_dlg )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpp_report  ***
+      !!
+      !! ** Purpose :   report use of mpp routines per time-setp
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*),           INTENT(in   ) ::   cdname      ! name of the calling subroutine
+      INTEGER         , OPTIONAL, INTENT(in   ) ::   kpk, kpl, kpf
+      LOGICAL         , OPTIONAL, INTENT(in   ) ::   ld_lbc, ld_glb, ld_dlg
+      !!
+      CHARACTER(len=128)                        ::   ccountname  ! name of a subroutine to count communications
+      CHARACTER(len=128)                        ::   ccom_report_name  ! name of communication report file
+      LOGICAL ::   ll_lbc, ll_glb, ll_dlg
+      INTEGER ::    ji,  jj,  jk,  jh, jf, jcount   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      ll_lbc = .FALSE.
+      IF( PRESENT(ld_lbc) ) ll_lbc = ld_lbc
+      ll_glb = .FALSE.
+      IF( PRESENT(ld_glb) ) ll_glb = ld_glb
+      ll_dlg = .FALSE.
+      IF( PRESENT(ld_dlg) ) ll_dlg = ld_dlg
+      !
+      ! find the smallest common frequency: default = frequency product, if multiple, choose the larger of the 2 frequency
+      IF( ncom_dttrc /= 1 )   CALL ctl_stop( 'STOP', 'mpp_report, ncom_dttrc /= 1 not coded...' ) 
+      ncom_freq = ncom_fsbc
+      !
+      IF ( ncom_stp == nit000+ncom_freq ) THEN   ! avoid to count extra communications in potential initializations at nit000
+         IF( ll_lbc ) THEN
+            IF( .NOT. ALLOCATED(ncomm_sequence) ) ALLOCATE( ncomm_sequence(ncom_rec_max,2) )
+            IF( .NOT. ALLOCATED(    crname_lbc) ) ALLOCATE(     crname_lbc(ncom_rec_max  ) )
+            n_sequence_lbc = n_sequence_lbc + 1
+            IF( n_sequence_lbc > ncom_rec_max ) CALL ctl_stop( 'STOP', 'lib_mpp, increase ncom_rec_max' )   ! deadlock
+            crname_lbc(n_sequence_lbc) = cdname     ! keep the name of the calling routine
+            ncomm_sequence(n_sequence_lbc,1) = kpk*kpl   ! size of 3rd and 4th dimensions
+            ncomm_sequence(n_sequence_lbc,2) = kpf       ! number of arrays to be treated (multi)
+         ENDIF
+         IF( ll_glb ) THEN
+            IF( .NOT. ALLOCATED(crname_glb) ) ALLOCATE( crname_glb(ncom_rec_max) )
+            n_sequence_glb = n_sequence_glb + 1
+            IF( n_sequence_glb > ncom_rec_max ) CALL ctl_stop( 'STOP', 'lib_mpp, increase ncom_rec_max' )   ! deadlock
+            crname_glb(n_sequence_glb) = cdname     ! keep the name of the calling routine
+         ENDIF
+         IF( ll_dlg ) THEN
+            IF( .NOT. ALLOCATED(crname_dlg) ) ALLOCATE( crname_dlg(ncom_rec_max) )
+            n_sequence_dlg = n_sequence_dlg + 1
+            IF( n_sequence_dlg > ncom_rec_max ) CALL ctl_stop( 'STOP', 'lib_mpp, increase ncom_rec_max' )   ! deadlock
+            crname_dlg(n_sequence_dlg) = cdname     ! keep the name of the calling routine
+         ENDIF
+      ELSE IF ( ncom_stp == nit000+2*ncom_freq ) THEN
+         ccom_report_name="communication_report.txt"
+         IF( TRIM(cxios_context) == 'sas' ) WRITE(ccom_report_name,'(A4,A24)') 'sas_', TRIM(ccom_report_name)
+         CALL ctl_opn( numcom, TRIM(ccom_report_name), 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
+         WRITE(numcom,*) ' '
+         WRITE(numcom,*) ' ------------------------------------------------------------'
+         WRITE(numcom,*) ' Communication pattern report (second oce+sbc+top time step):'
+         WRITE(numcom,*) ' ------------------------------------------------------------'
+         WRITE(numcom,*) ' '
+         WRITE(numcom,'(A,I4)') ' Exchanged halos : ', n_sequence_lbc
+         jj = 0; jk = 0; jf = 0; jh = 0
+         DO ji = 1, n_sequence_lbc
+            IF ( ncomm_sequence(ji,1) .GT. 1 ) jk = jk + 1
+            IF ( ncomm_sequence(ji,2) .GT. 1 ) jf = jf + 1
+            IF ( ncomm_sequence(ji,1) .GT. 1 .AND. ncomm_sequence(ji,2) .GT. 1 ) jj = jj + 1
+            jh = MAX (jh, ncomm_sequence(ji,1)*ncomm_sequence(ji,2))
+         END DO
+         WRITE(numcom,'(A,I3)') ' 3D Exchanged halos : ', jk
+         WRITE(numcom,'(A,I3)') ' Multi arrays exchanged halos : ', jf
+         WRITE(numcom,'(A,I3)') '   from which 3D : ', jj
+         WRITE(numcom,'(A,I10)') ' Array max size : ', jh*jpi*jpj
+         WRITE(numcom,*) ' '
+         WRITE(numcom,*) ' lbc_lnk called'
+         DO ji = 1, n_sequence_lbc - 1
+            IF ( crname_lbc(ji) /= 'already counted' ) THEN
+               ccountname = crname_lbc(ji)
+               crname_lbc(ji) = 'already counted'
+               jcount = 1
+               DO jj = ji + 1, n_sequence_lbc
+                  IF ( ccountname ==  crname_lbc(jj) ) THEN
+                     jcount = jcount + 1
+                     crname_lbc(jj) = 'already counted'
+                  END IF
+               END DO
+               WRITE(numcom,'(A, I4, A, A)') ' - ', jcount,' times by subroutine ', TRIM(ccountname)
+            END IF
+         END DO
+         IF ( crname_lbc(n_sequence_lbc) /= 'already counted' ) THEN
+            WRITE(numcom,'(A, I4, A, A)') ' - ', 1,' times by subroutine ', TRIM(crname_lbc(ncom_rec_max))
+         END IF
+         WRITE(numcom,*) ' '
+         IF ( n_sequence_glb > 0 ) THEN
+            WRITE(numcom,'(A,I4)') ' Global communications : ', n_sequence_glb
+            jj = 1
+            DO ji = 2, n_sequence_glb
+               IF( crname_glb(ji-1) /= crname_glb(ji) ) THEN
+                  WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_glb(ji-1))
+                  jj = 0
+               END IF
+               jj = jj + 1 
+            END DO
+            WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_glb(n_sequence_glb))
+            DEALLOCATE(crname_glb)
+         ELSE
+            WRITE(numcom,*) ' No MPI global communication '
+         ENDIF
+         WRITE(numcom,*) ' '
+         IF ( n_sequence_dlg > 0 ) THEN
+            WRITE(numcom,'(A,I4)') ' Delayed global communications : ', n_sequence_dlg
+            jj = 1
+            DO ji = 2, n_sequence_dlg
+               IF( crname_dlg(ji-1) /= crname_dlg(ji) ) THEN
+                  WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_dlg(ji-1))
+                  jj = 0
+               END IF
+               jj = jj + 1 
+            END DO
+            WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_dlg(n_sequence_dlg))
+            DEALLOCATE(crname_dlg)
+         ELSE
+            WRITE(numcom,*) ' No MPI delayed global communication '
+         ENDIF
+         WRITE(numcom,*) ' '
+         WRITE(numcom,*) ' -----------------------------------------------'
+         WRITE(numcom,*) ' '
+         DEALLOCATE(ncomm_sequence)
+         DEALLOCATE(crname_lbc)
+      ENDIF
+   END SUBROUTINE mpp_report
+
+   
+   SUBROUTINE tic_tac (ld_tic, ld_global)
+
+    LOGICAL,           INTENT(IN) :: ld_tic
+    LOGICAL, OPTIONAL, INTENT(IN) :: ld_global
+    REAL(wp), DIMENSION(2), SAVE :: tic_wt
+    REAL(wp),               SAVE :: tic_ct = 0._wp
+    INTEGER :: ii
+
+    IF( ncom_stp <= nit000 ) RETURN
+    IF( ncom_stp == nitend ) RETURN
+    ii = 1
+    IF( PRESENT( ld_global ) ) THEN
+       IF( ld_global ) ii = 2
+    END IF
+    
+    IF ( ld_tic ) THEN
+       tic_wt(ii) = MPI_Wtime()                                                    ! start count tic->tac (waiting time)
+       IF ( tic_ct > 0.0_wp ) compute_time = compute_time + MPI_Wtime() - tic_ct   ! cumulate count tac->tic
+    ELSE
+       waiting_time(ii) = waiting_time(ii) + MPI_Wtime() - tic_wt(ii)              ! cumulate count tic->tac
+       tic_ct = MPI_Wtime()                                                        ! start count tac->tic (waiting time)
+    ENDIF
+    
+   END SUBROUTINE tic_tac
+
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Default case:            Dummy module        share memory computing
+   !!----------------------------------------------------------------------
+   USE in_out_manager
+
+   INTERFACE mpp_sum
+      MODULE PROCEDURE mppsum_int, mppsum_a_int, mppsum_real, mppsum_a_real, mppsum_realdd, mppsum_a_realdd
+   END INTERFACE
+   INTERFACE mpp_max
+      MODULE PROCEDURE mppmax_a_int, mppmax_int, mppmax_a_real, mppmax_real
+   END INTERFACE
+   INTERFACE mpp_min
+      MODULE PROCEDURE mppmin_a_int, mppmin_int, mppmin_a_real, mppmin_real
+   END INTERFACE
+   INTERFACE mpp_minloc
+      MODULE PROCEDURE mpp_minloc2d ,mpp_minloc3d
+   END INTERFACE
+   INTERFACE mpp_maxloc
+      MODULE PROCEDURE mpp_maxloc2d ,mpp_maxloc3d
+   END INTERFACE
+
+   LOGICAL, PUBLIC, PARAMETER ::   lk_mpp = .FALSE.      !: mpp flag
+   LOGICAL, PUBLIC            ::   ln_nnogather          !: namelist control of northfold comms (needed here in case "key_mpp_mpi" is not used)
+   INTEGER, PUBLIC            ::   mpi_comm_oce          ! opa local communicator
+
+   INTEGER, PARAMETER, PUBLIC               ::   nbdelay = 0   ! make sure we don't enter loops: DO ji = 1, nbdelay
+   CHARACTER(len=32), DIMENSION(1), PUBLIC  ::   c_delaylist = 'empty'
+   CHARACTER(len=32), DIMENSION(1), PUBLIC  ::   c_delaycpnt = 'empty'
+   LOGICAL, PUBLIC                          ::   l_full_nf_update = .TRUE.
+   TYPE ::   DELAYARR
+      REAL(   wp), POINTER, DIMENSION(:) ::  z1d => NULL()
+      COMPLEX(wp), POINTER, DIMENSION(:) ::  y1d => NULL()
+   END TYPE DELAYARR
+   TYPE( DELAYARR ), DIMENSION(1), PUBLIC  ::   todelay              
+   INTEGER,  PUBLIC, DIMENSION(1)           ::   ndelayid = -1
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION lib_mpp_alloc(kumout)          ! Dummy function
+      INTEGER, INTENT(in) ::   kumout
+      lib_mpp_alloc = 0
+   END FUNCTION lib_mpp_alloc
+
+   FUNCTION mynode( ldtxt, ldname, kumnam_ref, knumnam_cfg,  kumond , kstop, localComm ) RESULT (function_value)
+      INTEGER, OPTIONAL            , INTENT(in   ) ::   localComm
+      CHARACTER(len=*),DIMENSION(:) ::   ldtxt
+      CHARACTER(len=*) ::   ldname
+      INTEGER ::   kumnam_ref, knumnam_cfg , kumond , kstop
+      IF( PRESENT( localComm ) ) mpi_comm_oce = localComm
+      function_value = 0
+      IF( .FALSE. )   ldtxt(:) = 'never done'
+      CALL ctl_opn( kumond, TRIM(ldname), 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. , 1 )
+   END FUNCTION mynode
+
+   SUBROUTINE mppsync                       ! Dummy routine
+   END SUBROUTINE mppsync
+
+   !!----------------------------------------------------------------------
+   !!    ***  mppmax_a_int, mppmax_int, mppmax_a_real, mppmax_real  ***
+   !!   
+   !!----------------------------------------------------------------------
+   !!
+#  define OPERATION_MAX
+#  define INTEGER_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppmax_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppmax_a_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef INTEGER_TYPE
+!
+#  define REAL_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppmax_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppmax_a_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef REAL_TYPE
+#  undef OPERATION_MAX
+   !!----------------------------------------------------------------------
+   !!    ***  mppmin_a_int, mppmin_int, mppmin_a_real, mppmin_real  ***
+   !!   
+   !!----------------------------------------------------------------------
+   !!
+#  define OPERATION_MIN
+#  define INTEGER_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppmin_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppmin_a_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef INTEGER_TYPE
+!
+#  define REAL_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppmin_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppmin_a_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef REAL_TYPE
+#  undef OPERATION_MIN
+
+   !!----------------------------------------------------------------------
+   !!    ***  mppsum_a_int, mppsum_int, mppsum_a_real, mppsum_real  ***
+   !!   
+   !!   Global sum of 1D array or a variable (integer, real or complex)
+   !!----------------------------------------------------------------------
+   !!
+#  define OPERATION_SUM
+#  define INTEGER_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppsum_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppsum_a_int
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef INTEGER_TYPE
+!
+#  define REAL_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppsum_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppsum_a_real
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef REAL_TYPE
+#  undef OPERATION_SUM
+
+#  define OPERATION_SUM_DD
+#  define COMPLEX_TYPE
+#  define DIM_0d
+#     define ROUTINE_ALLREDUCE           mppsum_realdd
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_0d
+#  define DIM_1d
+#     define ROUTINE_ALLREDUCE           mppsum_a_realdd
+#     include "mpp_allreduce_generic.h90"
+#     undef ROUTINE_ALLREDUCE
+#  undef DIM_1d
+#  undef COMPLEX_TYPE
+#  undef OPERATION_SUM_DD
+
+   !!----------------------------------------------------------------------
+   !!    ***  mpp_minloc2d, mpp_minloc3d, mpp_maxloc2d, mpp_maxloc3d
+   !!   
+   !!----------------------------------------------------------------------
+   !!
+#  define OPERATION_MINLOC
+#  define DIM_2d
+#     define ROUTINE_LOC           mpp_minloc2d
+#     include "mpp_loc_generic.h90"
+#     undef ROUTINE_LOC
+#  undef DIM_2d
+#  define DIM_3d
+#     define ROUTINE_LOC           mpp_minloc3d
+#     include "mpp_loc_generic.h90"
+#     undef ROUTINE_LOC
+#  undef DIM_3d
+#  undef OPERATION_MINLOC
+
+#  define OPERATION_MAXLOC
+#  define DIM_2d
+#     define ROUTINE_LOC           mpp_maxloc2d
+#     include "mpp_loc_generic.h90"
+#     undef ROUTINE_LOC
+#  undef DIM_2d
+#  define DIM_3d
+#     define ROUTINE_LOC           mpp_maxloc3d
+#     include "mpp_loc_generic.h90"
+#     undef ROUTINE_LOC
+#  undef DIM_3d
+#  undef OPERATION_MAXLOC
+
+   SUBROUTINE mpp_delay_sum( cdname, cdelay, y_in, pout, ldlast, kcom )
+      CHARACTER(len=*), INTENT(in   )               ::   cdname  ! name of the calling subroutine
+      CHARACTER(len=*), INTENT(in   )               ::   cdelay  ! name (used as id) of the delayed operation
+      COMPLEX(wp),      INTENT(in   ), DIMENSION(:) ::   y_in
+      REAL(wp),         INTENT(  out), DIMENSION(:) ::   pout
+      LOGICAL,          INTENT(in   )               ::   ldlast  ! true if this is the last time we call this routine
+      INTEGER,          INTENT(in   ), OPTIONAL     ::   kcom
+      !
+      pout(:) = REAL(y_in(:), wp)
+   END SUBROUTINE mpp_delay_sum
+
+   SUBROUTINE mpp_delay_max( cdname, cdelay, p_in, pout, ldlast, kcom )
+      CHARACTER(len=*), INTENT(in   )               ::   cdname  ! name of the calling subroutine
+      CHARACTER(len=*), INTENT(in   )               ::   cdelay  ! name (used as id) of the delayed operation
+      REAL(wp),         INTENT(in   ), DIMENSION(:) ::   p_in
+      REAL(wp),         INTENT(  out), DIMENSION(:) ::   pout
+      LOGICAL,          INTENT(in   )               ::   ldlast  ! true if this is the last time we call this routine
+      INTEGER,          INTENT(in   ), OPTIONAL     ::   kcom
+      !
+      pout(:) = p_in(:)
+   END SUBROUTINE mpp_delay_max
+
+   SUBROUTINE mpp_delay_rcv( kid )
+      INTEGER,INTENT(in   )      ::  kid 
+      WRITE(*,*) 'mpp_delay_rcv: You should not have seen this print! error?', kid
+   END SUBROUTINE mpp_delay_rcv
+   
+   SUBROUTINE mppstop( ldfinal, ld_force_abort )
+      LOGICAL, OPTIONAL, INTENT(in) :: ldfinal    ! source process number
+      LOGICAL, OPTIONAL, INTENT(in) :: ld_force_abort    ! source process number
+      STOP      ! non MPP case, just stop the run
+   END SUBROUTINE mppstop
+
+   SUBROUTINE mpp_ini_znl( knum )
+      INTEGER :: knum
+      WRITE(*,*) 'mpp_ini_znl: You should not have seen this print! error?', knum
+   END SUBROUTINE mpp_ini_znl
+
+   SUBROUTINE mpp_comm_free( kcom )
+      INTEGER :: kcom
+      WRITE(*,*) 'mpp_comm_free: You should not have seen this print! error?', kcom
+   END SUBROUTINE mpp_comm_free
+   
+#endif
+
+   !!----------------------------------------------------------------------
+   !!   All cases:         ctl_stop, ctl_warn, get_unit, ctl_opn, ctl_nam   routines
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE ctl_stop( cd1, cd2, cd3, cd4, cd5 ,   &
+      &                 cd6, cd7, cd8, cd9, cd10 )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE  stop_opa  ***
+      !!
+      !! ** Purpose :   print in ocean.outpput file a error message and
+      !!                increment the error number (nstop) by one.
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in), OPTIONAL ::  cd1, cd2, cd3, cd4, cd5
+      CHARACTER(len=*), INTENT(in), OPTIONAL ::  cd6, cd7, cd8, cd9, cd10
+      !!----------------------------------------------------------------------
+      !
+      nstop = nstop + 1
+
+      ! force to open ocean.output file
+      IF( numout == 6 ) CALL ctl_opn( numout, 'ocean.output', 'APPEND', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
+       
+      WRITE(numout,cform_err)
+      IF( PRESENT(cd1 ) )   WRITE(numout,*) TRIM(cd1)
+      IF( PRESENT(cd2 ) )   WRITE(numout,*) TRIM(cd2)
+      IF( PRESENT(cd3 ) )   WRITE(numout,*) TRIM(cd3)
+      IF( PRESENT(cd4 ) )   WRITE(numout,*) TRIM(cd4)
+      IF( PRESENT(cd5 ) )   WRITE(numout,*) TRIM(cd5)
+      IF( PRESENT(cd6 ) )   WRITE(numout,*) TRIM(cd6)
+      IF( PRESENT(cd7 ) )   WRITE(numout,*) TRIM(cd7)
+      IF( PRESENT(cd8 ) )   WRITE(numout,*) TRIM(cd8)
+      IF( PRESENT(cd9 ) )   WRITE(numout,*) TRIM(cd9)
+      IF( PRESENT(cd10) )   WRITE(numout,*) TRIM(cd10)
+
+                               CALL FLUSH(numout    )
+      IF( numstp     /= -1 )   CALL FLUSH(numstp    )
+      IF( numrun     /= -1 )   CALL FLUSH(numrun    )
+      IF( numevo_ice /= -1 )   CALL FLUSH(numevo_ice)
+      !
+      IF( cd1 == 'STOP' ) THEN
+         WRITE(numout,*)  'huge E-R-R-O-R : immediate stop'
+         CALL mppstop(ld_force_abort = .true.)
+      ENDIF
+      !
+   END SUBROUTINE ctl_stop
+
+
+   SUBROUTINE ctl_warn( cd1, cd2, cd3, cd4, cd5,   &
+      &                 cd6, cd7, cd8, cd9, cd10 )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE  stop_warn  ***
+      !!
+      !! ** Purpose :   print in ocean.outpput file a error message and
+      !!                increment the warning number (nwarn) by one.
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in), OPTIONAL ::  cd1, cd2, cd3, cd4, cd5
+      CHARACTER(len=*), INTENT(in), OPTIONAL ::  cd6, cd7, cd8, cd9, cd10
+      !!----------------------------------------------------------------------
+      !
+      nwarn = nwarn + 1
+      IF(lwp) THEN
+         WRITE(numout,cform_war)
+         IF( PRESENT(cd1 ) ) WRITE(numout,*) TRIM(cd1)
+         IF( PRESENT(cd2 ) ) WRITE(numout,*) TRIM(cd2)
+         IF( PRESENT(cd3 ) ) WRITE(numout,*) TRIM(cd3)
+         IF( PRESENT(cd4 ) ) WRITE(numout,*) TRIM(cd4)
+         IF( PRESENT(cd5 ) ) WRITE(numout,*) TRIM(cd5)
+         IF( PRESENT(cd6 ) ) WRITE(numout,*) TRIM(cd6)
+         IF( PRESENT(cd7 ) ) WRITE(numout,*) TRIM(cd7)
+         IF( PRESENT(cd8 ) ) WRITE(numout,*) TRIM(cd8)
+         IF( PRESENT(cd9 ) ) WRITE(numout,*) TRIM(cd9)
+         IF( PRESENT(cd10) ) WRITE(numout,*) TRIM(cd10)
+      ENDIF
+      CALL FLUSH(numout)
+      !
+   END SUBROUTINE ctl_warn
+
+
+   SUBROUTINE ctl_opn( knum, cdfile, cdstat, cdform, cdacce, klengh, kout, ldwp, karea )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ctl_opn  ***
+      !!
+      !! ** Purpose :   Open file and check if required file is available.
+      !!
+      !! ** Method  :   Fortan open
+      !!----------------------------------------------------------------------
+      INTEGER          , INTENT(  out) ::   knum      ! logical unit to open
+      CHARACTER(len=*) , INTENT(in   ) ::   cdfile    ! file name to open
+      CHARACTER(len=*) , INTENT(in   ) ::   cdstat    ! disposition specifier
+      CHARACTER(len=*) , INTENT(in   ) ::   cdform    ! formatting specifier
+      CHARACTER(len=*) , INTENT(in   ) ::   cdacce    ! access specifier
+      INTEGER          , INTENT(in   ) ::   klengh    ! record length
+      INTEGER          , INTENT(in   ) ::   kout      ! number of logical units for write
+      LOGICAL          , INTENT(in   ) ::   ldwp      ! boolean term for print
+      INTEGER, OPTIONAL, INTENT(in   ) ::   karea     ! proc number
+      !
+      CHARACTER(len=80) ::   clfile
+      INTEGER           ::   iost
+      !!----------------------------------------------------------------------
+      !
+      ! adapt filename
+      ! ----------------
+      clfile = TRIM(cdfile)
+      IF( PRESENT( karea ) ) THEN
+         IF( karea > 1 )   WRITE(clfile, "(a,'_',i6.6)") TRIM(clfile), karea-1
+      ENDIF
+#if defined key_agrif
+      IF( .NOT. Agrif_Root() )   clfile = TRIM(Agrif_CFixed())//'_'//TRIM(clfile)
+      knum=Agrif_Get_Unit()
+#else
+      knum=get_unit()
+#endif
+      IF( TRIM(cdfile) == '/dev/null' )   clfile = TRIM(cdfile)   ! force the use of /dev/null
+      !
+      iost=0
+      IF( cdacce(1:6) == 'DIRECT' )  THEN         ! cdacce has always more than 6 characters
+         OPEN( UNIT=knum, FILE=clfile, FORM=cdform, ACCESS=cdacce, STATUS=cdstat, RECL=klengh         , ERR=100, IOSTAT=iost )
+      ELSE IF( TRIM(cdstat) == 'APPEND' )  THEN   ! cdstat can have less than 6 characters
+         OPEN( UNIT=knum, FILE=clfile, FORM=cdform, ACCESS=cdacce, STATUS='UNKNOWN', POSITION='APPEND', ERR=100, IOSTAT=iost )
+      ELSE
+         OPEN( UNIT=knum, FILE=clfile, FORM=cdform, ACCESS=cdacce, STATUS=cdstat                      , ERR=100, IOSTAT=iost )
+      ENDIF
+      IF( iost /= 0 .AND. TRIM(clfile) == '/dev/null' ) &   ! for windows
+         &  OPEN(UNIT=knum,FILE='NUL', FORM=cdform, ACCESS=cdacce, STATUS=cdstat                      , ERR=100, IOSTAT=iost )   
+      IF( iost == 0 ) THEN
+         IF(ldwp) THEN
+            WRITE(kout,*) '     file   : ', TRIM(clfile),' open ok'
+            WRITE(kout,*) '     unit   = ', knum
+            WRITE(kout,*) '     status = ', cdstat
+            WRITE(kout,*) '     form   = ', cdform
+            WRITE(kout,*) '     access = ', cdacce
+            WRITE(kout,*)
+         ENDIF
+      ENDIF
+100   CONTINUE
+      IF( iost /= 0 ) THEN
+         IF(ldwp) THEN
+            WRITE(kout,*)
+            WRITE(kout,*) ' ===>>>> : bad opening file: ', TRIM(clfile)
+            WRITE(kout,*) ' =======   ===  '
+            WRITE(kout,*) '           unit   = ', knum
+            WRITE(kout,*) '           status = ', cdstat
+            WRITE(kout,*) '           form   = ', cdform
+            WRITE(kout,*) '           access = ', cdacce
+            WRITE(kout,*) '           iostat = ', iost
+            WRITE(kout,*) '           we stop. verify the file '
+            WRITE(kout,*)
+         ELSE  !!! Force writing to make sure we get the information - at least once - in this violent STOP!!
+            WRITE(*,*)
+            WRITE(*,*) ' ===>>>> : bad opening file: ', TRIM(clfile)
+            WRITE(*,*) ' =======   ===  '
+            WRITE(*,*) '           unit   = ', knum
+            WRITE(*,*) '           status = ', cdstat
+            WRITE(*,*) '           form   = ', cdform
+            WRITE(*,*) '           access = ', cdacce
+            WRITE(*,*) '           iostat = ', iost
+            WRITE(*,*) '           we stop. verify the file '
+            WRITE(*,*)
+         ENDIF
+         CALL FLUSH( kout ) 
+         STOP 'ctl_opn bad opening'
+      ENDIF
+      !
+   END SUBROUTINE ctl_opn
+
+
+   SUBROUTINE ctl_nam ( kios, cdnam, ldwp )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ctl_nam  ***
+      !!
+      !! ** Purpose :   Informations when error while reading a namelist
+      !!
+      !! ** Method  :   Fortan open
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(inout) ::   kios    ! IO status after reading the namelist
+      CHARACTER(len=*), INTENT(in   ) ::   cdnam   ! group name of namelist for which error occurs
+      CHARACTER(len=5)                ::   clios   ! string to convert iostat in character for print
+      LOGICAL         , INTENT(in   ) ::   ldwp    ! boolean term for print
+      !!----------------------------------------------------------------------
+      !
+      WRITE (clios, '(I5.0)')   kios
+      IF( kios < 0 ) THEN         
+         CALL ctl_warn( 'end of record or file while reading namelist '   &
+            &           // TRIM(cdnam) // ' iostat = ' // TRIM(clios) )
+      ENDIF
+      !
+      IF( kios > 0 ) THEN
+         CALL ctl_stop( 'misspelled variable in namelist '   &
+            &           // TRIM(cdnam) // ' iostat = ' // TRIM(clios) )
+      ENDIF
+      kios = 0
+      RETURN
+      !
+   END SUBROUTINE ctl_nam
+
+
+   INTEGER FUNCTION get_unit()
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  get_unit  ***
+      !!
+      !! ** Purpose :   return the index of an unused logical unit
+      !!----------------------------------------------------------------------
+      LOGICAL :: llopn
+      !!----------------------------------------------------------------------
+      !
+      get_unit = 15   ! choose a unit that is big enough then it is not already used in NEMO
+      llopn = .TRUE.
+      DO WHILE( (get_unit < 998) .AND. llopn )
+         get_unit = get_unit + 1
+         INQUIRE( unit = get_unit, opened = llopn )
+      END DO
+      IF( (get_unit == 999) .AND. llopn ) THEN
+         CALL ctl_stop( 'get_unit: All logical units until 999 are used...' )
+         get_unit = -1
+      ENDIF
+      !
+   END FUNCTION get_unit
+
+   !!----------------------------------------------------------------------
+END MODULE lib_mpp
Index: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/nemogcm.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/nemogcm.F90	(revision 11620)
+++ NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/nemogcm.F90	(revision 11646)
@@ -52,4 +52,7 @@
    CHARACTER(lc) ::   cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
 
+#if defined key_mpp_mpi
+   INCLUDE 'mpif.h'
+#endif
    !!----------------------------------------------------------------------
    !! NEMO/SAS 4.0 , NEMO Consortium (2018)
@@ -113,4 +116,9 @@
       !
       DO WHILE( istp <= nitend .AND. nstop == 0 )
+#if defined key_mpp_mpi
+            ncom_stp = istp
+            IF ( istp == ( nit000 + 1 ) ) elapsed_time = MPI_Wtime()
+            IF ( istp ==         nitend ) elapsed_time = MPI_Wtime() - elapsed_time
+#endif
          CALL stp
          istp = istp + 1
@@ -128,4 +136,9 @@
          !
          DO WHILE( istp <= nitend .AND. nstop == 0 )
+#if defined key_mpp_mpi
+            ncom_stp = istp
+            IF ( istp == ( nit000 + 1 ) ) elapsed_time = MPI_Wtime()
+            IF ( istp ==         nitend ) elapsed_time = MPI_Wtime() - elapsed_time
+#endif
             CALL stp        ( istp ) 
             istp = istp + 1
Index: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/sbcssm.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/sbcssm.F90	(revision 11620)
+++ NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/sbcssm.F90	(revision 11646)
@@ -178,5 +178,7 @@
       !!----------------------------------------------------------------------
       !
-      IF( ln_rstart .AND. nn_components == jp_iam_sas )   RETURN
+      !EM disable reading if not restart 
+      !EM IF( ln_rstart .AND. nn_components == jp_iam_sas )   RETURN
+      IF( nn_components == jp_iam_sas )   RETURN
       !
       IF(lwp) THEN
Index: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/stpctl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/stpctl.F90	(revision 11620)
+++ NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/stpctl.F90	(revision 11646)
@@ -60,4 +60,5 @@
       LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns
       CHARACTER(len=20) :: clname
+      CHARACTER(len=20) :: cfts_name, cfrs_name, cfrsnc_name
       !!----------------------------------------------------------------------
       !
@@ -70,12 +71,22 @@
          WRITE(numout,*) '~~~~~~~'
          !                                ! open time.step file
-         IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         cfts_name='time.step'
+         cfrs_name='run.stat'
+         cfrsnc_name='run.stat.nc'
+
+         IF( TRIM(cxios_context) == 'sas' ) THEN
+            WRITE(cfts_name,'(A4,A9)') 'sas_', TRIM(cfts_name)
+            WRITE(cfrs_name,'(A4,A8)') 'sas_', TRIM(cfrs_name)
+            WRITE(cfrsnc_name,'(A4,A11)') 'sas_', TRIM(cfrsnc_name)
+         ENDIF
+
+         IF( lwm ) CALL ctl_opn( numstp, TRIM(cfts_name), 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
          !                                ! open run.stat file(s) at start whatever
          !                                ! the value of sn_cfctl%ptimincr
          IF( lwm .AND. ( ln_ctl .OR. sn_cfctl%l_runstat ) ) THEN
-            CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
-            clname = 'run.stat.nc'
+            CALL ctl_opn( numrun, TRIM(cfrs_name), 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            clname = TRIM(cfrsnc_name)
             IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
-            istatus = NF90_CREATE( 'run.stat.nc', NF90_CLOBBER, idrun )
+            istatus = NF90_CREATE( TRIM(cfrsnc_name), NF90_CLOBBER, idrun )
             istatus = NF90_DEF_DIM( idrun, 'time'     , NF90_UNLIMITED, idtime )
             istatus = NF90_DEF_VAR( idrun, 'vt_i_max' , NF90_DOUBLE, (/ idtime /), idssh )
Index: NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/timing.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/timing.F90	(revision 11646)
+++ NEMO/branches/UKMO/NEMO_4.0_mirror_SI3_decoupled/src/SAS/timing.F90	(revision 11646)
@@ -0,0 +1,870 @@
+MODULE timing
+   !!========================================================================
+   !!                     ***  MODULE  timing  ***
+   !!========================================================================
+   !! History : 4.0  ! 2001-05  (R. Benshila)   
+   !!------------------------------------------------------------------------
+
+   !!------------------------------------------------------------------------
+   !!   timming_init    : initialize timing process 
+   !!   timing_start    : start Timer
+   !!   timing_stop     : stop  Timer
+   !!   timing_reset    : end timing variable creation
+   !!   timing_finalize : compute stats and write output in calling w*_info 
+   !!   timing_ini_var  : create timing variables 
+   !!   timing_listing  : print instumented subroutines in ocean.output
+   !!   wcurrent_info   : compute and print detailed stats on the current CPU
+   !!   wave_info       : compute and print averaged statson all processors
+   !!   wmpi_info       : compute and write global stats  
+   !!   supress         : suppress an element of the timing linked list  
+   !!   insert          : insert an element of the timing linked list  
+   !!------------------------------------------------------------------------
+   USE in_out_manager  ! I/O manager 
+   USE dom_oce         ! ocean domain
+   USE lib_mpp          
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   timing_init, timing_finalize   ! called in nemogcm module 
+   PUBLIC   timing_reset                   ! called in step module 
+   PUBLIC   timing_start, timing_stop      ! called in each routine to time 
+   
+#if defined key_mpp_mpi
+   INCLUDE 'mpif.h'
+#endif
+
+   ! Variables for fine grain timing
+   TYPE timer
+      CHARACTER(LEN=20)  :: cname
+      CHARACTER(LEN=20)  :: surname
+      INTEGER :: rank
+      REAL(wp)  :: t_cpu, t_clock, tsum_cpu, tsum_clock, tmax_cpu, tmax_clock, tmin_cpu, tmin_clock, tsub_cpu, tsub_clock
+      INTEGER :: ncount, ncount_max, ncount_rate  
+      INTEGER :: niter
+      LOGICAL :: l_tdone
+      TYPE(timer), POINTER :: next => NULL()
+      TYPE(timer), POINTER :: prev => NULL()
+      TYPE(timer), POINTER :: parent_section => NULL()
+   END TYPE timer
+    
+   TYPE alltimer
+      CHARACTER(LEN=20), DIMENSION(:), POINTER :: cname => NULL()
+      REAL(wp), DIMENSION(:), POINTER :: tsum_cpu   => NULL()
+      REAL(wp), DIMENSION(:), POINTER :: tsum_clock => NULL()
+      INTEGER, DIMENSION(:), POINTER :: niter => NULL()
+      TYPE(alltimer), POINTER :: next => NULL()
+      TYPE(alltimer), POINTER :: prev => NULL()
+   END TYPE alltimer 
+ 
+   TYPE(timer), POINTER :: s_timer_root => NULL()
+   TYPE(timer), POINTER :: s_timer      => NULL()
+   TYPE(timer), POINTER :: s_timer_old      => NULL()
+
+   TYPE(timer), POINTER :: s_wrk        => NULL()
+   REAL(wp) :: t_overclock, t_overcpu
+   LOGICAL :: l_initdone = .FALSE.
+   INTEGER :: nsize
+   
+   ! Variables for coarse grain timing
+   REAL(wp) :: tot_etime, tot_ctime
+   REAL(kind=wp), DIMENSION(2)     :: t_elaps, t_cpu
+   REAL(wp), ALLOCATABLE, DIMENSION(:) :: all_etime, all_ctime
+   INTEGER :: nfinal_count, ncount, ncount_rate, ncount_max
+   INTEGER, DIMENSION(8)           :: nvalues
+   CHARACTER(LEN=8), DIMENSION(2)  :: cdate
+   CHARACTER(LEN=10), DIMENSION(2) :: ctime
+   CHARACTER(LEN=5)                :: czone
+    
+   ! From of ouput file (1/proc or one global)   !RB to put in nammpp or namctl
+   LOGICAL :: ln_onefile = .TRUE. 
+   LOGICAL :: lwriter
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE timing_start(cdinfo)
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE timing_start  ***
+      !! ** Purpose :   collect execution time
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in) :: cdinfo
+      !
+       IF(ASSOCIATED(s_timer) ) s_timer_old => s_timer
+       !
+      ! Create timing structure at first call of the routine 
+       CALL timing_ini_var(cdinfo)
+   !   write(*,*) 'after inivar ', s_timer%cname
+
+      ! ici timing_ini_var a soit retrouve s_timer et fait return soit ajoute un maillon
+      ! maintenant on regarde si le call d'avant corrsspond a un parent ou si il est ferme
+      IF( .NOT. s_timer_old%l_tdone ) THEN      
+         s_timer%parent_section => s_timer_old
+      ELSE
+         s_timer%parent_section => NULL()
+      ENDIF    
+
+      s_timer%l_tdone = .FALSE.
+      s_timer%niter = s_timer%niter + 1
+      s_timer%t_cpu = 0.
+      s_timer%t_clock = 0.
+                  
+      ! CPU time collection
+      CALL CPU_TIME( s_timer%t_cpu  )
+      ! clock time collection
+#if defined key_mpp_mpi
+      s_timer%t_clock= MPI_Wtime()
+#else
+      CALL SYSTEM_CLOCK(COUNT_RATE=s_timer%ncount_rate, COUNT_MAX=s_timer%ncount_max)
+      CALL SYSTEM_CLOCK(COUNT = s_timer%ncount)
+#endif
+!      write(*,*) 'end of start ', s_timer%cname
+
+      !
+   END SUBROUTINE timing_start
+
+
+   SUBROUTINE timing_stop(cdinfo, csection)
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE timing_stop  ***
+      !! ** Purpose :   finalize timing and output
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in) :: cdinfo
+      CHARACTER(len=*), INTENT(in), OPTIONAL :: csection
+      !
+      INTEGER  :: ifinal_count, iperiods    
+      REAL(wp) :: zcpu_end, zmpitime,zcpu_raw,zclock_raw
+      !
+      s_wrk => NULL()
+
+      ! clock time collection
+#if defined key_mpp_mpi
+      zmpitime = MPI_Wtime()
+#else
+      CALL SYSTEM_CLOCK(COUNT = ifinal_count)
+#endif
+      ! CPU time collection
+      CALL CPU_TIME( zcpu_end )
+
+!!$      IF(associated(s_timer%parent_section))then
+!!$        write(*,*) s_timer%cname,' <-- ', s_timer%parent_section%cname
+!!$      ENDIF  
+
+ !     No need to search ... : s_timer has the last value defined in start
+ !     s_timer => s_timer_root
+ !     DO WHILE( TRIM(s_timer%cname) /= TRIM(cdinfo) ) 
+ !        IF( ASSOCIATED(s_timer%next) ) s_timer => s_timer%next
+ !     END DO
+ 
+      ! CPU time correction
+      zcpu_raw = zcpu_end - s_timer%t_cpu - t_overcpu ! total time including child
+      s_timer%t_cpu  = zcpu_raw - s_timer%tsub_cpu
+  !    IF(s_timer%cname==trim('lbc_lnk_2d'))  write(*,*) s_timer%tsub_cpu,zcpu_end
+
+      ! clock time correction
+#if defined key_mpp_mpi
+      zclock_raw = zmpitime - s_timer%t_clock - t_overclock ! total time including child
+      s_timer%t_clock = zclock_raw - t_overclock - s_timer%tsub_clock
+#else
+      iperiods = ifinal_count - s_timer%ncount
+      IF( ifinal_count < s_timer%ncount )  &
+         iperiods = iperiods + s_timer%ncount_max 
+         zclock_raw = REAL(iperiods) / s_timer%ncount_rate !- t_overclock   
+         s_timer%t_clock  = zclock_raw - s_timer%tsub_clock
+#endif
+ !     IF(s_timer%cname==trim('lbc_lnk_2d')) write(*,*) zclock_raw , s_timer%tsub_clock
+      
+      ! Correction of parent section
+      IF( .NOT. PRESENT(csection) ) THEN
+         IF ( ASSOCIATED(s_timer%parent_section ) ) THEN
+            s_timer%parent_section%tsub_cpu   = zcpu_raw   + s_timer%parent_section%tsub_cpu 
+            s_timer%parent_section%tsub_clock = zclock_raw + s_timer%parent_section%tsub_clock             
+         ENDIF
+      ENDIF
+            
+      ! time diagnostics 
+      s_timer%tsum_clock = s_timer%tsum_clock + s_timer%t_clock 
+      s_timer%tsum_cpu   = s_timer%tsum_cpu   + s_timer%t_cpu
+!RB to use to get min/max during a time integration
+!      IF( .NOT. l_initdone ) THEN
+!         s_timer%tmin_clock = s_timer%t_clock 
+!         s_timer%tmin_cpu   = s_timer%t_cpu 
+!      ELSE
+!         s_timer%tmin_clock = MIN( s_timer%tmin_clock, s_timer%t_clock ) 
+!         s_timer%tmin_cpu   = MIN( s_timer%tmin_cpu  , s_timer%t_cpu   ) 
+!      ENDIF   
+!      s_timer%tmax_clock = MAX( s_timer%tmax_clock, s_timer%t_clock ) 
+!      s_timer%tmax_cpu   = MAX( s_timer%tmax_cpu  , s_timer%t_cpu   )  
+      !
+      s_timer%tsub_clock = 0.
+      s_timer%tsub_cpu = 0.
+      s_timer%l_tdone = .TRUE.
+      !
+      !
+      ! we come back
+      IF ( ASSOCIATED(s_timer%parent_section ) ) s_timer => s_timer%parent_section
+     
+!      write(*,*) 'end of stop ', s_timer%cname
+
+   END SUBROUTINE timing_stop
+ 
+ 
+   SUBROUTINE timing_init
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE timing_init  ***
+      !! ** Purpose :   open timing output file
+      !!----------------------------------------------------------------------
+      INTEGER            :: iperiods, istart_count, ifinal_count
+      REAL(wp)           :: zdum
+      LOGICAL            :: ll_f
+      CHARACTER(len=128) :: ctiming_report_name  ! name of timing report file
+             
+      ctiming_report_name="timing.output"
+      IF( TRIM(cxios_context) == 'sas' ) WRITE(ctiming_report_name,'(A4,A13)') 'sas_', TRIM(ctiming_report_name)
+
+      IF( ln_onefile ) THEN
+         IF( lwp) CALL ctl_opn( numtime, TRIM(ctiming_report_name), 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,.TRUE., narea )
+         lwriter = lwp
+      ELSE
+         CALL ctl_opn( numtime, TRIM(ctiming_report_name), 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,.FALSE., narea )
+         lwriter = .TRUE.
+      ENDIF
+      
+      IF( lwriter) THEN      
+         WRITE(numtime,*)
+         WRITE(numtime,*) '      CNRS - NERC - Met OFFICE - MERCATOR-ocean - CMCC - INGV'
+         WRITE(numtime,*) '                             NEMO team'
+         WRITE(numtime,*) '                  Ocean General Circulation Model'
+         WRITE(numtime,*) '                        version 4.0  (2019) '
+         WRITE(numtime,*)
+         WRITE(numtime,*) '                        Timing Informations '
+         WRITE(numtime,*)
+         WRITE(numtime,*)
+      ENDIF   
+      
+      ! Compute clock function overhead
+#if defined key_mpp_mpi        
+      t_overclock = MPI_WTIME()
+      t_overclock = MPI_WTIME() - t_overclock
+#else        
+      CALL SYSTEM_CLOCK(COUNT_RATE=ncount_rate, COUNT_MAX=ncount_max)
+      CALL SYSTEM_CLOCK(COUNT = istart_count)
+      CALL SYSTEM_CLOCK(COUNT = ifinal_count)
+      iperiods = ifinal_count - istart_count
+      IF( ifinal_count < istart_count )  &
+          iperiods = iperiods + ncount_max 
+      t_overclock = REAL(iperiods) / ncount_rate
+#endif
+
+      ! Compute cpu_time function overhead
+      CALL CPU_TIME(zdum)
+      CALL CPU_TIME(t_overcpu)
+      
+      ! End overhead omputation  
+      t_overcpu = t_overcpu - zdum        
+      t_overclock = t_overcpu + t_overclock        
+
+      ! Timing on date and time
+      CALL DATE_AND_TIME(cdate(1),ctime(1),czone,nvalues)
+    
+      CALL CPU_TIME(t_cpu(1))      
+#if defined key_mpp_mpi        
+      ! Start elapsed and CPU time counters
+      t_elaps(1) = MPI_WTIME()
+#else
+      CALL SYSTEM_CLOCK(COUNT_RATE=ncount_rate, COUNT_MAX=ncount_max)
+      CALL SYSTEM_CLOCK(COUNT = ncount)
+#endif                 
+      !
+   END SUBROUTINE timing_init
+
+
+   SUBROUTINE timing_finalize
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE timing_finalize ***
+      !! ** Purpose :  compute average time 
+      !!               write timing output file
+      !!----------------------------------------------------------------------
+      TYPE(timer), POINTER :: s_temp
+      INTEGER :: idum, iperiods, icode
+      INTEGER :: ji
+      LOGICAL :: ll_ord, ll_averep
+      CHARACTER(len=120) :: clfmt            
+      REAL(wp), DIMENSION(:), ALLOCATABLE ::   timing_glob
+      REAL(wp) ::   zsypd   ! simulated years per day (Balaji 2017)
+      REAL(wp) ::   zperc, ztot
+
+      ll_averep = .TRUE.
+    
+      ! total CPU and elapse
+      CALL CPU_TIME(t_cpu(2))
+      t_cpu(2)   = t_cpu(2)    - t_cpu(1)   - t_overcpu
+#if defined key_mpp_mpi
+      t_elaps(2) = MPI_WTIME() - t_elaps(1) - t_overclock
+#else
+      CALL SYSTEM_CLOCK(COUNT = nfinal_count)
+      iperiods = nfinal_count - ncount
+      IF( nfinal_count < ncount )  &
+          iperiods = iperiods + ncount_max 
+      t_elaps(2) = REAL(iperiods) / ncount_rate - t_overclock
+#endif      
+
+      ! End of timings on date & time
+      CALL DATE_AND_TIME(cdate(2),ctime(2),czone,nvalues)
+       
+      ! Compute the numer of routines
+      nsize = 0 
+      s_timer => s_timer_root
+      DO WHILE( ASSOCIATED(s_timer) )
+         nsize = nsize + 1
+         s_timer => s_timer%next
+      END DO
+      idum = nsize
+      CALL mpp_sum('timing', idum)
+      IF( idum/jpnij /= nsize ) THEN
+         IF( lwriter ) WRITE(numtime,*) '        ===> W A R N I N G: '
+         IF( lwriter ) WRITE(numtime,*) ' Some CPU have different number of routines instrumented for timing'
+         IF( lwriter ) WRITE(numtime,*) ' No detailed report on averaged timing can be provided'
+         IF( lwriter ) WRITE(numtime,*) ' The following detailed report only deals with the current processor'
+         IF( lwriter ) WRITE(numtime,*)
+         ll_averep = .FALSE.
+      ENDIF   
+
+#if defined key_mpp_mpi      
+      ! in MPI gather some info
+      ALLOCATE( all_etime(jpnij), all_ctime(jpnij) )
+      CALL MPI_ALLGATHER(t_elaps(2), 1, MPI_DOUBLE_PRECISION,   &
+                         all_etime , 1, MPI_DOUBLE_PRECISION,   &
+                         MPI_COMM_OCE, icode)
+      CALL MPI_ALLGATHER(t_cpu(2) , 1, MPI_DOUBLE_PRECISION,   &
+                         all_ctime, 1, MPI_DOUBLE_PRECISION,   &
+                         MPI_COMM_OCE, icode)
+      tot_etime = SUM(all_etime(:))
+      tot_ctime = SUM(all_ctime(:))
+#else
+      tot_etime = t_elaps(2)
+      tot_ctime = t_cpu  (2)           
+#endif
+
+      ! write output file
+      IF( lwriter ) WRITE(numtime,*) 'Total timing (sum) :'
+      IF( lwriter ) WRITE(numtime,*) '--------------------'
+      IF( lwriter ) WRITE(numtime,"('Elapsed Time (s)  CPU Time (s)')")
+      IF( lwriter ) WRITE(numtime,'(5x,f12.3,1x,f12.3)')  tot_etime, tot_ctime
+      IF( lwriter ) WRITE(numtime,*) 
+#if defined key_mpp_mpi
+      IF( ll_averep ) CALL waver_info
+      CALL wmpi_info
+#endif      
+      IF( lwriter ) CALL wcurrent_info
+      
+      clfmt='(1X,"Timing started on ",2(A2,"/"),A4," at ",2(A2,":"),A2," MET ",A3,":",A2," from GMT")'
+      IF( lwriter ) WRITE(numtime, TRIM(clfmt)) &           
+      &       cdate(1)(7:8), cdate(1)(5:6), cdate(1)(1:4),   &
+      &       ctime(1)(1:2), ctime(1)(3:4), ctime(1)(5:6),   &
+      &       czone(1:3),    czone(4:5)                     
+      clfmt='(1X,  "Timing   ended on ",2(A2,"/"),A4," at ",2(A2,":"),A2," MET ",A3,":",A2," from GMT")'
+      IF( lwriter ) WRITE(numtime, TRIM(clfmt)) &           
+      &       cdate(2)(7:8), cdate(2)(5:6), cdate(2)(1:4),   &
+      &       ctime(2)(1:2), ctime(2)(3:4), ctime(2)(5:6),   &
+      &       czone(1:3),    czone(4:5)
+
+#if defined key_mpp_mpi
+      ALLOCATE(timing_glob(4*jpnij), stat=icode)
+      CALL MPI_GATHER( (/compute_time, waiting_time(1), waiting_time(2), elapsed_time/),   &
+         &             4, MPI_DOUBLE_PRECISION, timing_glob, 4, MPI_DOUBLE_PRECISION, 0, MPI_COMM_OCE, icode)
+      IF( narea == 1 ) THEN
+         WRITE(numtime,*) ' '
+         WRITE(numtime,*) ' Report on time spent on waiting MPI messages '
+         WRITE(numtime,*) '    total timing measured between nit000+1 and nitend-1 '
+         WRITE(numtime,*) '    warning: includes restarts writing time if output before nitend... '
+         WRITE(numtime,*) ' '
+         DO ji = 1, jpnij
+            ztot = SUM( timing_glob(4*ji-3:4*ji-1) )
+            WRITE(numtime,'(A28,F11.6,            A34,I8)') 'Computing       time : ',timing_glob(4*ji-3), ' on MPI rank : ', ji
+            IF ( ztot /= 0. ) zperc = timing_glob(4*ji-2) / ztot * 100.
+            WRITE(numtime,'(A28,F11.6,A2, F4.1,A3,A25,I8)') 'Waiting lbc_lnk time : ',timing_glob(4*ji-2)   &
+               &                                                         , ' (',      zperc,' %)',   ' on MPI rank : ', ji
+            IF ( ztot /= 0. ) zperc = timing_glob(4*ji-1) / ztot * 100.
+            WRITE(numtime,'(A28,F11.6,A2, F4.1,A3,A25,I8)') 'Waiting  global time : ',timing_glob(4*ji-1)   &
+               &                                                         , ' (',      zperc,' %)',   ' on MPI rank : ', ji
+            zsypd = rn_rdt * REAL(nitend-nit000-1, wp) / (timing_glob(4*ji) * 365.)
+            WRITE(numtime,'(A28,F11.6,A7,F10.3,A2,A15,I8)') 'Total           time : ',timing_glob(4*ji  )   &
+               &                                                         , ' (SYPD: ', zsypd, ')',   ' on MPI rank : ', ji
+         END DO
+      ENDIF
+      DEALLOCATE(timing_glob)
+#endif      
+
+      IF( lwriter ) CLOSE(numtime) 
+      !
+   END SUBROUTINE timing_finalize
+   
+
+   SUBROUTINE wcurrent_info
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE wcurrent_info ***
+      !! ** Purpose :  compute and write timing output file
+      !!----------------------------------------------------------------------
+      LOGICAL :: ll_ord
+      CHARACTER(len=2048) :: clfmt            
+   
+      ! reorder the current list by elapse time      
+      s_wrk => NULL()
+      s_timer => s_timer_root
+      DO
+         ll_ord = .TRUE.
+         s_timer => s_timer_root
+         DO WHILE ( ASSOCIATED( s_timer%next ) )
+         IF (.NOT. ASSOCIATED(s_timer%next)) EXIT
+            IF ( s_timer%tsum_clock < s_timer%next%tsum_clock ) THEN 
+               ALLOCATE(s_wrk)
+               s_wrk = s_timer%next
+               CALL insert  (s_timer, s_timer_root, s_wrk)
+               CALL suppress(s_timer%next)            
+               ll_ord = .FALSE.
+               CYCLE            
+            ENDIF           
+         IF( ASSOCIATED(s_timer%next) ) s_timer => s_timer%next
+         END DO         
+         IF( ll_ord ) EXIT
+      END DO
+            
+      ! write current info
+      WRITE(numtime,*) 'Detailed timing for proc :', narea-1
+      WRITE(numtime,*) '--------------------------'
+      WRITE(numtime,*) 'Section             ',            &
+      &   'Elapsed Time (s)  ','Elapsed Time (%)  ',   &
+      &   'CPU Time(s)  ','CPU Time (%)  ','CPU/Elapsed  ','Frequency' 
+      s_timer => s_timer_root  
+      clfmt = '(1x,a,4x,f12.3,6x,f12.3,x,f12.3,2x,f12.3,6x,f7.3,2x,i9)'
+      DO WHILE ( ASSOCIATED(s_timer) )
+         WRITE(numtime,TRIM(clfmt))   s_timer%cname,   &
+         &   s_timer%tsum_clock,s_timer%tsum_clock*100./t_elaps(2),            &
+         &   s_timer%tsum_cpu  ,s_timer%tsum_cpu*100./t_cpu(2)    ,            &
+         &   s_timer%tsum_cpu/s_timer%tsum_clock, s_timer%niter
+         s_timer => s_timer%next
+      END DO
+      WRITE(numtime,*)
+      !                  
+   END SUBROUTINE wcurrent_info
+
+#if defined key_mpp_mpi     
+   SUBROUTINE waver_info
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE wcurrent_info ***
+      !! ** Purpose :  compute and write averaged timing informations
+      !!----------------------------------------------------------------------
+      TYPE(alltimer), POINTER :: sl_timer_glob_root => NULL()
+      TYPE(alltimer), POINTER :: sl_timer_glob      => NULL()
+      TYPE(timer), POINTER :: sl_timer_ave_root => NULL()
+      TYPE(timer), POINTER :: sl_timer_ave      => NULL()
+      INTEGER :: icode
+      INTEGER :: ierr
+      LOGICAL :: ll_ord           
+      CHARACTER(len=200) :: clfmt              
+                 
+      ! Initialised the global strucutre   
+      ALLOCATE(sl_timer_glob_root, Stat=ierr)
+      IF(ierr /= 0)THEN
+         WRITE(numtime,*) 'Failed to allocate global timing structure in waver_info'
+         RETURN
+      END IF
+
+      ALLOCATE(sl_timer_glob_root%cname     (jpnij), &
+               sl_timer_glob_root%tsum_cpu  (jpnij), &
+               sl_timer_glob_root%tsum_clock(jpnij), &
+               sl_timer_glob_root%niter     (jpnij), Stat=ierr)
+      IF(ierr /= 0)THEN
+         WRITE(numtime,*) 'Failed to allocate global timing structure in waver_info'
+         RETURN
+      END IF
+      sl_timer_glob_root%cname(:)       = ''
+      sl_timer_glob_root%tsum_cpu(:)   = 0._wp
+      sl_timer_glob_root%tsum_clock(:) = 0._wp
+      sl_timer_glob_root%niter(:)      = 0
+      sl_timer_glob_root%next => NULL()
+      sl_timer_glob_root%prev => NULL()
+      !ARPDBG - don't need to allocate a pointer that's immediately then
+      !         set to point to some other object.
+      !ALLOCATE(sl_timer_glob)
+      !ALLOCATE(sl_timer_glob%cname     (jpnij))
+      !ALLOCATE(sl_timer_glob%tsum_cpu  (jpnij))
+      !ALLOCATE(sl_timer_glob%tsum_clock(jpnij))
+      !ALLOCATE(sl_timer_glob%niter     (jpnij))
+      sl_timer_glob => sl_timer_glob_root
+      !
+      IF( narea .EQ. 1 ) THEN
+         ALLOCATE(sl_timer_ave_root)
+         sl_timer_ave_root%cname       = ''
+         sl_timer_ave_root%t_cpu      = 0._wp
+         sl_timer_ave_root%t_clock    = 0._wp
+         sl_timer_ave_root%tsum_cpu   = 0._wp
+         sl_timer_ave_root%tsum_clock = 0._wp
+         sl_timer_ave_root%tmax_cpu   = 0._wp
+         sl_timer_ave_root%tmax_clock = 0._wp
+         sl_timer_ave_root%tmin_cpu   = 0._wp
+         sl_timer_ave_root%tmin_clock = 0._wp
+         sl_timer_ave_root%tsub_cpu   = 0._wp
+         sl_timer_ave_root%tsub_clock = 0._wp
+         sl_timer_ave_root%ncount      = 0
+         sl_timer_ave_root%ncount_rate = 0
+         sl_timer_ave_root%ncount_max  = 0
+         sl_timer_ave_root%niter       = 0
+         sl_timer_ave_root%l_tdone  = .FALSE.
+         sl_timer_ave_root%next => NULL()
+         sl_timer_ave_root%prev => NULL()
+         ALLOCATE(sl_timer_ave)
+         sl_timer_ave => sl_timer_ave_root            
+      ENDIF 
+
+      ! Gather info from all processors
+      s_timer => s_timer_root
+      DO WHILE ( ASSOCIATED(s_timer) )
+         CALL MPI_GATHER(s_timer%cname     , 20, MPI_CHARACTER,   &
+                         sl_timer_glob%cname, 20, MPI_CHARACTER,   &
+                         0, MPI_COMM_OCE, icode)
+         CALL MPI_GATHER(s_timer%tsum_clock     , 1, MPI_DOUBLE_PRECISION,   &
+                         sl_timer_glob%tsum_clock, 1, MPI_DOUBLE_PRECISION,   &
+                         0, MPI_COMM_OCE, icode)
+         CALL MPI_GATHER(s_timer%tsum_cpu     , 1, MPI_DOUBLE_PRECISION,   &
+                         sl_timer_glob%tsum_cpu, 1, MPI_DOUBLE_PRECISION,   &
+                         0, MPI_COMM_OCE, icode)
+         CALL MPI_GATHER(s_timer%niter     , 1, MPI_INTEGER,   &
+                         sl_timer_glob%niter, 1, MPI_INTEGER,   &
+                         0, MPI_COMM_OCE, icode)
+
+         IF( narea == 1 .AND. ASSOCIATED(s_timer%next) ) THEN
+            ALLOCATE(sl_timer_glob%next)
+            ALLOCATE(sl_timer_glob%next%cname     (jpnij))
+            ALLOCATE(sl_timer_glob%next%tsum_cpu  (jpnij))
+            ALLOCATE(sl_timer_glob%next%tsum_clock(jpnij))
+            ALLOCATE(sl_timer_glob%next%niter     (jpnij))
+            sl_timer_glob%next%prev => sl_timer_glob
+            sl_timer_glob%next%next => NULL()
+            sl_timer_glob           => sl_timer_glob%next
+         ENDIF              
+         s_timer => s_timer%next
+      END DO      
+      
+      IF( narea == 1 ) THEN    
+         ! Compute some stats
+         sl_timer_glob => sl_timer_glob_root
+         DO WHILE( ASSOCIATED(sl_timer_glob) )
+            sl_timer_ave%cname  = sl_timer_glob%cname(1)
+            sl_timer_ave%tsum_cpu   = SUM   (sl_timer_glob%tsum_cpu  (:)) / jpnij
+            sl_timer_ave%tsum_clock = SUM   (sl_timer_glob%tsum_clock(:)) / jpnij
+            sl_timer_ave%tmax_cpu   = MAXVAL(sl_timer_glob%tsum_cpu  (:))
+            sl_timer_ave%tmax_clock = MAXVAL(sl_timer_glob%tsum_clock(:))
+            sl_timer_ave%tmin_cpu   = MINVAL(sl_timer_glob%tsum_cpu  (:))
+            sl_timer_ave%tmin_clock = MINVAL(sl_timer_glob%tsum_clock(:))
+            sl_timer_ave%niter      = SUM   (sl_timer_glob%niter     (:))
+            !
+            IF( ASSOCIATED(sl_timer_glob%next) ) THEN
+               ALLOCATE(sl_timer_ave%next)          
+               sl_timer_ave%next%prev => sl_timer_ave
+               sl_timer_ave%next%next => NULL()           
+               sl_timer_ave           => sl_timer_ave%next
+            ENDIF
+            sl_timer_glob => sl_timer_glob%next                                
+         END DO
+      
+         ! reorder the averaged list by CPU time      
+         s_wrk => NULL()
+         sl_timer_ave => sl_timer_ave_root
+         DO
+            ll_ord = .TRUE.
+            sl_timer_ave => sl_timer_ave_root
+            DO WHILE( ASSOCIATED( sl_timer_ave%next ) )
+
+               IF( .NOT. ASSOCIATED(sl_timer_ave%next) ) EXIT
+
+               IF ( sl_timer_ave%tsum_clock < sl_timer_ave%next%tsum_clock ) THEN 
+                  ALLOCATE(s_wrk)
+                  ! Copy data into the new object pointed to by s_wrk
+                  s_wrk = sl_timer_ave%next
+                  ! Insert this new timer object before our current position
+                  CALL insert  (sl_timer_ave, sl_timer_ave_root, s_wrk)
+                  ! Remove the old object from the list
+                  CALL suppress(sl_timer_ave%next)            
+                  ll_ord = .FALSE.
+                  CYCLE            
+               ENDIF           
+               IF( ASSOCIATED(sl_timer_ave%next) ) sl_timer_ave => sl_timer_ave%next
+            END DO         
+            IF( ll_ord ) EXIT
+         END DO
+
+         ! write averaged info
+         WRITE(numtime,"('Averaged timing on all processors :')")
+         WRITE(numtime,"('-----------------------------------')")
+         WRITE(numtime,"('Section',13x,'Elap. Time(s)',2x,'Elap. Time(%)',2x, &
+         &   'CPU Time(s)',2x,'CPU Time(%)',2x,'CPU/Elap',1x,   &
+         &   'Max elap(%)',2x,'Min elap(%)',2x,            &           
+         &   'Freq')")
+         sl_timer_ave => sl_timer_ave_root  
+         clfmt = '((A),E15.7,2x,f6.2,5x,f12.2,5x,f6.2,5x,f7.2,2x,f12.2,4x,f6.2,2x,f9.2)'
+         DO WHILE ( ASSOCIATED(sl_timer_ave) )
+            WRITE(numtime,TRIM(clfmt))   sl_timer_ave%cname(1:18),                            &
+            &   sl_timer_ave%tsum_clock,sl_timer_ave%tsum_clock*100.*jpnij/tot_etime,   &
+            &   sl_timer_ave%tsum_cpu  ,sl_timer_ave%tsum_cpu*100.*jpnij/tot_ctime  ,   &
+            &   sl_timer_ave%tsum_cpu/sl_timer_ave%tsum_clock,                          &
+            &   sl_timer_ave%tmax_clock*100.*jpnij/tot_etime,                           &
+            &   sl_timer_ave%tmin_clock*100.*jpnij/tot_etime,                           &                                               
+            &   sl_timer_ave%niter/REAL(jpnij)
+            sl_timer_ave => sl_timer_ave%next
+         END DO
+         WRITE(numtime,*)
+         !
+         DEALLOCATE(sl_timer_ave_root)
+      ENDIF
+      !
+      DEALLOCATE(sl_timer_glob_root)
+      !                  
+   END SUBROUTINE waver_info
+  
+  
+   SUBROUTINE wmpi_info
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE wmpi_time  ***
+      !! ** Purpose :   compute and write a summary of MPI infos 
+      !!----------------------------------------------------------------------   
+      !   
+      INTEGER                            :: idum, icode
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: iall_rank
+      REAL(wp) :: ztot_ratio
+      REAL(wp) :: zmax_etime, zmax_ctime, zmax_ratio, zmin_etime, zmin_ctime, zmin_ratio
+      REAL(wp) :: zavg_etime, zavg_ctime, zavg_ratio
+      REAL(wp), ALLOCATABLE, DIMENSION(:) :: zall_ratio
+      CHARACTER(LEN=128), dimension(8) :: cllignes
+      CHARACTER(LEN=128)               :: clhline, clstart_date, clfinal_date
+      CHARACTER(LEN=2048)              :: clfmt    
+   
+      ! Gather all times
+      ALLOCATE( zall_ratio(jpnij), iall_rank(jpnij) )
+      IF( narea == 1 ) THEN
+         iall_rank(:) = (/ (idum,idum=0,jpnij-1) /)
+   
+         ! Compute elapse user time
+         zavg_etime = tot_etime/REAL(jpnij,wp)
+         zmax_etime = MAXVAL(all_etime(:))
+         zmin_etime = MINVAL(all_etime(:))
+
+         ! Compute CPU user time
+         zavg_ctime = tot_ctime/REAL(jpnij,wp)
+         zmax_ctime = MAXVAL(all_ctime(:))
+         zmin_ctime = MINVAL(all_ctime(:))
+   
+         ! Compute cpu/elapsed ratio
+         zall_ratio(:) = all_ctime(:) / all_etime(:)
+         ztot_ratio    = SUM(zall_ratio(:))
+         zavg_ratio    = ztot_ratio/REAL(jpnij,wp)
+         zmax_ratio    = MAXVAL(zall_ratio(:))
+         zmin_ratio    = MINVAL(zall_ratio(:))   
+   
+         ! Output Format
+         clhline    ='1x,13("-"),"|",18("-"),"|",14("-"),"|",18("-"),/,'
+         cllignes(1)='(1x,"MPI summary report :",/,'
+         cllignes(2)='1x,"--------------------",//,'
+         cllignes(3)='1x,"Process Rank |"," Elapsed Time (s) |"," CPU Time (s) |"," Ratio CPU/Elapsed",/,'
+         cllignes(4)='    (1x,i4,9x,"|",f12.3,6x,"|",f12.3,2x,"|",4x,f7.3,/),'
+         WRITE(cllignes(4)(1:4),'(I4)') jpnij
+         cllignes(5)='1x,"Total        |",f12.3,6x,"|",F12.3,2x,"|",4x,f7.3,/,'
+         cllignes(6)='1x,"Minimum      |",f12.3,6x,"|",F12.3,2x,"|",4x,f7.3,/,'
+         cllignes(7)='1x,"Maximum      |",f12.3,6x,"|",F12.3,2x,"|",4x,f7.3,/,'
+         cllignes(8)='1x,"Average      |",f12.3,6x,"|",F12.3,2x,"|",4x,f7.3)'
+         clfmt=TRIM(cllignes(1))// TRIM(cllignes(2))//TRIM(cllignes(3))//          &
+           & TRIM(clhline)//TRIM(cllignes(4))//TRIM(clhline)//TRIM(cllignes(5))//  &
+           & TRIM(clhline)//TRIM(cllignes(6))//TRIM(clhline)//TRIM(cllignes(7))//  &
+           & TRIM(clhline)//TRIM(cllignes(8))
+         WRITE(numtime, TRIM(clfmt)) &
+             (iall_rank(idum),all_etime(idum),all_ctime(idum),zall_ratio(idum),idum=1, jpnij), &
+             tot_etime,     tot_ctime,     ztot_ratio,   &
+             zmin_etime,    zmin_ctime,    zmin_ratio,   &
+             zmax_etime,    zmax_ctime,    zmax_ratio,   &
+             zavg_etime,    zavg_ctime,    zavg_ratio
+         WRITE(numtime,*)    
+      END IF
+      !
+      DEALLOCATE(zall_ratio, iall_rank)
+      !
+   END SUBROUTINE wmpi_info
+#endif   
+
+
+   SUBROUTINE timing_ini_var(cdinfo)
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE timing_ini_var  ***
+      !! ** Purpose :   create timing structure 
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in) :: cdinfo
+      LOGICAL :: ll_section
+       
+      !
+      IF( .NOT. ASSOCIATED(s_timer_root) ) THEN
+         ALLOCATE(s_timer_root)
+         s_timer_root%cname       = cdinfo
+         s_timer_root%t_cpu      = 0._wp
+         s_timer_root%t_clock    = 0._wp
+         s_timer_root%tsum_cpu   = 0._wp
+         s_timer_root%tsum_clock = 0._wp
+         s_timer_root%tmax_cpu   = 0._wp
+         s_timer_root%tmax_clock = 0._wp
+         s_timer_root%tmin_cpu   = 0._wp
+         s_timer_root%tmin_clock = 0._wp
+         s_timer_root%tsub_cpu   = 0._wp
+         s_timer_root%tsub_clock = 0._wp
+         s_timer_root%ncount      = 0
+         s_timer_root%ncount_rate = 0
+         s_timer_root%ncount_max  = 0
+         s_timer_root%niter       = 0
+         s_timer_root%l_tdone  = .FALSE.
+         s_timer_root%next => NULL()
+         s_timer_root%prev => NULL()
+         s_timer => s_timer_root
+         !
+         ALLOCATE(s_wrk)
+         s_wrk => NULL()
+         !
+         ALLOCATE(s_timer_old)
+         s_timer_old%cname       = cdinfo
+         s_timer_old%t_cpu      = 0._wp
+         s_timer_old%t_clock    = 0._wp
+         s_timer_old%tsum_cpu   = 0._wp
+         s_timer_old%tsum_clock = 0._wp
+         s_timer_old%tmax_cpu   = 0._wp
+         s_timer_old%tmax_clock = 0._wp
+         s_timer_old%tmin_cpu   = 0._wp
+         s_timer_old%tmin_clock = 0._wp
+         s_timer_old%tsub_cpu   = 0._wp
+         s_timer_old%tsub_clock = 0._wp
+         s_timer_old%ncount      = 0
+         s_timer_old%ncount_rate = 0
+         s_timer_old%ncount_max  = 0
+         s_timer_old%niter       = 0
+         s_timer_old%l_tdone  = .TRUE.
+         s_timer_old%next => NULL()
+         s_timer_old%prev => NULL()
+
+      ELSE
+         s_timer => s_timer_root
+         ! case of already existing area (typically inside a loop)
+   !         write(*,*) 'in ini_var for routine : ', cdinfo
+         DO WHILE( ASSOCIATED(s_timer) ) 
+            IF( TRIM(s_timer%cname) .EQ. TRIM(cdinfo) ) THEN
+ !             write(*,*) 'in ini_var for routine : ', cdinfo,' we return'           
+               RETURN ! cdinfo is already in the chain
+            ENDIF
+            s_timer => s_timer%next
+         END DO
+
+         ! end of the chain
+         s_timer => s_timer_root
+         DO WHILE( ASSOCIATED(s_timer%next) )
+            s_timer => s_timer%next
+         END DO
+
+    !     write(*,*) 'after search', s_timer%cname
+         ! cdinfo is not part of the chain so we add it with initialisation          
+          ALLOCATE(s_timer%next)
+    !     write(*,*) 'after allocation of next'
+  
+         s_timer%next%cname       = cdinfo
+         s_timer%next%t_cpu      = 0._wp
+         s_timer%next%t_clock    = 0._wp
+         s_timer%next%tsum_cpu   = 0._wp
+         s_timer%next%tsum_clock = 0._wp  
+         s_timer%next%tmax_cpu   = 0._wp
+         s_timer%next%tmax_clock = 0._wp
+         s_timer%next%tmin_cpu   = 0._wp
+         s_timer%next%tmin_clock = 0._wp
+         s_timer%next%tsub_cpu   = 0._wp
+         s_timer%next%tsub_clock = 0._wp
+         s_timer%next%ncount      = 0
+         s_timer%next%ncount_rate = 0
+         s_timer%next%ncount_max  = 0
+         s_timer%next%niter       = 0
+         s_timer%next%l_tdone  = .FALSE.
+         s_timer%next%parent_section => NULL()
+         s_timer%next%prev => s_timer
+         s_timer%next%next => NULL()
+         s_timer => s_timer%next
+      ENDIF 
+      !    write(*,*) 'after allocation'
+     !
+   END SUBROUTINE timing_ini_var
+
+
+   SUBROUTINE timing_reset
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE timing_reset  ***
+      !! ** Purpose :   go to root of timing tree 
+      !!----------------------------------------------------------------------
+      l_initdone = .TRUE. 
+!      IF(lwp) WRITE(numout,*)
+!      IF(lwp) WRITE(numout,*) 'timing_reset : instrumented routines for timing'
+!      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      CALL timing_list(s_timer_root)
+!      WRITE(numout,*)
+      !
+   END SUBROUTINE timing_reset
+
+
+   RECURSIVE SUBROUTINE timing_list(ptr)
+   
+      TYPE(timer), POINTER, INTENT(inout) :: ptr
+      !
+      IF( ASSOCIATED(ptr%next) ) CALL timing_list(ptr%next)
+      IF(lwp) WRITE(numout,*)'   ', ptr%cname   
+      !
+   END SUBROUTINE timing_list
+
+
+   SUBROUTINE insert(sd_current, sd_root ,sd_ptr)
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE insert  ***
+      !! ** Purpose :   insert an element in timer structure
+      !!----------------------------------------------------------------------
+      TYPE(timer), POINTER, INTENT(inout) :: sd_current, sd_root, sd_ptr
+      !
+     
+      IF( ASSOCIATED( sd_current, sd_root ) ) THEN
+         ! If our current element is the root element then
+         ! replace it with the one being inserted
+         sd_root => sd_ptr
+      ELSE
+         sd_current%prev%next => sd_ptr
+      END IF
+      sd_ptr%next     => sd_current
+      sd_ptr%prev     => sd_current%prev
+      sd_current%prev => sd_ptr
+      ! Nullify the pointer to the new element now that it is held
+      ! within the list. If we don't do this then a subsequent call
+      ! to ALLOCATE memory to this pointer will fail.
+      sd_ptr => NULL()
+      !    
+   END SUBROUTINE insert
+  
+  
+   SUBROUTINE suppress(sd_ptr)
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE suppress  ***
+      !! ** Purpose :   supress an element in timer structure
+      !!----------------------------------------------------------------------
+      TYPE(timer), POINTER, INTENT(inout) :: sd_ptr
+      !
+      TYPE(timer), POINTER :: sl_temp
+    
+      sl_temp => sd_ptr
+      sd_ptr => sd_ptr%next    
+      IF ( ASSOCIATED(sl_temp%next) ) sl_temp%next%prev => sl_temp%prev
+      DEALLOCATE(sl_temp)
+      sl_temp => NULL()
+      !
+    END SUBROUTINE suppress
+
+   !!=====================================================================
+END MODULE timing
