MODULE bdyini
   !!======================================================================
   !!                       ***  MODULE  bdyini  ***
   !! Unstructured open boundaries : initialisation
   !!======================================================================
   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
   !!             -   !  2007-01  (D. Storkey) Update to use IOM module
   !!             -   !  2007-01  (D. Storkey) Tidal forcing
   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
   !!            3.3  !  2010-09  (E.O'Dea) updates for Shelf configurations
   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
   !!----------------------------------------------------------------------
#if defined key_bdy
   !!----------------------------------------------------------------------
   !!   'key_bdy'                     Unstructured Open Boundary Conditions
   !!----------------------------------------------------------------------
   !!   bdy_init       : Initialization of unstructured open boundaries
   !!----------------------------------------------------------------------
   USE timing          ! Timing
   USE oce             ! ocean dynamics and tracers variables
   USE dom_oce         ! ocean space and time domain
   USE bdy_oce         ! unstructured open boundary conditions
   USE in_out_manager  ! I/O units
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
   USE lib_mpp         ! for mpp_sum  
   USE iom             ! I/O

   IMPLICIT NONE
   PRIVATE

   PUBLIC   bdy_init   ! routine called in nemo_init

   !!----------------------------------------------------------------------
   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
   !! $Id: bdyini.F90 3424 2012-07-05 09:35:12Z davestorkey $ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS
   
   SUBROUTINE bdy_init
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE bdy_init  ***
      !!         
      !! ** Purpose :   Initialization of the dynamics and tracer fields with 
      !!              unstructured open boundaries.
      !!
      !! ** Method  :   Read initialization arrays (mask, indices) to identify 
      !!              an unstructured open boundary
      !!
      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
      !!----------------------------------------------------------------------      
      ! namelist variables
      !-------------------
      INTEGER, PARAMETER          :: jp_nseg = 100
      INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
      INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft
      INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft
      INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft
      INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft

      ! local variables
      !-------------------
      INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers
      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       -
      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       -
      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
      REAL   , POINTER  ::  flagu, flagv                   !    -   -
      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
      INTEGER, DIMENSION (2)                ::   kdimsz
      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays 
      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
      REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask
      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
      CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid
      !!
      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             &
         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &  
#if defined key_lim2
         &             nn_ice_lim2, nn_ice_lim2_dta,                       &
#endif
         &             ln_vol, nn_volctl, nn_rimwidth
      !!
      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             &
                             nbdysegw, jpiwob, jpjwdt, jpjwft,             &
                             nbdysegn, jpjnob, jpindt, jpinft,             &
                             nbdysegs, jpjsob, jpisdt, jpisft

      !!----------------------------------------------------------------------

      IF( nn_timing == 1 ) CALL timing_start('bdy_init')

      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )

      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
      IF(lwp) WRITE(numout,*) '~~~~~~~~'
      !

      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
         &                               ' and general open boundary condition are not compatible' )

      cgrid= (/'t','u','v'/)

      ! -----------------------------------------
      ! Initialise and read namelist parameters
      ! -----------------------------------------

      nb_bdy            = 0
      ln_coords_file(:) = .false.
      cn_coords_file(:) = ''
      ln_mask_file      = .false.
      cn_mask_file(:)   = ''
      nn_dyn2d(:)       = 0
      nn_dyn2d_dta(:)   = -1  ! uninitialised flag
      nn_dyn3d(:)       = 0
      nn_dyn3d_dta(:)   = -1  ! uninitialised flag
      nn_tra(:)         = 0
      nn_tra_dta(:)     = -1  ! uninitialised flag
#if defined key_lim2
      nn_ice_lim2(:)    = 0
      nn_ice_lim2_dta(:)= -1  ! uninitialised flag
#endif
      ln_vol            = .false.
      nn_volctl         = -1  ! uninitialised flag
      nn_rimwidth(:)    = -1  ! uninitialised flag

      REWIND( numnam )                    
      READ  ( numnam, nambdy )

      ! -----------------------------------------
      ! Check and write out namelist parameters
      ! -----------------------------------------

      !                                   ! control prints
      IF(lwp) WRITE(numout,*) '         nambdy'

      IF( nb_bdy .eq. 0 ) THEN 
        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
      ELSE
        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
      ENDIF

      DO ib_bdy = 1,nb_bdy
        IF(lwp) WRITE(numout,*) ' ' 
        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 

        IF( ln_coords_file(ib_bdy) ) THEN
           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
        ELSE
           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
        ENDIF
        IF(lwp) WRITE(numout,*)

        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
        SELECT CASE( nn_dyn2d(ib_bdy) )                  
          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
        END SELECT
        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   ! 
              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file'
              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
           END SELECT
        ENDIF
        IF(lwp) WRITE(numout,*)

        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
        SELECT CASE( nn_dyn3d(ib_bdy) )                  
          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
        END SELECT
        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN
           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   ! 
              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
           END SELECT
        ENDIF
        IF(lwp) WRITE(numout,*)

        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
        SELECT CASE( nn_tra(ib_bdy) )                  
          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
        END SELECT
        IF( nn_tra(ib_bdy) .gt. 0 ) THEN
           SELECT CASE( nn_tra_dta(ib_bdy) )                   ! 
              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
           END SELECT
        ENDIF
        IF(lwp) WRITE(numout,*)

#if defined key_lim2
        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
        SELECT CASE( nn_ice_lim2(ib_bdy) )                  
          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
        END SELECT
        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN 
           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   ! 
              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
           END SELECT
        ENDIF
        IF(lwp) WRITE(numout,*)
#endif

        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy)
        IF(lwp) WRITE(numout,*)

      ENDDO

     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
       IF(lwp) WRITE(numout,*)
       SELECT CASE ( nn_volctl )
         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
         CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
       END SELECT
       IF(lwp) WRITE(numout,*)
     ELSE
       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
       IF(lwp) WRITE(numout,*)
     ENDIF

      ! -------------------------------------------------
      ! Initialise indices arrays for open boundaries
      ! -------------------------------------------------

      ! Work out global dimensions of boundary data
      ! ---------------------------------------------
      REWIND( numnam )                    
      DO ib_bdy = 1, nb_bdy

         jpbdta = 1
         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
 
            ! No REWIND here because may need to read more than one nambdy_index namelist.
            READ  ( numnam, nambdy_index )

            ! Automatic boundary definition: if nbdysegX = -1
            ! set boundary to whole side of model domain.
            IF( nbdysege == -1 ) THEN
               nbdysege = 1
               jpieob(1) = jpiglo - 1
               jpjedt(1) = 2
               jpjeft(1) = jpjglo - 1
            ENDIF
            IF( nbdysegw == -1 ) THEN
               nbdysegw = 1
               jpiwob(1) = 2
               jpjwdt(1) = 2
               jpjwft(1) = jpjglo - 1
            ENDIF
            IF( nbdysegn == -1 ) THEN
               nbdysegn = 1
               jpjnob(1) = jpjglo - 1
               jpindt(1) = 2
               jpinft(1) = jpiglo - 1
            ENDIF
            IF( nbdysegs == -1 ) THEN
               nbdysegs = 1
               jpjsob(1) = 2
               jpisdt(1) = 2
               jpisft(1) = jpiglo - 1
            ENDIF

            nblendta(:,ib_bdy) = 0
            DO iseg = 1, nbdysege
               igrd = 1
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
               igrd = 2
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
               igrd = 3
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)               
            ENDDO
            DO iseg = 1, nbdysegw
               igrd = 1
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
               igrd = 2
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
               igrd = 3
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)               
            ENDDO
            DO iseg = 1, nbdysegn
               igrd = 1
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1               
               igrd = 2
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)               
               igrd = 3
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1
            ENDDO
            DO iseg = 1, nbdysegs
               igrd = 1
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
               igrd = 2
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg)
               igrd = 3
               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
            ENDDO

            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy)
            jpbdta = MAXVAL(nblendta(:,ib_bdy))               


         ELSE            ! Read size of arrays in boundary coordinates file.


            CALL iom_open( cn_coords_file(ib_bdy), inum )
            jpbdta = 1
            DO igrd = 1, jpbgrd
               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )  
               nblendta(igrd,ib_bdy) = kdimsz(1)
               jpbdta = MAX(jpbdta, kdimsz(1))
            ENDDO

         ENDIF 

      ENDDO ! ib_bdy

      ! Allocate arrays
      !---------------
      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
         &      nbrdta(jpbdta, jpbgrd, nb_bdy) )

      ALLOCATE( dta_global(jpbdta, 1, jpk) )

      ! Calculate global boundary index arrays or read in from file
      !------------------------------------------------------------
      REWIND( numnam )                    
      DO ib_bdy = 1, nb_bdy

         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters

            ! No REWIND here because may need to read more than one nambdy_index namelist.
            READ  ( numnam, nambdy_index )

            ! Automatic boundary definition: if nbdysegX = -1
            ! set boundary to whole side of model domain.
            IF( nbdysege == -1 ) THEN
               nbdysege = 1
               jpieob(1) = jpiglo - 1
               jpjedt(1) = 2
               jpjeft(1) = jpjglo - 1
            ENDIF
            IF( nbdysegw == -1 ) THEN
               nbdysegw = 1
               jpiwob(1) = 2
               jpjwdt(1) = 2
               jpjwft(1) = jpjglo - 1
            ENDIF
            IF( nbdysegn == -1 ) THEN
               nbdysegn = 1
               jpjnob(1) = jpjglo - 1
               jpindt(1) = 2
               jpinft(1) = jpiglo - 1
            ENDIF
            IF( nbdysegs == -1 ) THEN
               nbdysegs = 1
               jpjsob(1) = 2
               jpisdt(1) = 2
               jpisft(1) = jpiglo - 1
            ENDIF

            ! ------------ T points -------------
            igrd = 1  
            icount = 0
            DO ir = 1, nn_rimwidth(ib_bdy)
               ! east
               DO iseg = 1, nbdysege
                  DO ij = jpjedt(iseg), jpjeft(iseg)
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
                     nbjdta(icount, igrd, ib_bdy) = ij
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! west
               DO iseg = 1, nbdysegw
                  DO ij = jpjwdt(iseg), jpjwft(iseg)
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
                     nbjdta(icount, igrd, ib_bdy) = ij
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! north
               DO iseg = 1, nbdysegn
                  DO ii = jpindt(iseg), jpinft(iseg)
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = ii
                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! south
               DO iseg = 1, nbdysegs
                  DO ii = jpisdt(iseg), jpisft(iseg)
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = ii
                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
            ENDDO

            ! ------------ U points -------------
            igrd = 2  
            icount = 0
            DO ir = 1, nn_rimwidth(ib_bdy)
               ! east
               DO iseg = 1, nbdysege
                  DO ij = jpjedt(iseg), jpjeft(iseg)
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir
                     nbjdta(icount, igrd, ib_bdy) = ij
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! west
               DO iseg = 1, nbdysegw
                  DO ij = jpjwdt(iseg), jpjwft(iseg)
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
                     nbjdta(icount, igrd, ib_bdy) = ij
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! north
               DO iseg = 1, nbdysegn
                  DO ii = jpindt(iseg), jpinft(iseg) - 1
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = ii
                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! south
               DO iseg = 1, nbdysegs
                  DO ii = jpisdt(iseg), jpisft(iseg) - 1
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = ii
                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
            ENDDO

            ! ------------ V points -------------
            igrd = 3  
            icount = 0
            DO ir = 1, nn_rimwidth(ib_bdy)
               ! east
               DO iseg = 1, nbdysege
                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
                     nbjdta(icount, igrd, ib_bdy) = ij
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! west
               DO iseg = 1, nbdysegw
                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
                     nbjdta(icount, igrd, ib_bdy) = ij
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! north
               DO iseg = 1, nbdysegn
                  DO ii = jpindt(iseg), jpinft(iseg)
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = ii
                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
               ! south
               DO iseg = 1, nbdysegs
                  DO ii = jpisdt(iseg), jpisft(iseg)
                     icount = icount + 1
                     nbidta(icount, igrd, ib_bdy) = ii
                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
                     nbrdta(icount, igrd, ib_bdy) = ir
                  ENDDO
               ENDDO
            ENDDO

         ELSE            ! Read global index arrays from boundary coordinates file.

            DO igrd = 1, jpbgrd
               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
               DO ii = 1,nblendta(igrd,ib_bdy)
                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
               END DO
               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
               DO ii = 1,nblendta(igrd,ib_bdy)
                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
               END DO
               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
               DO ii = 1,nblendta(igrd,ib_bdy)
                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
               END DO

               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
               IF(lwp) WRITE(numout,*)
               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
               IF (ibr_max < nn_rimwidth(ib_bdy))   &
                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )

            END DO
            CALL iom_close( inum )

         ENDIF 

      ENDDO 

      ! Work out dimensions of boundary data on each processor
      ! ------------------------------------------------------
     
      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1

      DO ib_bdy = 1, nb_bdy
         DO igrd = 1, jpbgrd
            icount  = 0
            icountr = 0
            idx_bdy(ib_bdy)%nblen(igrd)    = 0
            idx_bdy(ib_bdy)%nblenrim(igrd) = 0
            DO ib = 1, nblendta(igrd,ib_bdy)
               ! check that data is in correct order in file
               ibm1 = MAX(1,ib-1)
               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
                    'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory')
                  ENDIF    
               ENDIF
               ! check if point is in local domain
               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
                  !
                  icount = icount  + 1
                  !
                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1
               ENDIF
            ENDDO
            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc        
         ENDDO  ! igrd

         ! Allocate index arrays for this boundary set
         !--------------------------------------------
         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )

         ! Dispatch mapping indices and discrete distances on each processor
         ! -----------------------------------------------------------------

         DO igrd = 1, jpbgrd
            icount  = 0
            ! Loop on rimwidth to ensure outermost points come first in the local arrays.
            DO ir=1, nn_rimwidth(ib_bdy)
               DO ib = 1, nblendta(igrd,ib_bdy)
                  ! check if point is in local domain and equals ir
                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                     !
                     icount = icount  + 1
                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
                  ENDIF
               ENDDO
            ENDDO
         ENDDO 

         ! Compute rim weights for FRS scheme
         ! ----------------------------------
         DO igrd = 1, jpbgrd
            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
!              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic
!              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear
            END DO
         END DO 

      ENDDO

      ! ------------------------------------------------------
      ! Initialise masks and find normal/tangential directions
      ! ------------------------------------------------------

      ! Read global 2D mask at T-points: bdytmask
      ! -----------------------------------------
      ! bdytmask = 1  on the computational domain AND on open boundaries
      !          = 0  elsewhere   
 
      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution
         zmask(         :                ,:) = 0.e0
         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0          
      ELSE IF( ln_mask_file ) THEN
         CALL iom_open( cn_mask_file, inum )
         CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
         CALL iom_close( inum )
      ELSE
         zmask(:,:) = 1.e0
      ENDIF

      DO ij = 1, nlcj      ! Save mask over local domain      
         DO ii = 1, nlci
            bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
         END DO
      END DO

      ! Derive mask on U and V grid from mask on T grid
      bdyumask(:,:) = 0.e0
      bdyvmask(:,:) = 0.e0
      DO ij=1, jpjm1
         DO ii=1, jpim1
            bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
            bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)  
         END DO
      END DO
      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.


      ! Mask corrections
      ! ----------------
      DO ik = 1, jpkm1
         DO ij = 1, jpj
            DO ii = 1, jpi
               tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
               umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
               vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
               bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
            END DO      
         END DO
      END DO

      DO ik = 1, jpkm1
         DO ij = 2, jpjm1
            DO ii = 2, jpim1
               fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
                  &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
            END DO      
         END DO
      END DO

      tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)             
      bdytmask(:,:) = tmask(:,:,1)

      ! bdy masks and bmask are now set to zero on boundary points:
      igrd = 1       ! In the free surface case, bmask is at T-points
      DO ib_bdy = 1, nb_bdy
        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
        ENDDO
      ENDDO
      !
      igrd = 1
      DO ib_bdy = 1, nb_bdy
        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)      
          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
        ENDDO
      ENDDO
      !
      igrd = 2
      DO ib_bdy = 1, nb_bdy
        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
        ENDDO
      ENDDO
      !
      igrd = 3
      DO ib_bdy = 1, nb_bdy
        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
        ENDDO
      ENDDO

      ! Lateral boundary conditions
      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )

      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components

         idx_bdy(ib_bdy)%flagu(:) = 0.e0
         idx_bdy(ib_bdy)%flagv(:) = 0.e0
         icount = 0 

         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
         !flagu =  0 : u is tangential
         !flagu =  1 : u is normal to the boundary and is direction is inward
  
         igrd = 2      ! u-component 
         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
            zefl = bdytmask(nbi  ,nbj)
            zwfl = bdytmask(nbi+1,nbj)
            IF( zefl + zwfl == 2 ) THEN
               icount = icount + 1
            ELSE
               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
            ENDIF
         END DO

         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
         !flagv =  0 : u is tangential
         !flagv =  1 : u is normal to the boundary and is direction is inward

         igrd = 3      ! v-component
         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
            znfl = bdytmask(nbi,nbj  )
            zsfl = bdytmask(nbi,nbj+1)
            IF( znfl + zsfl == 2 ) THEN
               icount = icount + 1
            ELSE
               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
            END IF
         END DO
 
         IF( icount /= 0 ) THEN
            IF(lwp) WRITE(numout,*)
            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
            IF(lwp) WRITE(numout,*) ' ========== '
            IF(lwp) WRITE(numout,*)
            nstop = nstop + 1
         ENDIF 
    
      ENDDO

      ! Compute total lateral surface for volume correction:
      ! ----------------------------------------------------
      bdysurftot = 0.e0 
      IF( ln_vol ) THEN  
         igrd = 2      ! Lateral surface at U-points
         DO ib_bdy = 1, nb_bdy
            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
               flagu => idx_bdy(ib_bdy)%flagu(ib)
               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
                  &                    * tmask_i(nbi  , nbj)                           &
                  &                    * tmask_i(nbi+1, nbj)                   
            ENDDO
         ENDDO

         igrd=3 ! Add lateral surface at V-points
         DO ib_bdy = 1, nb_bdy
            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
               flagv => idx_bdy(ib_bdy)%flagv(ib)
               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
                  &                    * tmask_i(nbi, nbj  )                           &
                  &                    * tmask_i(nbi, nbj+1)
            ENDDO
         ENDDO
         !
         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
      END IF   
      !
      ! Tidy up
      !--------
      DEALLOCATE(nbidta, nbjdta, nbrdta)

      IF( nn_timing == 1 ) CALL timing_stop('bdy_init')

   END SUBROUTINE bdy_init

#else
   !!---------------------------------------------------------------------------------
   !!   Dummy module                                   NO open boundaries
   !!---------------------------------------------------------------------------------
CONTAINS
   SUBROUTINE bdy_init      ! Dummy routine
   END SUBROUTINE bdy_init
#endif

   !!=================================================================================
END MODULE bdyini
