SUBROUTINE mpp_nfd_/**/PRECISION( ptab, cd_nat, psgn, kfillmode, pfillval, khls, kfld )
      TYPE(PTR_4d_/**/PRECISION),  DIMENSION(:), INTENT(inout) ::   ptab        ! pointer of arrays on which apply the b.c.
      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   cd_nat      ! nature of array grid-points
      REAL(wp),  DIMENSION(:), INTENT(in   ) ::   psgn        ! sign used across the north fold boundary
      INTEGER                       , INTENT(in   ) ::   kfillmode   ! filling method for halo over land 
      REAL(PRECISION)               , INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries)
      INTEGER                       , INTENT(in   ) ::   khls        ! halo size, default = nn_hls
      INTEGER                       , INTENT(in   ) ::   kfld        ! number of pt3d arrays
      !
      LOGICAL  ::   ll_add_line
      INTEGER  ::   ji,  jj,  jk,  jl, jh, jf, jr   ! dummy loop indices
      INTEGER  ::   ipi, ipj, ipj2, ipk, ipl, ipf   ! dimension of the input array
      INTEGER  ::   imigr, iihom, ijhom             ! local integers
      INTEGER  ::   ierr, ibuffsize, iis0, iie0, impp
      INTEGER  ::   ii1, ii2, ij1, ij2
      INTEGER  ::   ipimax, i0max
      INTEGER  ::   ij, iproc, ipni, ijnr
      INTEGER, DIMENSION (jpmaxngh)       ::   ml_req_nf   ! for mpi_isend when avoiding mpi_allgather
      INTEGER                             ::   ml_err      ! for mpi_isend when avoiding mpi_allgather
      !                                                    ! Workspace for message transfers avoiding mpi_allgather
      INTEGER                             ::   ipj_b       ! sum of lines for all multi fields
      INTEGER                             ::   i012        ! 0, 1 or 2
      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   jj_s  ! position of sent lines
      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   jj_b  ! position of buffer lines
      INTEGER , DIMENSION(:)          , ALLOCATABLE ::   ipj_s ! number of sent lines
      REAL(PRECISION), DIMENSION(:,:,:,:)    , ALLOCATABLE ::   ztabb, ztabr, ztabw  ! buffer, receive and work arrays
      REAL(PRECISION), DIMENSION(:,:,:,:,:)  , ALLOCATABLE ::   znorthloc
      REAL(PRECISION), DIMENSION(:,:,:,:,:,:), ALLOCATABLE ::   znorthglo
      TYPE(PTR_4D_/**/PRECISION), DIMENSION(:), ALLOCATABLE ::   ztabglo        ! array or pointer of arrays on which apply the b.c.
      !!----------------------------------------------------------------------
      !
      ipk = SIZE(ptab(1)%pt4d,3)
      ipl = SIZE(ptab(1)%pt4d,4)
      ipf = kfld
      !
      IF( ln_nnogather ) THEN      !==  no allgather exchanges  ==!

         !   ---   define number of exchanged lines   ---
         !
         ! In theory we should exchange only nn_hls lines.
         !
         ! However, some other points are duplicated in the north pole folding:
         !  - c_NFtype='T', grid=T : half of the last line (jpiglo/2+2:jpiglo-nn_hls)
         !  - c_NFtype='T', grid=U : half of the last line (jpiglo/2+1:jpiglo-nn_hls)
         !  - c_NFtype='T', grid=V : all the last line nn_hls+1 and (nn_hls+2:jpiglo-nn_hls)
         !  - c_NFtype='T', grid=F : all the last line (nn_hls+1:jpiglo-nn_hls)
         !  - c_NFtype='F', grid=T : 2 points of the last line (jpiglo/2+1 and jpglo-nn_hls)
         !  - c_NFtype='F', grid=U : no points are duplicated
         !  - c_NFtype='F', grid=V : half of the last line (jpiglo/2+1:jpiglo-nn_hls)
         !  - c_NFtype='F', grid=F : half of the last line (jpiglo/2+1:jpiglo-nn_hls-1)
         ! The order of the calculations may differ for these duplicated points (as, for example jj+1 becomes jj-1)
         ! This explain why these duplicated points may have different values even if they are at the exact same location.
         ! In consequence, we may want to force the folding on these points by setting l_full_nf_update = .TRUE.
         ! This is slightly slower but necessary to avoid different values on identical grid points!!
         !
         !!!!!!!!!           temporary switch off this optimisation ==> force TRUE           !!!!!!!!
         !!!!!!!!!  needed to get the same results without agrif and with agrif and no zoom  !!!!!!!!
         !!!!!!!!!                    I don't know why we must do that...                    !!!!!!!!
         l_full_nf_update = .TRUE.
         ! also force it if not restart during the first 2 steps (leap frog?)
         ll_add_line = l_full_nf_update .OR. ( ncom_stp <= nit000+1 .AND. .NOT. ln_rstart )
         
         ALLOCATE(ipj_s(ipf))                ! how many lines do we exchange?
         IF( ll_add_line ) THEN
            DO jf = 1, ipf                      ! Loop over the number of arrays to be processed
               ipj_s(jf) = khls + COUNT( (/ c_NFtype == 'T' .OR. cd_nat(jf) == 'V' .OR. cd_nat(jf) == 'F' /) ) 
            END DO
         ELSE
            ipj_s(:) = khls
         ENDIF
         
         ipj   = MAXVAL(ipj_s(:))            ! Max 2nd dimension of message transfers
         ipj_b = SUM(   ipj_s(:))            ! Total number of lines to be exchanged
         ALLOCATE( jj_s(ipj, ipf), jj_b(ipj, ipf) )

         ! Index of modifying lines in input
         ij1 = 0
         DO jf = 1, ipf                      ! Loop over the number of arrays to be processed
            !
            IF( c_NFtype == 'T' ) THEN          ! *  North fold  T-point pivot
               SELECT CASE ( cd_nat(jf) )
               CASE ( 'T', 'W', 'U' )   ;   i012 = 1   ! T-, U-, W-point
               CASE ( 'V', 'F'      )   ;   i012 = 2   ! V-, F-point
               END SELECT
            ENDIF
            IF( c_NFtype == 'F' ) THEN          ! *  North fold  F-point pivot
               SELECT CASE ( cd_nat(jf) )
               CASE ( 'T', 'W', 'U' )   ;   i012 = 0   ! T-, U-, W-point
               CASE ( 'V', 'F'      )   ;   i012 = 1   ! V-, F-point
               END SELECT
            ENDIF
               !
            DO jj = 1, ipj_s(jf)
               ij1 = ij1 + 1
               jj_b(jj,jf) = ij1
               jj_s(jj,jf) = jpj - 2*khls + jj - i012
            END DO
            !
         END DO
         !
         ALLOCATE( ztabb(jpimax,ipj_b,ipk,ipl) )   ! store all the data to be sent in a buffer array
         ibuffsize = jpimax * ipj_b * ipk * ipl
         !
         DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
            DO jj = 1, ipj_s(jf)
               ij1 = jj_b(jj,jf)
               ij2 = jj_s(jj,jf)
               DO ji = 1, jpi
                  ztabb(ji,ij1,jk,jl) = ptab(jf)%pt4d(ji,ij2,jk,jl)
               END DO
               DO ji = jpi+1, jpimax
                  ztabb(ji,ij1,jk,jl) = HUGE(0._/**/PRECISION)   ! avoid sending uninitialized values (make sure we don't use it)
               END DO
            END DO
         END DO   ;   END DO   ;   END DO
         !
         ! start waiting time measurement
         IF( ln_timing ) CALL tic_tac(.TRUE.)
         !
         ! send the data as soon as possible
         DO jr = 1, nsndto
            iproc = nfproc(isendto(jr))
            IF( iproc /= narea-1 .AND. iproc /= -1 ) THEN
#if ! defined key_mpi_off
               CALL MPI_ISEND( ztabb, ibuffsize, MPI_TYPE, iproc, 5, mpi_comm_oce, ml_req_nf(jr), ierr )
#endif
            ENDIF
         END DO
         !
         ipimax = jpimax * jpmaxngh
         ALLOCATE( ztabw(jpimax,ipj_b,ipk,ipl), ztabr(ipimax,ipj_b,ipk,ipl) ) 
         !
         DO jr = 1, nsndto
            !
            ipni  = isendto(jr)
            iproc = nfproc(ipni)
            ipi   = nfjpi (ipni)
            !
            IF( ipni ==   1  ) THEN   ;   iis0 =   1          ! domain  left side: as e-w comm already done -> from 1st column
            ELSE                      ;   iis0 =   1 + khls   ! default: -> from inner domain 
            ENDIF
            IF( ipni == jpni ) THEN   ;   iie0 = ipi          ! domain right side: as e-w comm already done -> until last column
            ELSE                      ;   iie0 = ipi - khls   ! default: -> until inner domain 
            ENDIF
            impp = nfimpp(ipni) - nfimpp(isendto(1))
            !
            IF(           iproc == -1 ) THEN   ! No neighbour (land proc that was suppressed)
               !
               SELECT CASE ( kfillmode )
               CASE ( jpfillnothing )               ! no filling 
               CASE ( jpfillcopy    )               ! filling with inner domain values
                  DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
                     DO jj = 1, ipj_s(jf)
                        ij1 = jj_b(jj,jf)
                        ij2 = jj_s(jj,jf)
                        DO ji = iis0, iie0
                           ztabr(impp+ji,ij1,jk,jl) = ptab(jf)%pt4d(Nis0,ij2,jk,jl)   ! chose to take the 1st iner domain point
                        END DO
                     END DO
                  END DO   ;   END DO   ;   END DO
               CASE ( jpfillcst     )               ! filling with constant value
                  DO jl = 1, ipl   ;   DO jk = 1, ipk
                     DO jj = 1, ipj_b
                        DO ji = iis0, iie0
                           ztabr(impp+ji,jj,jk,jl) = pfillval
                        END DO
                     END DO
                  END DO   ;   END DO
               END SELECT
               !
            ELSE IF( iproc == narea-1 ) THEN   ! get data from myself!
               !
               DO jf = 1, ipf   ;   DO jl = 1, ipl  ;   DO jk = 1, ipk
                  DO jj = 1, ipj_s(jf)
                     ij1 = jj_b(jj,jf)
                     ij2 = jj_s(jj,jf)
                     DO ji = iis0, iie0
                        ztabr(impp+ji,ij1,jk,jl) = ptab(jf)%pt4d(ji,ij2,jk,jl)
                     END DO
                  END DO
               END DO   ;   END DO   ;   END DO
               !
            ELSE                               ! get data from a neighbour trough communication
               !  
#if ! defined key_mpi_off
               CALL MPI_RECV( ztabw, ibuffsize, MPI_TYPE, iproc, 5, mpi_comm_oce, MPI_STATUS_IGNORE, ierr )
#endif
               DO jl = 1, ipl   ;   DO jk = 1, ipk
                  DO jj = 1, ipj_b
                     DO ji = iis0, iie0
                        ztabr(impp+ji,jj,jk,jl) = ztabw(ji,jj,jk,jl)
                     END DO
                  END DO
               END DO   ;   END DO
               
            ENDIF
            !
         END DO   ! nsndto
         !
         IF( ln_timing ) CALL tic_tac(.FALSE.)
         !
         ! North fold boundary condition
         !
         DO jf = 1, ipf
            ij1 = jj_b(       1 ,jf)
            ij2 = jj_b(ipj_s(jf),jf)
            CALL lbc_nfd_nogather( ptab(jf)%pt4d(:,:,:,:), ztabr(:,ij1:ij2,:,:), cd_nat(jf), psgn(jf), khls )
         END DO
         !
         DEALLOCATE( ztabr, ztabw, jj_s, jj_b, ipj_s )
         !
         DO jr = 1,nsndto
            iproc = nfproc(isendto(jr))
            IF( iproc /= narea-1 .AND. iproc /= -1 ) THEN
               CALL mpi_wait( ml_req_nf(jr), MPI_STATUS_IGNORE, ml_err )   ! put the wait at the very end just before the deallocate
            ENDIF
         END DO
         DEALLOCATE( ztabb )
         !
      ELSE                             !==  allgather exchanges  ==!
         !
         ! how many lines do we exchange at max? -> ipj    (no further optimizations in this case...)
         ipj =      khls + 2
         ! how many lines do we     need at max? -> ipj2   (no further optimizations in this case...)
         ipj2 = 2 * khls + 2
         !
         i0max = jpimax - 2 * khls
         ibuffsize = i0max * ipj * ipk * ipl * ipf
         ALLOCATE( znorthloc(i0max,ipj,ipk,ipl,ipf), znorthglo(i0max,ipj,ipk,ipl,ipf,ndim_rank_north) )
         !
         DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk               ! put in znorthloc ipj j-lines of ptab
            DO jj = 1, ipj
               ij2 = jpj - ipj2 + jj                        ! the first ipj lines of the last ipj2 lines
               DO ji = 1, Ni_0
                  ii2 = Nis0 - 1 + ji                       ! inner domain: Nis0 to Nie0
                  znorthloc(ji,jj,jk,jl,jf) = ptab(jf)%pt4d(ii2,ij2,jk,jl)
               END DO
               DO ji = Ni_0+1, i0max
                  znorthloc(ji,jj,jk,jl,jf) = HUGE(0._/**/PRECISION)   ! avoid sending uninitialized values (make sure we don't use it)
               END DO
            END DO
         END DO   ;   END DO   ;   END DO
         !
         ! start waiting time measurement
         IF( ln_timing ) CALL tic_tac(.TRUE.)
#if ! defined key_mpi_off
         CALL MPI_ALLGATHER( znorthloc, ibuffsize, MPI_TYPE, znorthglo, ibuffsize, MPI_TYPE, ncomm_north, ierr )
#endif
         ! stop waiting time measurement
         IF( ln_timing ) CALL tic_tac(.FALSE.)
         DEALLOCATE( znorthloc )
         ALLOCATE( ztabglo(ipf) )
         DO jf = 1, ipf
            ALLOCATE( ztabglo(jf)%pt4d(jpiglo,ipj2,ipk,ipl) )
         END DO
         !
         ! need to fill only the first ipj lines of ztabglo as lbc_nfd don't use the last khls lines
         ijnr = 0
         DO jr = 1, jpni                                                        ! recover the global north array
            iproc = nfproc(jr)
            impp  = nfimpp(jr)
            ipi   = nfjpi( jr) - 2 * khls                       ! corresponds to Ni_0 but for subdomain iproc
            IF( iproc == -1 ) THEN   ! No neighbour (land proc that was suppressed)
              !
               SELECT CASE ( kfillmode )
               CASE ( jpfillnothing )               ! no filling 
               CASE ( jpfillcopy    )               ! filling with inner domain values
                  DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
                     DO jj = 1, ipj
                        ij2 = jpj - ipj2 + jj                    ! the first ipj lines of the last ipj2 lines
                        DO ji = 1, ipi
                           ii1 = impp + khls + ji - 1            ! corresponds to mig(khls + ji) but for subdomain iproc
                           ztabglo(jf)%pt4d(ii1,jj,jk,jl) = ptab(jf)%pt4d(Nis0,ij2,jk,jl) ! chose to take the 1st inner domain point
                        END DO
                     END DO
                  END DO   ;   END DO   ;   END DO
               CASE ( jpfillcst     )               ! filling with constant value
                  DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
                     DO jj = 1, ipj
                        DO ji = 1, ipi
                           ii1 = impp + khls + ji - 1            ! corresponds to mig(khls + ji) but for subdomain iproc
                           ztabglo(jf)%pt4d(ii1,jj,jk,jl) = pfillval
                        END DO
                     END DO
                 END DO   ;   END DO   ;   END DO
               END SELECT
               !
            ELSE
               ijnr = ijnr + 1
               DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
                  DO jj = 1, ipj
                     DO ji = 1, ipi
                        ii1 = impp + khls + ji - 1               ! corresponds to mig(khls + ji) but for subdomain iproc
                        ztabglo(jf)%pt4d(ii1,jj,jk,jl) = znorthglo(ji,jj,jk,jl,jf,ijnr)
                     END DO
                  END DO
               END DO   ;   END DO   ;   END DO
            ENDIF
            !
         END DO   ! jpni
         DEALLOCATE( znorthglo )
         !
         DO jf = 1, ipf
            CALL lbc_nfd( ztabglo(jf:jf), cd_nat(jf:jf), psgn(jf:jf), khls, 1 )   ! North fold boundary condition
            DO jl = 1, ipl   ;   DO jk = 1, ipk                  ! e-w periodicity
               DO jj = 1, khls + 1
                  ij1 = ipj2 - (khls + 1) + jj                   ! need only the last khls + 1 lines until ipj2
                  ztabglo(jf)%pt4d(            1:  khls,ij1,jk,jl) = ztabglo(jf)%pt4d(jpiglo-2*khls+1:jpiglo-khls,ij1,jk,jl)
                  ztabglo(jf)%pt4d(jpiglo-khls+1:jpiglo,ij1,jk,jl) = ztabglo(jf)%pt4d(         khls+1:     2*khls,ij1,jk,jl)
               END DO
            END DO   ;   END DO
         END DO     
         !
         DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk               ! Scatter back to ARRAY_IN
            DO jj = 1, khls + 1
               ij1 = jpj  - (khls + 1) + jj   ! last khls + 1 lines until jpj
               ij2 = ipj2 - (khls + 1) + jj   ! last khls + 1 lines until ipj2
               DO ji= 1, jpi
                  ii2 = mig(ji)
                  ptab(jf)%pt4d(ji,ij1,jk,jl) = ztabglo(jf)%pt4d(ii2,ij2,jk,jl)
               END DO
            END DO
         END DO   ;   END DO   ;   END DO
         !
         DO jf = 1, ipf
            DEALLOCATE( ztabglo(jf)%pt4d )
         END DO
         DEALLOCATE( ztabglo )
         !
      ENDIF   ! l_north_nogather
      !
   END SUBROUTINE mpp_nfd_/**/PRECISION