Ignore:
Timestamp:
2019-09-11T15:54:18+02:00 (15 months ago)
Author:
smasson
Message:

trunk: merge dev_r10984_HPC-13 into the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/BDY/bdyini.F90

    r10983 r11536  
    3333   PRIVATE 
    3434 
    35    PUBLIC   bdy_init   ! routine called in nemo_init 
     35   PUBLIC   bdy_init    ! routine called in nemo_init 
     36   PUBLIC   find_neib   ! routine called in bdy_nmn 
    3637 
    3738   INTEGER, PARAMETER ::   jp_nseg = 100   !  
    38    INTEGER, PARAMETER ::   nrimmax =  20   ! maximum rimwidth in structured 
    39                                                ! open boundary data files 
    4039   ! Straight open boundary segment parameters: 
    4140   INTEGER  ::   nbdysege, nbdysegw, nbdysegn, nbdysegs  
     
    6867         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
    6968         &             cn_ice, nn_ice_dta,                                     & 
    70          &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
    71          &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
     69         &             ln_vol, nn_volctl, nn_rimwidth 
    7270         ! 
    7371      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     
    7977      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries 
    8078      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901) 
    81 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp ) 
     79901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist' ) 
     80      ! make sur that all elements of the namelist variables have a default definition from namelist_ref 
     81      ln_coords_file (2:jp_bdy) = ln_coords_file (1) 
     82      cn_coords_file (2:jp_bdy) = cn_coords_file (1) 
     83      cn_dyn2d       (2:jp_bdy) = cn_dyn2d       (1) 
     84      nn_dyn2d_dta   (2:jp_bdy) = nn_dyn2d_dta   (1) 
     85      cn_dyn3d       (2:jp_bdy) = cn_dyn3d       (1) 
     86      nn_dyn3d_dta   (2:jp_bdy) = nn_dyn3d_dta   (1) 
     87      cn_tra         (2:jp_bdy) = cn_tra         (1) 
     88      nn_tra_dta     (2:jp_bdy) = nn_tra_dta     (1)     
     89      ln_tra_dmp     (2:jp_bdy) = ln_tra_dmp     (1) 
     90      ln_dyn3d_dmp   (2:jp_bdy) = ln_dyn3d_dmp   (1) 
     91      rn_time_dmp    (2:jp_bdy) = rn_time_dmp    (1) 
     92      rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1) 
     93      cn_ice         (2:jp_bdy) = cn_ice         (1) 
     94      nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1) 
    8295      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
    8396      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 ) 
    84 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp ) 
     97902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' ) 
    8598      IF(lwm) WRITE ( numond, nambdy ) 
    8699 
    87100      IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE.   ! forced for Agrif children 
     101 
     102      IF( nb_bdy == 0 ) ln_bdy = .FALSE. 
    88103       
    89104      ! ----------------------------------------- 
     
    96111         ! 
    97112         ! Open boundaries definition (arrays and masks) 
    98          CALL bdy_segs 
     113         CALL bdy_def 
     114         IF( ln_meshmask )   CALL bdy_meshwri() 
    99115         ! 
    100116         ! Open boundaries initialisation of external data arrays 
     
    114130 
    115131 
    116    SUBROUTINE bdy_segs 
     132   SUBROUTINE bdy_def 
    117133      !!---------------------------------------------------------------------- 
    118134      !!                 ***  ROUTINE bdy_init  *** 
     
    125141      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries 
    126142      !!----------------------------------------------------------------------       
    127       INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices 
    128       INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers 
     143      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, ir, iseg     ! dummy loop indices 
     144      INTEGER  ::   icount, icountr, icountr0, ibr_max     ! local integers 
     145      INTEGER  ::   ilen1                                  !   -       - 
    129146      INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       - 
    130       INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       - 
    131       INTEGER  ::   jpbdtau, jpbdtas                       !   -       - 
     147      INTEGER  ::   jpbdta                                 !   -       - 
    132148      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       - 
    133       INTEGER  ::   i_offset, j_offset                     !   -       - 
    134       INTEGER , POINTER  ::  nbi, nbj, nbr                 ! short cuts 
    135       REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields 
    136       REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    137       INTEGER, DIMENSION (2)                  ::   kdimsz 
    138       INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays  
    139       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
    140       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points 
    141       CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid 
    142       INTEGER :: com_east, com_west, com_south, com_north, jpk_max  ! Flags for boundaries sending 
    143       INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving 
    144       INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates 
    145       REAL(wp), TARGET, DIMENSION(jpi,jpj) ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat) 
    146       !! 
    147       CHARACTER(LEN=1)                     ::   ctypebdy   !     -        -  
    148       INTEGER                              ::   nbdyind, nbdybeg, nbdyend 
    149       !! 
    150       NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
    151       INTEGER  ::   ios                 ! Local integer output status for namelist read 
     149      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3           !   -       - 
     150      INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       - 
     151      INTEGER  ::   flagu, flagv                           ! short cuts 
     152      INTEGER  ::   nbdyind, nbdybeg, nbdyend 
     153      INTEGER              , DIMENSION(4)             ::   kdimsz 
     154      INTEGER              , DIMENSION(jpbgrd,jp_bdy) ::   nblendta          ! Length of index arrays  
     155      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbidta, nbjdta    ! Index arrays: i and j indices of bdy dta 
     156      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbrdta            ! Discrete distance from rim points 
     157      CHARACTER(LEN=1)     , DIMENSION(jpbgrd)        ::   cgrid 
     158      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zz_read                 ! work space for 2D global boundary data 
     159      REAL(wp), POINTER    , DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields 
     160      REAL(wp)             , DIMENSION(jpi,jpj) ::   zfmask   ! temporary fmask array excluding coastal boundary condition (shlat) 
     161      REAL(wp)             , DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array 
    152162      !!---------------------------------------------------------------------- 
    153163      ! 
     
    160170         &                               ' and general open boundary condition are not compatible' ) 
    161171 
    162       IF( nb_bdy == 0 ) THEN  
    163         IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.' 
    164       ELSE 
    165         IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 
     172      IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 
     173 
     174      DO ib_bdy = 1,nb_bdy 
     175 
     176         IF(lwp) THEN 
     177            WRITE(numout,*) ' '  
     178            WRITE(numout,*) '------ Open boundary data set ',ib_bdy,' ------'  
     179            IF( ln_coords_file(ib_bdy) ) THEN 
     180               WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 
     181            ELSE 
     182               WRITE(numout,*) 'Boundary defined in namelist.' 
     183            ENDIF 
     184            WRITE(numout,*) 
     185         ENDIF 
     186 
     187         ! barotropic bdy 
     188         !---------------- 
     189         IF(lwp) THEN 
     190            WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
     191            SELECT CASE( cn_dyn2d(ib_bdy) )                   
     192            CASE( 'none' )           ;   WRITE(numout,*) '      no open boundary condition'         
     193            CASE( 'frs' )            ;   WRITE(numout,*) '      Flow Relaxation Scheme' 
     194            CASE( 'flather' )        ;   WRITE(numout,*) '      Flather radiation condition' 
     195            CASE( 'orlanski' )       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     196            CASE( 'orlanski_npo' )   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     197            CASE DEFAULT             ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 
     198            END SELECT 
     199         ENDIF 
     200 
     201         dta_bdy(ib_bdy)%lneed_ssh   = cn_dyn2d(ib_bdy) == 'flather' 
     202         dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none' 
     203 
     204         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN 
     205            SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
     206            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'         
     207            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file' 
     208            CASE( 2 )      ;   WRITE(numout,*) '      tidal harmonic forcing taken from file' 
     209            CASE( 3 )      ;   WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files' 
     210            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
     211            END SELECT 
     212         ENDIF 
     213         IF ( dta_bdy(ib_bdy)%lneed_dyn2d .AND. nn_dyn2d_dta(ib_bdy) .GE. 2  .AND. .NOT.ln_tide ) THEN 
     214            CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' ) 
     215         ENDIF 
     216         IF(lwp) WRITE(numout,*) 
     217 
     218         ! baroclinic bdy 
     219         !---------------- 
     220         IF(lwp) THEN 
     221            WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
     222            SELECT CASE( cn_dyn3d(ib_bdy) )                   
     223            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'         
     224            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme' 
     225            CASE('specified')      ;   WRITE(numout,*) '      Specified value' 
     226            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions' 
     227            CASE('zerograd')       ;   WRITE(numout,*) '      Zero gradient for baroclinic velocities' 
     228            CASE('zero')           ;   WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     229            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     230            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     231            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' ) 
     232            END SELECT 
     233         ENDIF 
     234 
     235         dta_bdy(ib_bdy)%lneed_dyn3d = cn_dyn3d(ib_bdy) == 'frs'      .OR. cn_dyn3d(ib_bdy) == 'specified'   & 
     236            &                     .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' 
     237 
     238         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn3d ) THEN 
     239            SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !  
     240            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'         
     241            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file' 
     242            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 
     243            END SELECT 
     244         END IF 
     245 
     246         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
     247            IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 
     248               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
     249               ln_dyn3d_dmp(ib_bdy) = .false. 
     250            ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 
     251               CALL ctl_stop( 'Use FRS OR relaxation' ) 
     252            ELSE 
     253               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
     254               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     255               IF(rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     256               dta_bdy(ib_bdy)%lneed_dyn3d = .TRUE. 
     257            ENDIF 
     258         ELSE 
     259            IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities' 
     260         ENDIF 
     261         IF(lwp) WRITE(numout,*) 
     262 
     263         !    tra bdy 
     264         !---------------- 
     265         IF(lwp) THEN 
     266            WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
     267            SELECT CASE( cn_tra(ib_bdy) )                   
     268            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'         
     269            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme' 
     270            CASE('specified')      ;   WRITE(numout,*) '      Specified value' 
     271            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions' 
     272            CASE('runoff')         ;   WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
     273            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
     274            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
     275            CASE DEFAULT           ;   CALL ctl_stop( 'unrecognised value for cn_tra' ) 
     276            END SELECT 
     277         ENDIF 
     278 
     279         dta_bdy(ib_bdy)%lneed_tra = cn_tra(ib_bdy) == 'frs'       .OR. cn_tra(ib_bdy) == 'specified'   & 
     280            &                   .OR. cn_tra(ib_bdy) == 'orlanski'  .OR. cn_tra(ib_bdy) == 'orlanski_npo'  
     281 
     282         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_tra ) THEN 
     283            SELECT CASE( nn_tra_dta(ib_bdy) )                   !  
     284            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'         
     285            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file' 
     286            CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
     287            END SELECT 
     288         ENDIF 
     289 
     290         IF ( ln_tra_dmp(ib_bdy) ) THEN 
     291            IF ( cn_tra(ib_bdy) == 'none' ) THEN 
     292               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
     293               ln_tra_dmp(ib_bdy) = .false. 
     294            ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 
     295               CALL ctl_stop( 'Use FRS OR relaxation' ) 
     296            ELSE 
     297               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
     298               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     299               IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days' 
     300               IF(lwp.AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     301               dta_bdy(ib_bdy)%lneed_tra = .TRUE. 
     302            ENDIF 
     303         ELSE 
     304            IF(lwp) WRITE(numout,*) '      NO T/S relaxation' 
     305         ENDIF 
     306         IF(lwp) WRITE(numout,*) 
     307 
     308#if defined key_si3 
     309         IF(lwp) THEN 
     310            WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
     311            SELECT CASE( cn_ice(ib_bdy) )                   
     312            CASE('none')   ;   WRITE(numout,*) '      no open boundary condition'         
     313            CASE('frs')    ;   WRITE(numout,*) '      Flow Relaxation Scheme' 
     314            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' ) 
     315            END SELECT 
     316         ENDIF 
     317 
     318         dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none' 
     319 
     320         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN  
     321            SELECT CASE( nn_ice_dta(ib_bdy) )                   !  
     322            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'         
     323            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file' 
     324            CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' ) 
     325            END SELECT 
     326         ENDIF 
     327#else 
     328         dta_bdy(ib_bdy)%lneed_ice = .FALSE. 
     329#endif 
     330         ! 
     331         IF(lwp) WRITE(numout,*) 
     332         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
     333         IF(lwp) WRITE(numout,*) 
     334         ! 
     335      END DO   ! nb_bdy 
     336 
     337      IF( lwp ) THEN 
     338         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
     339            WRITE(numout,*) 'Volume correction applied at open boundaries' 
     340            WRITE(numout,*) 
     341            SELECT CASE ( nn_volctl ) 
     342            CASE( 1 )      ;   WRITE(numout,*) '      The total volume will be constant' 
     343            CASE( 0 )      ;   WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
     344            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
     345            END SELECT 
     346            WRITE(numout,*) 
     347            ! 
     348            ! sanity check if used with tides         
     349            IF( ln_tide ) THEN  
     350               WRITE(numout,*) ' The total volume correction is not working with tides. ' 
     351               WRITE(numout,*) ' Set ln_vol to .FALSE. ' 
     352               WRITE(numout,*) ' or ' 
     353               WRITE(numout,*) ' equilibriate your bdy input files ' 
     354               CALL ctl_stop( 'The total volume correction is not working with tides.' ) 
     355            END IF 
     356         ELSE 
     357            WRITE(numout,*) 'No volume correction applied at open boundaries' 
     358            WRITE(numout,*) 
     359         ENDIF 
    166360      ENDIF 
    167  
    168       DO ib_bdy = 1,nb_bdy 
    169         IF(lwp) WRITE(numout,*) ' '  
    170         IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'  
    171  
    172         IF( ln_coords_file(ib_bdy) ) THEN 
    173            IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy)) 
    174         ELSE 
    175            IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.' 
    176         ENDIF 
    177         IF(lwp) WRITE(numout,*) 
    178  
    179         IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
    180         SELECT CASE( cn_dyn2d(ib_bdy) )                   
    181           CASE( 'none' )          
    182              IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    183              dta_bdy(ib_bdy)%ll_ssh = .false. 
    184              dta_bdy(ib_bdy)%ll_u2d = .false. 
    185              dta_bdy(ib_bdy)%ll_v2d = .false. 
    186           CASE( 'frs' )           
    187              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    188              dta_bdy(ib_bdy)%ll_ssh = .false. 
    189              dta_bdy(ib_bdy)%ll_u2d = .true. 
    190              dta_bdy(ib_bdy)%ll_v2d = .true. 
    191           CASE( 'flather' )       
    192              IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    193              dta_bdy(ib_bdy)%ll_ssh = .true. 
    194              dta_bdy(ib_bdy)%ll_u2d = .true. 
    195              dta_bdy(ib_bdy)%ll_v2d = .true. 
    196           CASE( 'orlanski' )      
    197              IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
    198              dta_bdy(ib_bdy)%ll_ssh = .false. 
    199              dta_bdy(ib_bdy)%ll_u2d = .true. 
    200              dta_bdy(ib_bdy)%ll_v2d = .true. 
    201           CASE( 'orlanski_npo' )  
    202              IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
    203              dta_bdy(ib_bdy)%ll_ssh = .false. 
    204              dta_bdy(ib_bdy)%ll_u2d = .true. 
    205              dta_bdy(ib_bdy)%ll_v2d = .true. 
    206           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 
    207         END SELECT 
    208         IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    209            SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
    210               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    211               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    212               CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file' 
    213               CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files' 
    214               CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
    215            END SELECT 
    216            IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.ln_tide)) THEN 
    217              CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' ) 
    218            ENDIF 
    219         ENDIF 
    220         IF(lwp) WRITE(numout,*) 
    221  
    222         IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
    223         SELECT CASE( cn_dyn3d(ib_bdy) )                   
    224           CASE('none') 
    225              IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    226              dta_bdy(ib_bdy)%ll_u3d = .false. 
    227              dta_bdy(ib_bdy)%ll_v3d = .false. 
    228           CASE('frs')        
    229              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    230              dta_bdy(ib_bdy)%ll_u3d = .true. 
    231              dta_bdy(ib_bdy)%ll_v3d = .true. 
    232           CASE('specified') 
    233              IF(lwp) WRITE(numout,*) '      Specified value' 
    234              dta_bdy(ib_bdy)%ll_u3d = .true. 
    235              dta_bdy(ib_bdy)%ll_v3d = .true. 
    236           CASE('neumann') 
    237              IF(lwp) WRITE(numout,*) '      Neumann conditions' 
    238              dta_bdy(ib_bdy)%ll_u3d = .false. 
    239              dta_bdy(ib_bdy)%ll_v3d = .false. 
    240           CASE('zerograd') 
    241              IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities' 
    242              dta_bdy(ib_bdy)%ll_u3d = .false. 
    243              dta_bdy(ib_bdy)%ll_v3d = .false. 
    244           CASE('zero') 
    245              IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
    246              dta_bdy(ib_bdy)%ll_u3d = .false. 
    247              dta_bdy(ib_bdy)%ll_v3d = .false. 
    248           CASE('orlanski') 
    249              IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
    250              dta_bdy(ib_bdy)%ll_u3d = .true. 
    251              dta_bdy(ib_bdy)%ll_v3d = .true. 
    252           CASE('orlanski_npo') 
    253              IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
    254              dta_bdy(ib_bdy)%ll_u3d = .true. 
    255              dta_bdy(ib_bdy)%ll_v3d = .true. 
    256           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' ) 
    257         END SELECT 
    258         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    259            SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !  
    260               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    261               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    262               CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' ) 
    263            END SELECT 
    264         ENDIF 
    265  
    266         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
    267            IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 
    268               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
    269               ln_dyn3d_dmp(ib_bdy)=.false. 
    270            ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 
    271               CALL ctl_stop( 'Use FRS OR relaxation' ) 
    272            ELSE 
    273               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
    274               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    275               IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
    276               dta_bdy(ib_bdy)%ll_u3d = .true. 
    277               dta_bdy(ib_bdy)%ll_v3d = .true. 
    278            ENDIF 
    279         ELSE 
    280            IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities' 
    281         ENDIF 
    282         IF(lwp) WRITE(numout,*) 
    283  
    284         IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
    285         SELECT CASE( cn_tra(ib_bdy) )                   
    286           CASE('none') 
    287              IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    288              dta_bdy(ib_bdy)%ll_tem = .false. 
    289              dta_bdy(ib_bdy)%ll_sal = .false. 
    290           CASE('frs') 
    291              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    292              dta_bdy(ib_bdy)%ll_tem = .true. 
    293              dta_bdy(ib_bdy)%ll_sal = .true. 
    294           CASE('specified') 
    295              IF(lwp) WRITE(numout,*) '      Specified value' 
    296              dta_bdy(ib_bdy)%ll_tem = .true. 
    297              dta_bdy(ib_bdy)%ll_sal = .true. 
    298           CASE('neumann') 
    299              IF(lwp) WRITE(numout,*) '      Neumann conditions' 
    300              dta_bdy(ib_bdy)%ll_tem = .false. 
    301              dta_bdy(ib_bdy)%ll_sal = .false. 
    302           CASE('runoff') 
    303              IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
    304              dta_bdy(ib_bdy)%ll_tem = .false. 
    305              dta_bdy(ib_bdy)%ll_sal = .false. 
    306           CASE('orlanski') 
    307              IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging' 
    308              dta_bdy(ib_bdy)%ll_tem = .true. 
    309              dta_bdy(ib_bdy)%ll_sal = .true. 
    310           CASE('orlanski_npo') 
    311              IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging' 
    312              dta_bdy(ib_bdy)%ll_tem = .true. 
    313              dta_bdy(ib_bdy)%ll_sal = .true. 
    314           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_tra' ) 
    315         END SELECT 
    316         IF( cn_tra(ib_bdy) /= 'none' ) THEN 
    317            SELECT CASE( nn_tra_dta(ib_bdy) )                   !  
    318               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    319               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    320               CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
    321            END SELECT 
    322         ENDIF 
    323  
    324         IF ( ln_tra_dmp(ib_bdy) ) THEN 
    325            IF ( cn_tra(ib_bdy) == 'none' ) THEN 
    326               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
    327               ln_tra_dmp(ib_bdy)=.false. 
    328            ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 
    329               CALL ctl_stop( 'Use FRS OR relaxation' ) 
    330            ELSE 
    331               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
    332               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    333               IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days' 
    334               IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
    335               dta_bdy(ib_bdy)%ll_tem = .true. 
    336               dta_bdy(ib_bdy)%ll_sal = .true. 
    337            ENDIF 
    338         ELSE 
    339            IF(lwp) WRITE(numout,*) '      NO T/S relaxation' 
    340         ENDIF 
    341         IF(lwp) WRITE(numout,*) 
    342  
    343 #if defined key_si3 
    344          IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
    345          SELECT CASE( cn_ice(ib_bdy) )                   
    346          CASE('none') 
    347              IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    348              dta_bdy(ib_bdy)%ll_a_i = .false. 
    349              dta_bdy(ib_bdy)%ll_h_i = .false. 
    350              dta_bdy(ib_bdy)%ll_h_s = .false. 
    351          CASE('frs') 
    352              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    353              dta_bdy(ib_bdy)%ll_a_i = .true. 
    354              dta_bdy(ib_bdy)%ll_h_i = .true. 
    355              dta_bdy(ib_bdy)%ll_h_s = .true. 
    356          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' ) 
    357          END SELECT 
    358         IF( cn_ice(ib_bdy) /= 'none' ) THEN  
    359            SELECT CASE( nn_ice_dta(ib_bdy) )                   !  
    360               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    361               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    362               CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' ) 
    363            END SELECT 
    364         ENDIF 
    365         IF(lwp) WRITE(numout,*) 
    366         IF(lwp) WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)          
    367         IF(lwp) WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)          
    368         IF(lwp) WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)          
    369 #endif 
    370  
    371         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
    372         IF(lwp) WRITE(numout,*) 
    373          ! 
    374       END DO 
    375  
    376      IF( nb_bdy > 0 ) THEN 
    377         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
    378           IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
    379           IF(lwp) WRITE(numout,*) 
    380           SELECT CASE ( nn_volctl ) 
    381             CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    382             CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    383             CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
    384           END SELECT 
    385           IF(lwp) WRITE(numout,*) 
    386           ! 
    387           ! sanity check if used with tides         
    388           IF( ln_tide ) THEN  
    389              IF(lwp) WRITE(numout,*) ' The total volume correction is not working with tides. ' 
    390              IF(lwp) WRITE(numout,*) ' Set ln_vol to .FALSE. ' 
    391              IF(lwp) WRITE(numout,*) ' or ' 
    392              IF(lwp) WRITE(numout,*) ' equilibriate your bdy input files ' 
    393              CALL ctl_stop( 'The total volume correction is not working with tides.' ) 
    394           END IF 
    395         ELSE 
    396           IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
    397           IF(lwp) WRITE(numout,*) 
    398         ENDIF 
    399         IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN 
    400            IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***' 
    401         ELSE 
    402            IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***' 
    403         ENDIF 
    404      ENDIF 
    405361 
    406362      ! ------------------------------------------------- 
     
    408364      ! ------------------------------------------------- 
    409365 
    410       ! Work out global dimensions of boundary data 
    411       ! --------------------------------------------- 
    412366      REWIND( numnam_cfg )      
    413  
    414367      nblendta(:,:) = 0 
    415368      nbdysege = 0 
     
    417370      nbdysegn = 0 
    418371      nbdysegs = 0 
    419       icount   = 0 ! count user defined segments 
    420       ! Dimensions below are used to allocate arrays to read external data 
    421       jpbdtas = 1 ! Maximum size of boundary data (structured case) 
    422       jpbdtau = 1 ! Maximum size of boundary data (unstructured case) 
    423  
     372 
     373      ! Define all boundaries  
     374      ! --------------------- 
    424375      DO ib_bdy = 1, nb_bdy 
    425  
    426          IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
    427   
    428             icount = icount + 1 
    429             ! No REWIND here because may need to read more than one nambdy_index namelist. 
    430             ! Read only namelist_cfg to avoid unseccessfull overwrite  
    431             ! keep full control of the configuration namelist 
    432             READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 ) 
    433 904         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist', lwp ) 
    434             IF(lwm) WRITE ( numond, nambdy_index ) 
    435  
    436             SELECT CASE ( TRIM(ctypebdy) ) 
    437               CASE( 'N' ) 
    438                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    439                     nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain. 
    440                     nbdybeg  = 2 
    441                     nbdyend  = jpiglo - 1 
    442                  ENDIF 
    443                  nbdysegn = nbdysegn + 1 
    444                  npckgn(nbdysegn) = ib_bdy ! Save bdy package number 
    445                  jpjnob(nbdysegn) = nbdyind 
    446                  jpindt(nbdysegn) = nbdybeg 
    447                  jpinft(nbdysegn) = nbdyend 
    448                  ! 
    449               CASE( 'S' ) 
    450                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    451                     nbdyind  = 2           ! set boundary to whole side of model domain. 
    452                     nbdybeg  = 2 
    453                     nbdyend  = jpiglo - 1 
    454                  ENDIF 
    455                  nbdysegs = nbdysegs + 1 
    456                  npckgs(nbdysegs) = ib_bdy ! Save bdy package number 
    457                  jpjsob(nbdysegs) = nbdyind 
    458                  jpisdt(nbdysegs) = nbdybeg 
    459                  jpisft(nbdysegs) = nbdyend 
    460                  ! 
    461               CASE( 'E' ) 
    462                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    463                     nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain. 
    464                     nbdybeg  = 2 
    465                     nbdyend  = jpjglo - 1 
    466                  ENDIF 
    467                  nbdysege = nbdysege + 1  
    468                  npckge(nbdysege) = ib_bdy ! Save bdy package number 
    469                  jpieob(nbdysege) = nbdyind 
    470                  jpjedt(nbdysege) = nbdybeg 
    471                  jpjeft(nbdysege) = nbdyend 
    472                  ! 
    473               CASE( 'W' ) 
    474                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    475                     nbdyind  = 2           ! set boundary to whole side of model domain. 
    476                     nbdybeg  = 2 
    477                     nbdyend  = jpjglo - 1 
    478                  ENDIF 
    479                  nbdysegw = nbdysegw + 1 
    480                  npckgw(nbdysegw) = ib_bdy ! Save bdy package number 
    481                  jpiwob(nbdysegw) = nbdyind 
    482                  jpjwdt(nbdysegw) = nbdybeg 
    483                  jpjwft(nbdysegw) = nbdyend 
    484                  ! 
    485               CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 
    486             END SELECT 
    487  
    488             ! For simplicity we assume that in case of straight bdy, arrays have the same length 
    489             ! (even if it is true that last tangential velocity points 
    490             ! are useless). This simplifies a little bit boundary data format (and agrees with format 
    491             ! used so far in obc package) 
    492  
    493             nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy) 
    494             jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1)) 
    495             IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) & 
    496             & CALL ctl_stop( 'rimwidth must be lower than nrimmax' ) 
    497  
    498          ELSE            ! Read size of arrays in boundary coordinates file. 
     376         ! 
     377         IF( .NOT. ln_coords_file(ib_bdy) ) THEN     ! build bdy coordinates with segments defined in namelist 
     378 
     379            CALL bdy_read_seg( ib_bdy, nblendta(:,ib_bdy) ) 
     380 
     381         ELSE                                        ! Read size of arrays in boundary coordinates file. 
     382             
    499383            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    500384            DO igrd = 1, jpbgrd 
    501385               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    502386               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz) 
    503                jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz)) 
    504387            END DO 
    505388            CALL iom_close( inum ) 
    506             ! 
    507          ENDIF  
     389         ENDIF 
    508390         ! 
    509391      END DO ! ib_bdy 
    510392 
    511       IF (nb_bdy>0) THEN 
    512          jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
    513  
    514          ! Allocate arrays 
    515          !--------------- 
    516          ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
    517             &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
    518           
    519          jpk_max = MAXVAL(nb_jpk_bdy) 
    520          jpk_max = MAX(jpk_max, jpk) 
    521  
    522          ALLOCATE( dta_global(jpbdtau, 1, jpk_max) ) 
    523          ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO 
    524          ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO 
    525  
    526          IF ( icount>0 ) THEN 
    527             ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) ) 
    528             ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO 
    529             ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO   
    530          ENDIF 
    531          !  
    532       ENDIF 
    533  
    534393      ! Now look for crossings in user (namelist) defined open boundary segments: 
    535       !-------------------------------------------------------------------------- 
    536       IF( icount>0 )   CALL bdy_ctl_seg 
    537  
     394      IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0)   CALL bdy_ctl_seg 
     395       
     396      ! Allocate arrays 
     397      !--------------- 
     398      jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
     399      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
     400     
    538401      ! Calculate global boundary index arrays or read in from file 
    539402      !------------------------------------------------------------                
     
    543406         IF( ln_coords_file(ib_bdy) ) THEN 
    544407            ! 
     408            ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) )           
    545409            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
     410            ! 
    546411            DO igrd = 1, jpbgrd 
    547                CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     412               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 
    548413               DO ii = 1,nblendta(igrd,ib_bdy) 
    549                   nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     414                  nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 
    550415               END DO 
    551                CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     416               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 
    552417               DO ii = 1,nblendta(igrd,ib_bdy) 
    553                   nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     418                  nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 
    554419               END DO 
    555                CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     420               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) ) 
    556421               DO ii = 1,nblendta(igrd,ib_bdy) 
    557                   nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) ) 
     422                  nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) ) 
    558423               END DO 
    559424               ! 
     
    563428               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy) 
    564429               IF (ibr_max < nn_rimwidth(ib_bdy))   & 
    565                      CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
    566             END DO 
     430                  CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
     431            END DO 
     432            ! 
    567433            CALL iom_close( inum ) 
     434            DEALLOCATE( zz_read ) 
    568435            ! 
    569          ENDIF  
    570          ! 
    571       END DO       
    572      
     436         ENDIF 
     437         ! 
     438      END DO 
     439 
    573440      ! 2. Now fill indices corresponding to straight open boundary arrays: 
    574       ! East 
    575       !----- 
    576       DO iseg = 1, nbdysege 
    577          ib_bdy = npckge(iseg) 
    578          ! 
    579          ! ------------ T points ------------- 
    580          igrd=1 
    581          icount=0 
    582          DO ir = 1, nn_rimwidth(ib_bdy) 
    583             DO ij = jpjedt(iseg), jpjeft(iseg) 
    584                icount = icount + 1 
    585                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
    586                nbjdta(icount, igrd, ib_bdy) = ij 
    587                nbrdta(icount, igrd, ib_bdy) = ir 
    588             ENDDO 
    589          ENDDO 
    590          ! 
    591          ! ------------ U points ------------- 
    592          igrd=2 
    593          icount=0 
    594          DO ir = 1, nn_rimwidth(ib_bdy) 
    595             DO ij = jpjedt(iseg), jpjeft(iseg) 
    596                icount = icount + 1 
    597                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 
    598                nbjdta(icount, igrd, ib_bdy) = ij 
    599                nbrdta(icount, igrd, ib_bdy) = ir 
    600             ENDDO 
    601          ENDDO 
    602          ! 
    603          ! ------------ V points ------------- 
    604          igrd=3 
    605          icount=0 
    606          DO ir = 1, nn_rimwidth(ib_bdy) 
    607 !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
    608             DO ij = jpjedt(iseg), jpjeft(iseg) 
    609                icount = icount + 1 
    610                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
    611                nbjdta(icount, igrd, ib_bdy) = ij 
    612                nbrdta(icount, igrd, ib_bdy) = ir 
    613             ENDDO 
    614             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    615             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    616          ENDDO 
    617       ENDDO 
    618       ! 
    619       ! West 
    620       !----- 
    621       DO iseg = 1, nbdysegw 
    622          ib_bdy = npckgw(iseg) 
    623          ! 
    624          ! ------------ T points ------------- 
    625          igrd=1 
    626          icount=0 
    627          DO ir = 1, nn_rimwidth(ib_bdy) 
    628             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    629                icount = icount + 1 
    630                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    631                nbjdta(icount, igrd, ib_bdy) = ij 
    632                nbrdta(icount, igrd, ib_bdy) = ir 
    633             ENDDO 
    634          ENDDO 
    635          ! 
    636          ! ------------ U points ------------- 
    637          igrd=2 
    638          icount=0 
    639          DO ir = 1, nn_rimwidth(ib_bdy) 
    640             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    641                icount = icount + 1 
    642                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    643                nbjdta(icount, igrd, ib_bdy) = ij 
    644                nbrdta(icount, igrd, ib_bdy) = ir 
    645             ENDDO 
    646          ENDDO 
    647          ! 
    648          ! ------------ V points ------------- 
    649          igrd=3 
    650          icount=0 
    651          DO ir = 1, nn_rimwidth(ib_bdy) 
    652 !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
    653             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    654                icount = icount + 1 
    655                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    656                nbjdta(icount, igrd, ib_bdy) = ij 
    657                nbrdta(icount, igrd, ib_bdy) = ir 
    658             ENDDO 
    659             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    660             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    661          ENDDO 
    662       ENDDO 
    663       ! 
    664       ! North 
    665       !----- 
    666       DO iseg = 1, nbdysegn 
    667          ib_bdy = npckgn(iseg) 
    668          ! 
    669          ! ------------ T points ------------- 
    670          igrd=1 
    671          icount=0 
    672          DO ir = 1, nn_rimwidth(ib_bdy) 
    673             DO ii = jpindt(iseg), jpinft(iseg) 
    674                icount = icount + 1 
    675                nbidta(icount, igrd, ib_bdy) = ii 
    676                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir  
    677                nbrdta(icount, igrd, ib_bdy) = ir 
    678             ENDDO 
    679          ENDDO 
    680          ! 
    681          ! ------------ U points ------------- 
    682          igrd=2 
    683          icount=0 
    684          DO ir = 1, nn_rimwidth(ib_bdy) 
    685 !            DO ii = jpindt(iseg), jpinft(iseg) - 1 
    686             DO ii = jpindt(iseg), jpinft(iseg) 
    687                icount = icount + 1 
    688                nbidta(icount, igrd, ib_bdy) = ii 
    689                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
    690                nbrdta(icount, igrd, ib_bdy) = ir 
    691             ENDDO 
    692             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    693             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    694          ENDDO 
    695          ! 
    696          ! ------------ V points ------------- 
    697          igrd=3 
    698          icount=0 
    699          DO ir = 1, nn_rimwidth(ib_bdy) 
    700             DO ii = jpindt(iseg), jpinft(iseg) 
    701                icount = icount + 1 
    702                nbidta(icount, igrd, ib_bdy) = ii 
    703                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 
    704                nbrdta(icount, igrd, ib_bdy) = ir 
    705             ENDDO 
    706          ENDDO 
    707       ENDDO 
    708       ! 
    709       ! South 
    710       !----- 
    711       DO iseg = 1, nbdysegs 
    712          ib_bdy = npckgs(iseg) 
    713          ! 
    714          ! ------------ T points ------------- 
    715          igrd=1 
    716          icount=0 
    717          DO ir = 1, nn_rimwidth(ib_bdy) 
    718             DO ii = jpisdt(iseg), jpisft(iseg) 
    719                icount = icount + 1 
    720                nbidta(icount, igrd, ib_bdy) = ii 
    721                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    722                nbrdta(icount, igrd, ib_bdy) = ir 
    723             ENDDO 
    724          ENDDO 
    725          ! 
    726          ! ------------ U points ------------- 
    727          igrd=2 
    728          icount=0 
    729          DO ir = 1, nn_rimwidth(ib_bdy) 
    730 !            DO ii = jpisdt(iseg), jpisft(iseg) - 1 
    731             DO ii = jpisdt(iseg), jpisft(iseg) 
    732                icount = icount + 1 
    733                nbidta(icount, igrd, ib_bdy) = ii 
    734                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    735                nbrdta(icount, igrd, ib_bdy) = ir 
    736             ENDDO 
    737             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    738             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    739          ENDDO 
    740          ! 
    741          ! ------------ V points ------------- 
    742          igrd=3 
    743          icount=0 
    744          DO ir = 1, nn_rimwidth(ib_bdy) 
    745             DO ii = jpisdt(iseg), jpisft(iseg) 
    746                icount = icount + 1 
    747                nbidta(icount, igrd, ib_bdy) = ii 
    748                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    749                nbrdta(icount, igrd, ib_bdy) = ir 
    750             ENDDO 
    751          ENDDO 
    752       ENDDO 
     441      CALL bdy_coords_seg( nbidta, nbjdta, nbrdta ) 
    753442 
    754443      !  Deal with duplicated points 
     
    764453                     DO ib2 = 1, nblendta(igrd,ib_bdy2) 
    765454                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. & 
    766                         &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 
    767 !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &  
    768 !                                                       &              nbidta(ib1, igrd, ib_bdy1),      &  
    769 !                                                       &              nbjdta(ib2, igrd, ib_bdy2) 
     455                           &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 
     456                           !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &  
     457                           !                                                       &              nbidta(ib1, igrd, ib_bdy1),      &  
     458                           !                                                       &              nbjdta(ib2, igrd, ib_bdy2) 
    770459                           ! keep only points with the lowest distance to boundary: 
    771460                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN 
    772                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 
    773                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 
     461                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 
     462                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 
    774463                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN 
    775                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    776                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    777                            ! Arbitrary choice if distances are the same: 
     464                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     465                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     466                              ! Arbitrary choice if distances are the same: 
    778467                           ELSE 
    779                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    780                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     468                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     469                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    781470                           ENDIF 
    782471                        END IF 
     
    787476         END DO 
    788477      END DO 
    789  
    790       ! Work out dimensions of boundary data on each processor 
    791       ! ------------------------------------------------------ 
    792  
    793       ! Rather assume that boundary data indices are given on global domain 
    794       ! TO BE DISCUSSED ? 
    795 !      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
    796 !      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
    797 !      is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
    798 !      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1       
    799       iwe = mig(1) - 1 + 2         ! if monotasking and no zoom, iw=2 
    800       ies = mig(1) + nlci-1 - 1  ! if monotasking and no zoom, ie=jpim1 
    801       iso = mjg(1) - 1 + 2         ! if monotasking and no zoom, is=2 
    802       ino = mjg(1) + nlcj-1 - 1  ! if monotasking and no zoom, in=jpjm1 
    803  
    804       ALLOCATE( nbondi_bdy(nb_bdy)) 
    805       ALLOCATE( nbondj_bdy(nb_bdy)) 
    806       nbondi_bdy(:)=2 
    807       nbondj_bdy(:)=2 
    808       ALLOCATE( nbondi_bdy_b(nb_bdy)) 
    809       ALLOCATE( nbondj_bdy_b(nb_bdy)) 
    810       nbondi_bdy_b(:)=2 
    811       nbondj_bdy_b(:)=2 
    812  
    813       ! Work out dimensions of boundary data on each neighbour process 
    814       IF(nbondi == 0) THEN 
    815          iw_b(1) = 1 + nimppt(nowe+1) 
    816          ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3 
    817          is_b(1) = 1 + njmppt(nowe+1) 
    818          in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3 
    819  
    820          iw_b(2) = 1 + nimppt(noea+1) 
    821          ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3 
    822          is_b(2) = 1 + njmppt(noea+1) 
    823          in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3 
    824       ELSEIF(nbondi == 1) THEN 
    825          iw_b(1) = 1 + nimppt(nowe+1) 
    826          ie_b(1) = 1 + nimppt(nowe+1)+nlcit(nowe+1)-3 
    827          is_b(1) = 1 + njmppt(nowe+1) 
    828          in_b(1) = 1 + njmppt(nowe+1)+nlcjt(nowe+1)-3 
    829       ELSEIF(nbondi == -1) THEN 
    830          iw_b(2) = 1 + nimppt(noea+1) 
    831          ie_b(2) = 1 + nimppt(noea+1)+nlcit(noea+1)-3 
    832          is_b(2) = 1 + njmppt(noea+1) 
    833          in_b(2) = 1 + njmppt(noea+1)+nlcjt(noea+1)-3 
    834       ENDIF 
    835  
    836       IF(nbondj == 0) THEN 
    837          iw_b(3) = 1 + nimppt(noso+1) 
    838          ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3 
    839          is_b(3) = 1 + njmppt(noso+1) 
    840          in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3 
    841  
    842          iw_b(4) = 1 + nimppt(nono+1) 
    843          ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3 
    844          is_b(4) = 1 + njmppt(nono+1) 
    845          in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3 
    846       ELSEIF(nbondj == 1) THEN 
    847          iw_b(3) = 1 + nimppt(noso+1) 
    848          ie_b(3) = 1 + nimppt(noso+1)+nlcit(noso+1)-3 
    849          is_b(3) = 1 + njmppt(noso+1) 
    850          in_b(3) = 1 + njmppt(noso+1)+nlcjt(noso+1)-3 
    851       ELSEIF(nbondj == -1) THEN 
    852          iw_b(4) = 1 + nimppt(nono+1) 
    853          ie_b(4) = 1 + nimppt(nono+1)+nlcit(nono+1)-3 
    854          is_b(4) = 1 + njmppt(nono+1) 
    855          in_b(4) = 1 + njmppt(nono+1)+nlcjt(nono+1)-3 
    856       ENDIF 
    857  
     478      ! 
     479      ! Find lenght of boundaries and rim on local mpi domain 
     480      !------------------------------------------------------ 
     481      ! 
     482      iwe = mig(1) 
     483      ies = mig(jpi) 
     484      iso = mjg(1)  
     485      ino = mjg(jpj)  
     486      ! 
    858487      DO ib_bdy = 1, nb_bdy 
    859488         DO igrd = 1, jpbgrd 
    860             icount  = 0 
    861             icountr = 0 
    862             idx_bdy(ib_bdy)%nblen(igrd)    = 0 
    863             idx_bdy(ib_bdy)%nblenrim(igrd) = 0 
     489            icount   = 0   ! initialization of local bdy length 
     490            icountr  = 0   ! initialization of local rim 0 and rim 1 bdy length 
     491            icountr0 = 0   ! initialization of local rim 0 bdy length 
     492            idx_bdy(ib_bdy)%nblen(igrd)     = 0 
     493            idx_bdy(ib_bdy)%nblenrim(igrd)  = 0 
     494            idx_bdy(ib_bdy)%nblenrim0(igrd) = 0 
    864495            DO ib = 1, nblendta(igrd,ib_bdy) 
    865496               ! check that data is in correct order in file 
    866                ibm1 = MAX(1,ib-1) 
    867                IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc... 
    868                   IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 
     497               IF( ib > 1 ) THEN 
     498                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN 
    869499                     CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', & 
    870                           &        ' in order of distance from edge nbr A utility for re-ordering ', & 
    871                           &        ' boundary coordinates and data files exists in the TOOLS/OBC directory') 
    872                   ENDIF     
     500                        &        ' in order of distance from edge nbr A utility for re-ordering ', & 
     501                        &        ' boundary coordinates and data files exists in the TOOLS/OBC directory') 
     502                  ENDIF 
    873503               ENDIF 
    874504               ! check if point is in local domain 
     
    876506                  & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino      ) THEN 
    877507                  ! 
    878                   icount = icount  + 1 
    879                   ! 
    880                   IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1 
     508                  icount = icount + 1 
     509                  IF( nbrdta(ib,igrd,ib_bdy) == 1 .OR. nbrdta(ib,igrd,ib_bdy) == 0 )   icountr = icountr + 1 
     510                  IF( nbrdta(ib,igrd,ib_bdy) == 0 )   icountr0 = icountr0 + 1 
    881511               ENDIF 
    882512            END DO 
    883             idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc 
    884             idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc         
    885          END DO  ! igrd 
     513            idx_bdy(ib_bdy)%nblen    (igrd) = icount   !: length of boundary data on each proc 
     514            idx_bdy(ib_bdy)%nblenrim (igrd) = icountr  !: length of rim 0 and rim 1 boundary data on each proc    
     515            idx_bdy(ib_bdy)%nblenrim0(igrd) = icountr0 !: length of rim 0 boundary data on each proc      
     516         END DO   ! igrd 
    886517 
    887518         ! Allocate index arrays for this boundary set 
     
    893524            &      idx_bdy(ib_bdy)%nbd   (ilen1,jpbgrd) ,   & 
    894525            &      idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) ,   & 
     526            &      idx_bdy(ib_bdy)%ntreat(ilen1,jpbgrd) ,   & 
    895527            &      idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) ,   & 
    896528            &      idx_bdy(ib_bdy)%nbw   (ilen1,jpbgrd) ,   & 
     
    900532         ! Dispatch mapping indices and discrete distances on each processor 
    901533         ! ----------------------------------------------------------------- 
    902  
    903          com_east  = 0 
    904          com_west  = 0 
    905          com_south = 0 
    906          com_north = 0 
    907  
    908          com_east_b  = 0 
    909          com_west_b  = 0 
    910          com_south_b = 0 
    911          com_north_b = 0 
    912  
    913534         DO igrd = 1, jpbgrd 
    914535            icount  = 0 
    915             ! Loop on rimwidth to ensure outermost points come first in the local arrays. 
    916             DO ir=1, nn_rimwidth(ib_bdy) 
     536            ! Outer loop on rimwidth to ensure outermost points come first in the local arrays. 
     537            DO ir = 0, nn_rimwidth(ib_bdy) 
    917538               DO ib = 1, nblendta(igrd,ib_bdy) 
    918539                  ! check if point is in local domain and equals ir 
     
    922543                     ! 
    923544                     icount = icount  + 1 
    924  
    925                      ! Rather assume that boundary data indices are given on global domain 
    926                      ! TO BE DISCUSSED ? 
    927 !                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
    928 !                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
    929                      idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
    930                      idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
    931                      ! check if point has to be sent 
    932                      ii = idx_bdy(ib_bdy)%nbi(icount,igrd) 
    933                      ij = idx_bdy(ib_bdy)%nbj(icount,igrd) 
    934                      if((com_east .ne. 1) .and. (ii == (nlci-1)) .and. (nbondi .le. 0)) then 
    935                         com_east = 1 
    936                      elseif((com_west .ne. 1) .and. (ii == 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then 
    937                         com_west = 1 
    938                      endif  
    939                      if((com_south .ne. 1) .and. (ij == 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then 
    940                         com_south = 1 
    941                      elseif((com_north .ne. 1) .and. (ij == (nlcj-1)) .and. (nbondj .le. 0)) then 
    942                         com_north = 1 
    943                      endif  
     545                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1   ! global to local indexes 
     546                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1   ! global to local indexes 
    944547                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy) 
    945548                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 
    946549                  ENDIF 
    947                   ! check if point has to be received from a neighbour 
    948                   IF(nbondi == 0) THEN 
    949                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   & 
    950                        & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   & 
    951                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    952                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2 
    953                        if( ii == (nlcit(nowe+1)-1) ) then 
    954                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2 
    955                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 
    956                             com_south = 1 
    957                           elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then 
    958                             com_north = 1 
    959                           endif 
    960                           com_west_b = 1 
    961                        endif  
    962                      ENDIF 
    963                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   & 
    964                        & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   & 
    965                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    966                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2 
    967                        if( ii == 2 ) then 
    968                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2 
    969                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 
    970                             com_south = 1 
    971                           elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then 
    972                             com_north = 1 
    973                           endif 
    974                           com_east_b = 1 
    975                        endif  
    976                      ENDIF 
    977                   ELSEIF(nbondi == 1) THEN 
    978                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   & 
    979                        & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   & 
    980                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    981                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2 
    982                        if( ii == (nlcit(nowe+1)-1) ) then 
    983                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2 
    984                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 
    985                             com_south = 1 
    986                           elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then 
    987                             com_north = 1 
    988                           endif 
    989                           com_west_b = 1 
    990                        endif  
    991                      ENDIF 
    992                   ELSEIF(nbondi == -1) THEN 
    993                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   & 
    994                        & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   & 
    995                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    996                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2 
    997                        if( ii == 2 ) then 
    998                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2 
    999                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then 
    1000                             com_south = 1 
    1001                           elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then 
    1002                             com_north = 1 
    1003                           endif 
    1004                           com_east_b = 1 
    1005                        endif  
    1006                      ENDIF 
    1007                   ENDIF 
    1008                   IF(nbondj == 0) THEN 
    1009                      IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  & 
    1010                        & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. & 
    1011                        & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 
    1012                        com_north_b = 1  
    1013                      ENDIF 
    1014                      IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1  & 
    1015                        &.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. & 
    1016                        & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 
    1017                        com_south_b = 1  
    1018                      ENDIF 
    1019                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   & 
    1020                        & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   & 
    1021                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    1022                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2 
    1023                        if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then 
    1024                           com_south_b = 1 
    1025                        endif  
    1026                      ENDIF 
    1027                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   & 
    1028                        & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   & 
    1029                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    1030                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2 
    1031                        if((com_north_b .ne. 1) .and. (ij == 2)) then 
    1032                           com_north_b = 1 
    1033                        endif  
    1034                      ENDIF 
    1035                   ELSEIF(nbondj == 1) THEN 
    1036                      IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. & 
    1037                        & nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. & 
    1038                        & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 
    1039                        com_south_b = 1  
    1040                      ENDIF 
    1041                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   & 
    1042                        & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   & 
    1043                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    1044                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2 
    1045                        if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then 
    1046                           com_south_b = 1 
    1047                        endif  
    1048                      ENDIF 
    1049                   ELSEIF(nbondj == -1) THEN 
    1050                      IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  & 
    1051                        & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. & 
    1052                        & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 
    1053                        com_north_b = 1  
    1054                      ENDIF 
    1055                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   & 
    1056                        & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   & 
    1057                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    1058                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2 
    1059                        if((com_north_b .ne. 1) .and. (ij == 2)) then 
    1060                           com_north_b = 1 
    1061                        endif  
    1062                      ENDIF 
    1063                   ENDIF 
    1064                ENDDO 
    1065             ENDDO 
    1066          ENDDO  
    1067  
    1068          ! definition of the i- and j- direction local boundaries arrays used for sending the boundaries 
    1069          IF(     (com_east  == 1) .and. (com_west  == 1) ) THEN   ;   nbondi_bdy(ib_bdy) =  0 
    1070          ELSEIF( (com_east  == 1) .and. (com_west  == 0) ) THEN   ;   nbondi_bdy(ib_bdy) = -1 
    1071          ELSEIF( (com_east  == 0) .and. (com_west  == 1) ) THEN   ;   nbondi_bdy(ib_bdy) =  1 
    1072          ENDIF 
    1073          IF(     (com_north == 1) .and. (com_south == 1) ) THEN   ;   nbondj_bdy(ib_bdy) =  0 
    1074          ELSEIF( (com_north == 1) .and. (com_south == 0) ) THEN   ;   nbondj_bdy(ib_bdy) = -1 
    1075          ELSEIF( (com_north == 0) .and. (com_south == 1) ) THEN   ;   nbondj_bdy(ib_bdy) =  1 
    1076          ENDIF 
    1077  
    1078          ! definition of the i- and j- direction local boundaries arrays used for receiving the boundaries 
    1079          IF(     (com_east_b  == 1) .and. (com_west_b  == 1) ) THEN   ;   nbondi_bdy_b(ib_bdy) =  0 
    1080          ELSEIF( (com_east_b  == 1) .and. (com_west_b  == 0) ) THEN   ;   nbondi_bdy_b(ib_bdy) = -1 
    1081          ELSEIF( (com_east_b  == 0) .and. (com_west_b  == 1) ) THEN   ;   nbondi_bdy_b(ib_bdy) =  1 
    1082          ENDIF 
    1083          IF(     (com_north_b == 1) .and. (com_south_b == 1) ) THEN   ;   nbondj_bdy_b(ib_bdy) =  0 
    1084          ELSEIF( (com_north_b == 1) .and. (com_south_b == 0) ) THEN   ;   nbondj_bdy_b(ib_bdy) = -1 
    1085          ELSEIF( (com_north_b == 0) .and. (com_south_b == 1) ) THEN   ;   nbondj_bdy_b(ib_bdy) =  1 
    1086          ENDIF 
     550               END DO 
     551            END DO 
     552         END DO   ! igrd 
     553 
     554      END DO   ! ib_bdy 
     555 
     556      ! Initialize array indicating communications in bdy 
     557      ! ------------------------------------------------- 
     558      ALLOCATE( lsend_bdy(nb_bdy,jpbgrd,4,0:1), lrecv_bdy(nb_bdy,jpbgrd,4,0:1) ) 
     559      lsend_bdy(:,:,:,:) = .false. 
     560      lrecv_bdy(:,:,:,:) = .false.  
     561 
     562      DO ib_bdy = 1, nb_bdy 
     563         DO igrd = 1, jpbgrd 
     564            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! only the rim triggers communications, see bdy routines 
     565               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     566               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     567               IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN   ;   ir = 0 
     568               ELSE                                                 ;   ir = 1 
     569               END IF 
     570               ! 
     571               ! check if point has to be sent     to   a neighbour 
     572               ! W neighbour and on the inner left  side 
     573               IF( ii == 2     .and. (nbondi == 0 .or. nbondi ==  1) )   lsend_bdy(ib_bdy,igrd,1,ir) = .true. 
     574               ! E neighbour and on the inner right side 
     575               IF( ii == jpi-1 .and. (nbondi == 0 .or. nbondi == -1) )   lsend_bdy(ib_bdy,igrd,2,ir) = .true. 
     576               ! S neighbour and on the inner down side 
     577               IF( ij == 2     .and. (nbondj == 0 .or. nbondj ==  1) )   lsend_bdy(ib_bdy,igrd,3,ir) = .true. 
     578               ! N neighbour and on the inner up   side 
     579               IF( ij == jpj-1 .and. (nbondj == 0 .or. nbondj == -1) )   lsend_bdy(ib_bdy,igrd,4,ir) = .true. 
     580               ! 
     581               ! check if point has to be received from a neighbour 
     582               ! W neighbour and on the outter left  side 
     583               IF( ii == 1     .and. (nbondi == 0 .or. nbondi ==  1) )   lrecv_bdy(ib_bdy,igrd,1,ir) = .true. 
     584               ! E neighbour and on the outter right side 
     585               IF( ii == jpi   .and. (nbondi == 0 .or. nbondi == -1) )   lrecv_bdy(ib_bdy,igrd,2,ir) = .true. 
     586               ! S neighbour and on the outter down side 
     587               IF( ij == 1     .and. (nbondj == 0 .or. nbondj ==  1) )   lrecv_bdy(ib_bdy,igrd,3,ir) = .true. 
     588               ! N neighbour and on the outter up   side 
     589               IF( ij == jpj   .and. (nbondj == 0 .or. nbondj == -1) )   lrecv_bdy(ib_bdy,igrd,4,ir) = .true. 
     590               ! 
     591            END DO 
     592         END DO  ! igrd 
    1087593 
    1088594         ! Compute rim weights for FRS scheme 
     
    1090596         DO igrd = 1, jpbgrd 
    1091597            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    1092                nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    1093                idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( nbr - 1 ) *0.5 )      ! tanh formulation 
    1094 !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
    1095 !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy))       ! linear 
    1096             END DO 
    1097          END DO  
     598               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same weights 
     599               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation 
     600               !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
     601               !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear 
     602            END DO 
     603         END DO 
    1098604 
    1099605         ! Compute damping coefficients 
     
    1101607         DO igrd = 1, jpbgrd 
    1102608            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    1103                nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
     609               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same damping coefficients 
    1104610               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &  
    1105                & *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
     611                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
    1106612               idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) &  
    1107                & *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
    1108             END DO 
    1109          END DO  
    1110  
    1111       END DO 
     613                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
     614            END DO 
     615         END DO 
     616 
     617      END DO ! ib_bdy 
    1112618 
    1113619      ! ------------------------------------------------------ 
    1114620      ! Initialise masks and find normal/tangential directions 
    1115621      ! ------------------------------------------------------ 
     622 
     623      ! ------------------------------------------ 
     624      ! handle rim0, do as if rim 1 was free ocean 
     625      ! ------------------------------------------ 
     626 
     627      ztmask(:,:) = tmask(:,:,1)   ;   zumask(:,:) = umask(:,:,1)   ;   zvmask(:,:) = vmask(:,:,1) 
     628      ! For the flagu/flagv calculation below we require a version of fmask without 
     629      ! the land boundary condition (shlat) included: 
     630      DO ij = 1, jpjm1 
     631         DO ii = 1, jpim1 
     632            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   & 
     633               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
     634         END DO 
     635      END DO 
     636      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) 
    1116637 
    1117638      ! Read global 2D mask at T-points: bdytmask 
     
    1119640      ! bdytmask = 1  on the computational domain AND on open boundaries 
    1120641      !          = 0  elsewhere    
    1121   
     642 
    1122643      bdytmask(:,:) = ssmask(:,:) 
    1123644 
    1124645      ! Derive mask on U and V grid from mask on T grid 
    1125  
    1126       bdyumask(:,:) = 0._wp 
    1127       bdyvmask(:,:) = 0._wp 
    1128646      DO ij = 1, jpjm1 
    1129647         DO ii = 1, jpim1 
    1130             bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij ) 
     648            bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1,ij ) 
    1131649            bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1)   
    1132650         END DO 
    1133651      END DO 
    1134       CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1. , bdyvmask, 'V', 1. )   ! Lateral boundary cond.  
    1135  
    1136       ! bdy masks are now set to zero on boundary points: 
    1137       ! 
    1138       igrd = 1 
     652      CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1., bdyvmask, 'V', 1. )   ! Lateral boundary cond. 
     653 
     654      ! bdy masks are now set to zero on rim 0 points: 
    1139655      DO ib_bdy = 1, nb_bdy 
    1140         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)       
    1141           bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 
    1142         END DO 
    1143       END DO 
    1144       ! 
    1145       igrd = 2 
     656         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0 
     657            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 
     658         END DO 
     659         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0 
     660            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 
     661         END DO 
     662         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0 
     663            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 
     664         END DO 
     665      END DO 
     666 
     667      CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. )   ! compute flagu, flagv, ntreat on rim 0 
     668 
     669      ! ------------------------------------ 
     670      ! handle rim1, do as if rim 0 was land 
     671      ! ------------------------------------ 
     672       
     673      ! z[tuv]mask are now set to zero on rim 0 points: 
    1146674      DO ib_bdy = 1, nb_bdy 
    1147         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    1148           bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 
    1149         END DO 
    1150       END DO 
    1151       ! 
    1152       igrd = 3 
     675         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0 
     676            ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 
     677         END DO 
     678         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0 
     679            zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 
     680         END DO 
     681         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0 
     682            zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 
     683         END DO 
     684      END DO 
     685 
     686      ! Recompute zfmask 
     687      DO ij = 1, jpjm1 
     688         DO ii = 1, jpim1 
     689            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   & 
     690               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 
     691         END DO 
     692      END DO 
     693      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) 
     694 
     695      ! bdy masks are now set to zero on rim1 points: 
    1153696      DO ib_bdy = 1, nb_bdy 
    1154         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    1155           bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 
    1156         END DO 
    1157       END DO 
    1158  
    1159       ! For the flagu/flagv calculation below we require a version of fmask without 
    1160       ! the land boundary condition (shlat) included: 
    1161       zfmask(:,:) = 0 
    1162       DO ij = 2, jpjm1 
    1163          DO ii = 2, jpim1 
    1164             zfmask(ii,ij) = tmask(ii,ij  ,1) * tmask(ii+1,ij  ,1)   & 
    1165            &              * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1) 
    1166          END DO       
    1167       END DO 
    1168  
    1169       ! Lateral boundary conditions 
    1170       CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )  
    1171       CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1. , bdyvmask, 'V', 1., bdytmask, 'T', 1. ) 
     697         DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1,  idx_bdy(ib_bdy)%nblenrim(1)   ! extent of rim 1 
     698            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 
     699         END DO 
     700         DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1,  idx_bdy(ib_bdy)%nblenrim(2)   ! extent of rim 1 
     701            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 
     702         END DO 
     703         DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1,  idx_bdy(ib_bdy)%nblenrim(3)   ! extent of rim 1 
     704            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 
     705         END DO 
     706      END DO 
     707 
     708      CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1 
     709      ! 
     710      ! Check which boundaries might need communication 
     711      ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,4,0:1), lrecv_bdyint(nb_bdy,jpbgrd,4,0:1) ) 
     712      lsend_bdyint(:,:,:,:) = .false. 
     713      lrecv_bdyint(:,:,:,:) = .false.  
     714      ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,4,0:1), lrecv_bdyext(nb_bdy,jpbgrd,4,0:1) ) 
     715      lsend_bdyext(:,:,:,:) = .false. 
     716      lrecv_bdyext(:,:,:,:) = .false. 
     717      ! 
     718      DO igrd = 1, jpbgrd 
     719         DO ib_bdy = 1, nb_bdy 
     720            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     721               IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE 
     722               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     723               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     724               ir = idx_bdy(ib_bdy)%nbr(ib,igrd) 
     725               flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd)) 
     726               flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd)) 
     727               iibe = ii - flagu   ! neighbouring point towards the exterior of the computational domain 
     728               ijbe = ij - flagv 
     729               iibi = ii + flagu   ! neighbouring point towards the interior of the computational domain 
     730               ijbi = ij + flagv 
     731               CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 )   ! free ocean neighbours 
     732               ! 
     733               ! search neighbour in the  west/east  direction 
     734               ! Rim is on the halo and computed ocean is towards exterior of mpi domain   
     735               !      <--    (o exterior)     -->   
     736               ! (1)  o|x         OR    (2)   x|o 
     737               !       |___                 ___|  
     738               IF( iibi == 0     .OR. ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1,ir) = .true. 
     739               IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2,ir) = .true.   
     740               IF( iibe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,1,ir) = .true. 
     741               IF( iibe == jpi+1                                                       )   lrecv_bdyext(ib_bdy,igrd,2,ir) = .true.   
     742               ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo 
     743               ! :¨¨¨¨¨|¨¨-->    |                                             |    <--¨¨|¨¨¨¨¨:  
     744               ! :     |  x:o    |    neighbour limited by ... would need o    |    o:x  |     : 
     745               ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....: 
     746               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. & 
     747                  & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1,ir)=.true. 
     748               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. & 
     749                  & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2,ir)=.true. 
     750               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3     )   lsend_bdyext(ib_bdy,igrd,1,ir)=.true. 
     751               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 )   lsend_bdyext(ib_bdy,igrd,2,ir)=.true. 
     752               ! 
     753               ! search neighbour in the north/south direction    
     754               ! Rim is on the halo and computed ocean is towards exterior of mpi domain 
     755               !(3)   |       |         ^   ___o___      
     756               !  |   |___x___|   OR    |  |   x   | 
     757               !  v       o           (4)  |       | 
     758               IF( ijbi == 0     .OR. ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3,ir) = .true. 
     759               IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4,ir) = .true. 
     760               IF( ijbe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,3,ir) = .true. 
     761               IF( ijbe == jpj+1                                                       )   lrecv_bdyext(ib_bdy,igrd,4,ir) = .true. 
     762               ! Check if neighbour has its rim parallel to its mpi subdomain     _________  border and next to its halo 
     763               !   ^  |    o    |                                                :         :  
     764               !   |  |¨¨¨¨x¨¨¨¨|   neighbour limited by ... would need o     |  |....x....| 
     765               !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |    
     766               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. & 
     767                  & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3,ir)=.true. 
     768               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. & 
     769                  & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4,ir)=.true. 
     770               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3     )   lsend_bdyext(ib_bdy,igrd,3,ir)=.true. 
     771               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 )   lsend_bdyext(ib_bdy,igrd,4,ir)=.true. 
     772            END DO 
     773         END DO 
     774      END DO 
     775 
     776      DO ib_bdy = 1,nb_bdy 
     777         IF(  cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 
     778            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 
     779            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN 
     780            DO igrd = 1, jpbgrd 
     781               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     782                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     783                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     784                  IF(  mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2  ) THEN 
     785                     WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain' 
     786                     CALL ctl_stop( ctmp1 ) 
     787                  END IF 
     788               END DO 
     789            END DO 
     790         END IF 
     791      END DO 
     792      ! 
     793      DEALLOCATE( nbidta, nbjdta, nbrdta ) 
     794      ! 
     795   END SUBROUTINE bdy_def 
     796 
     797 
     798   SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 ) 
     799      !!---------------------------------------------------------------------- 
     800      !!                 ***  ROUTINE bdy_rim_treat  *** 
     801      !! 
     802      !! ** Purpose :   Initialize structures ( flagu, flagv, ntreat ) indicating how rim points 
     803      !!                  are to be handled in the boundary condition treatment 
     804      !! 
     805      !! ** Method  :   - to handle rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior) 
     806      !!                            and bdymasks must be set at 0 on rim 0       (set at one on rim 1 and interior) 
     807      !!                    (as if rim 1 was free ocean) 
     808      !!                - to handle rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior) 
     809      !!                            and bdymasks must indicate free ocean points (set at one on interior) 
     810      !!                    (as if rim 0 was land) 
     811      !!                - we can then check in which direction the interior of the computational domain is with the difference 
     812      !!                         mask array values on both sides to compute flagu and flagv 
     813      !!                - and look at the ocean neighbours to compute ntreat 
     814      !!---------------------------------------------------------------------- 
     815      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask   ! temporary fmask excluding coastal boundary condition (shlat) 
     816      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary t/u/v mask array 
     817      LOGICAL                             , INTENT (in   ) :: lrim0    ! .true. -> rim 0   .false. -> rim 1 
     818      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices 
     819      INTEGER  ::   i_offset, j_offset, inn                ! local integer 
     820      INTEGER  ::   ibeg, iend                             ! local integer 
     821      LOGICAL  ::   llnon, llson, llean, llwen             ! local logicals indicating the presence of a ocean neighbour 
     822      REAL(wp), POINTER, DIMENSION(:,:)       ::   zmask   ! pointer to 2D mask fields 
     823      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
     824      CHARACTER(LEN=1), DIMENSION(jpbgrd)     ::   cgrid 
     825      REAL(wp)        , DIMENSION(jpi,jpj)    ::   ztmp 
     826      !!---------------------------------------------------------------------- 
     827 
     828      cgrid = (/'t','u','v'/) 
     829 
    1172830      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components 
    1173  
    1174          idx_bdy(ib_bdy)%flagu(:,:) = 0._wp 
    1175          idx_bdy(ib_bdy)%flagv(:,:) = 0._wp 
    1176          icount = 0  
    1177831 
    1178832         ! Calculate relationship of U direction to the local orientation of the boundary 
     
    1180834         ! flagu =  0 : u is tangential 
    1181835         ! flagu =  1 : u is normal to the boundary and is direction is inward 
    1182    
    1183836         DO igrd = 1, jpbgrd  
    1184837            SELECT CASE( igrd ) 
    1185                CASE( 1 )   ;   pmask => umask   (:,:,1)   ;   i_offset = 0 
    1186                CASE( 2 )   ;   pmask => bdytmask(:,:)     ;   i_offset = 1 
    1187                CASE( 3 )   ;   pmask => zfmask  (:,:)     ;   i_offset = 0 
     838               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0 
     839               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1 
     840               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0 
    1188841            END SELECT  
    1189842            icount = 0 
    1190             DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
    1191                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    1192                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1193                zefl = pmask(nbi+i_offset-1,nbj) 
    1194                zwfl = pmask(nbi+i_offset,nbj) 
     843            ztmp(:,:) = -999._wp 
     844            IF( lrim0 ) THEN   ! extent of rim 0 
     845               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 
     846            ELSE               ! extent of rim 1 
     847               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd) 
     848            END IF 
     849            DO ib = ibeg, iend  
     850               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     851               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     852               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE 
     853               zwfl = zmask(ii+i_offset-1,ij) 
     854               zefl = zmask(ii+i_offset  ,ij) 
    1195855               ! This error check only works if you are using the bdyXmask arrays 
    1196                IF( i_offset == 1 .and. zefl + zwfl == 2 ) THEN 
     856               IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN 
    1197857                  icount = icount + 1 
    1198                   IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj) 
     858                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij) 
    1199859               ELSE 
    1200                   idx_bdy(ib_bdy)%flagu(ib,igrd) = -zefl + zwfl 
     860                  ztmp(ii,ij) = -zwfl + zefl 
    1201861               ENDIF 
    1202862            END DO 
    1203863            IF( icount /= 0 ) THEN 
    1204                WRITE(ctmp1,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   & 
     864               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   & 
    1205865                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 
    1206                WRITE(ctmp2,*) ' ========== ' 
    1207                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     866               CALL ctl_stop( ctmp1 ) 
    1208867            ENDIF  
     868            SELECT CASE( igrd ) 
     869               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )  
     870               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )  
     871               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )  
     872            END SELECT  
     873            DO ib = ibeg, iend 
     874               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     875               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     876               idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij) 
     877            END DO 
    1209878         END DO 
    1210879 
     
    1213882         ! flagv =  0 : v is tangential 
    1214883         ! flagv =  1 : v is normal to the boundary and is direction is inward 
    1215  
    1216884         DO igrd = 1, jpbgrd  
    1217885            SELECT CASE( igrd ) 
    1218                CASE( 1 )   ;   pmask => vmask (:,:,1)   ;   j_offset = 0 
    1219                CASE( 2 )   ;   pmask => zfmask(:,:)     ;   j_offset = 0 
    1220                CASE( 3 )   ;   pmask => bdytmask        ;   j_offset = 1 
     886               CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0 
     887               CASE( 2 )   ;   zmask => pfmask     ;   j_offset = 0 
     888               CASE( 3 )   ;   zmask => bdytmask   ;   j_offset = 1 
    1221889            END SELECT  
    1222890            icount = 0 
    1223             DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   
    1224                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    1225                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1226                znfl = pmask(nbi,nbj+j_offset-1) 
    1227                zsfl = pmask(nbi,nbj+j_offset  ) 
     891            ztmp(:,:) = -999._wp 
     892            IF( lrim0 ) THEN   ! extent of rim 0 
     893               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 
     894            ELSE               ! extent of rim 1 
     895               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd) 
     896            END IF 
     897            DO ib = ibeg, iend 
     898               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     899               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     900               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE 
     901               zsfl = zmask(ii,ij+j_offset-1) 
     902               znfl = zmask(ii,ij+j_offset  ) 
    1228903               ! This error check only works if you are using the bdyXmask arrays 
    1229                IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN 
    1230                   IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj) 
     904               IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN 
     905                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij) 
    1231906                  icount = icount + 1 
    1232907               ELSE 
    1233                   idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl 
     908                  ztmp(ii,ij) = -zsfl + znfl 
    1234909               END IF 
    1235910            END DO 
    1236911            IF( icount /= 0 ) THEN 
    1237                WRITE(ctmp1,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   & 
     912               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   & 
    1238913                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 
    1239                WRITE(ctmp2,*) ' ========== ' 
    1240                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
    1241             ENDIF  
    1242          END DO 
    1243          ! 
    1244       END DO 
    1245       ! 
    1246       ! Tidy up 
    1247       !-------- 
    1248       IF( nb_bdy>0 )   DEALLOCATE( nbidta, nbjdta, nbrdta ) 
    1249       ! 
    1250    END SUBROUTINE bdy_segs 
    1251  
     914               CALL ctl_stop( ctmp1 ) 
     915            ENDIF 
     916            SELECT CASE( igrd ) 
     917               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )  
     918               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )  
     919               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )  
     920            END SELECT  
     921            DO ib = ibeg, iend 
     922               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     923               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     924               idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij) 
     925            END DO 
     926         END DO 
     927         ! 
     928      END DO ! ib_bdy 
     929       
     930      DO ib_bdy = 1, nb_bdy 
     931         DO igrd = 1, jpbgrd 
     932            SELECT CASE( igrd ) 
     933               CASE( 1 )   ;   zmask => bdytmask  
     934               CASE( 2 )   ;   zmask => bdyumask  
     935               CASE( 3 )   ;   zmask => bdyvmask  
     936            END SELECT 
     937            ztmp(:,:) = -999._wp 
     938            IF( lrim0 ) THEN   ! extent of rim 0 
     939               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 
     940            ELSE               ! extent of rim 1 
     941               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd) 
     942            END IF 
     943            DO ib = ibeg, iend 
     944               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     945               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     946               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE 
     947               llnon = zmask(ii  ,ij+1) == 1.   
     948               llson = zmask(ii  ,ij-1) == 1.  
     949               llean = zmask(ii+1,ij  ) == 1.  
     950               llwen = zmask(ii-1,ij  ) == 1.  
     951               inn  = COUNT( (/ llnon, llson, llean, llwen /) ) 
     952               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points 
     953                  !               !              !     _____     !     _____    !    __     __ 
     954                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error 
     955                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x| 
     956                  IF(     zmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1. 
     957                  ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2. 
     958                  ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3. 
     959                  ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4. 
     960                  ELSE                                    ;   ztmp(ii,ij) = -1. 
     961                     WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   & 
     962                       ' on boundary set ', ib_bdy, ' has no free ocean neighbour' 
     963                     IF( lrim0 ) THEN 
     964                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.' 
     965                     ELSE 
     966                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.' 
     967                     END IF 
     968                     CALL ctl_warn( ctmp1, ctmp2 ) 
     969                  END IF 
     970               END IF 
     971               IF( inn == 1 ) THEN   ! middle of linear bdy  or incomplete corner  ! ___ o 
     972                  !    |         !         |   !      o     !    ______            !    |x___ 
     973                  ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x 
     974                  !    |         !         |   !            !      o 
     975                  IF( llean )   ztmp(ii,ij) = 5. 
     976                  IF( llwen )   ztmp(ii,ij) = 6. 
     977                  IF( llnon )   ztmp(ii,ij) = 7. 
     978                  IF( llson )   ztmp(ii,ij) = 8. 
     979               END IF 
     980               IF( inn == 2 ) THEN   ! exterior of a corner 
     981                  !        o      !        o      !    _____|       !       |_____   
     982                  !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x       
     983                  !         |     !       |       !        o        !        o  
     984                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9. 
     985                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10. 
     986                  IF( llson .AND. llean )   ztmp(ii,ij) = 11. 
     987                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12. 
     988               END IF 
     989               IF( inn == 3 ) THEN   ! 3 neighbours     __   __ 
     990                  !    |_  o      !        o  _|  !       |_|     !       o          
     991                  ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o        
     992                  !    |   o      !        o   |  !        o      !    __|¨|__     
     993                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13. 
     994                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14. 
     995                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15. 
     996                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16. 
     997               END IF 
     998               IF( inn == 4 ) THEN 
     999                  WRITE(ctmp1,*)  'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   & 
     1000                       ' on boundary set ', ib_bdy, ' have 4 neighbours' 
     1001                  CALL ctl_stop( ctmp1 ) 
     1002               END IF 
     1003            END DO 
     1004            SELECT CASE( igrd ) 
     1005               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )  
     1006               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. )  
     1007               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. )  
     1008            END SELECT  
     1009            DO ib = ibeg, iend 
     1010               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1011               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1012               idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij)) 
     1013            END DO 
     1014         END DO 
     1015      END DO 
     1016 
     1017    END SUBROUTINE bdy_rim_treat 
     1018 
     1019    
     1020    SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 ) 
     1021      !!---------------------------------------------------------------------- 
     1022      !!                 ***  ROUTINE find_neib  *** 
     1023      !! 
     1024      !! ** Purpose :   get ii1, ij1, ii2, ij2, ii3, ij3, the indices of 
     1025      !!               the free ocean neighbours of (ii,ij) for bdy treatment 
     1026      !! 
     1027      !! ** Method  :  use itreat input to select a case 
     1028      !!               N.B. ntreat is defined for all bdy points in routine bdy_rim_treat 
     1029      !! 
     1030      !!---------------------------------------------------------------------- 
     1031      INTEGER, INTENT(in   )      ::   ii, ij, itreat 
     1032      INTEGER, INTENT(  out)      ::   ii1, ij1, ii2, ij2, ii3, ij3 
     1033      !!---------------------------------------------------------------------- 
     1034      SELECT CASE( itreat )   ! points that will be used by bdy routines, -1 will be discarded 
     1035         !               !               !     _____     !     _____      
     1036         !  1 |   o      !  2  o   |     !  3 | x        !  4     x |     
     1037         !    |_x_ _     !    _ _x_|     !    |   o      !      o   | 
     1038      CASE( 1 )    ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1039      CASE( 2 )    ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1040      CASE( 3 )    ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1041      CASE( 4 )    ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1042         !    |          !         |     !      o        !    ______                   ! or incomplete corner 
     1043         ! 5  | x o      ! 6   o x |     ! 7  __x__      ! 8    x                      !  7  ____ o 
     1044         !    |          !         |     !               !      o                      !         |x___ 
     1045      CASE( 5 )    ;   ii1 = ii+1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1046      CASE( 6 )    ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1047      CASE( 7 )    ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1048      CASE( 8 )    ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1 
     1049         !        o      !        o      !    _____|     !       |_____   
     1050         !  9 ____x o    ! 10   o x___   ! 11     x o    ! 12   o x       
     1051         !         |     !       |       !        o      !        o       
     1052      CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1  
     1053      CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
     1054      CASE( 11 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
     1055      CASE( 12 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
     1056         !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o          
     1057         ! 13  _| x o    !  14  o x |_   !  15  o x o    ! 16  o x o        
     1058         !    |   o      !        o   |  !        o      !    __|¨|__  
     1059      CASE( 13 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1    
     1060      CASE( 14 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1  
     1061      CASE( 15 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij    
     1062      CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij 
     1063      END SELECT 
     1064   END SUBROUTINE find_neib 
     1065     
     1066 
     1067   SUBROUTINE bdy_read_seg( kb_bdy, knblendta )  
     1068      !!---------------------------------------------------------------------- 
     1069      !!                 ***  ROUTINE bdy_coords_seg  *** 
     1070      !! 
     1071      !! ** Purpose :  build bdy coordinates with segments defined in namelist 
     1072      !! 
     1073      !! ** Method  :  read namelist nambdy_index blocks 
     1074      !! 
     1075      !!---------------------------------------------------------------------- 
     1076      INTEGER                   , INTENT (in   ) ::   kb_bdy           ! bdy number 
     1077      INTEGER, DIMENSION(jpbgrd), INTENT (  out) ::   knblendta        ! length of index arrays  
     1078      !! 
     1079      INTEGER          ::   ios                 ! Local integer output status for namelist read 
     1080      INTEGER          ::   nbdyind, nbdybeg, nbdyend 
     1081      CHARACTER(LEN=1) ::   ctypebdy   !     -        -  
     1082      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
     1083      !!---------------------------------------------------------------------- 
     1084 
     1085      ! No REWIND here because may need to read more than one nambdy_index namelist. 
     1086      ! Read only namelist_cfg to avoid unseccessfull overwrite  
     1087      ! keep full control of the configuration namelist 
     1088      READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 ) 
     1089904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' ) 
     1090      IF(lwm) WRITE ( numond, nambdy_index ) 
     1091       
     1092      SELECT CASE ( TRIM(ctypebdy) ) 
     1093      CASE( 'N' ) 
     1094         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     1095            nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain. 
     1096            nbdybeg  = 2 
     1097            nbdyend  = jpiglo - 1 
     1098         ENDIF 
     1099         nbdysegn = nbdysegn + 1 
     1100         npckgn(nbdysegn) = kb_bdy ! Save bdy package number 
     1101         jpjnob(nbdysegn) = nbdyind 
     1102         jpindt(nbdysegn) = nbdybeg 
     1103         jpinft(nbdysegn) = nbdyend 
     1104         ! 
     1105      CASE( 'S' ) 
     1106         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     1107            nbdyind  = 2           ! set boundary to whole side of model domain. 
     1108            nbdybeg  = 2 
     1109            nbdyend  = jpiglo - 1 
     1110         ENDIF 
     1111         nbdysegs = nbdysegs + 1 
     1112         npckgs(nbdysegs) = kb_bdy ! Save bdy package number 
     1113         jpjsob(nbdysegs) = nbdyind 
     1114         jpisdt(nbdysegs) = nbdybeg 
     1115         jpisft(nbdysegs) = nbdyend 
     1116         ! 
     1117      CASE( 'E' ) 
     1118         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     1119            nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain. 
     1120            nbdybeg  = 2 
     1121            nbdyend  = jpjglo - 1 
     1122         ENDIF 
     1123         nbdysege = nbdysege + 1  
     1124         npckge(nbdysege) = kb_bdy ! Save bdy package number 
     1125         jpieob(nbdysege) = nbdyind 
     1126         jpjedt(nbdysege) = nbdybeg 
     1127         jpjeft(nbdysege) = nbdyend 
     1128         ! 
     1129      CASE( 'W' ) 
     1130         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     1131            nbdyind  = 2           ! set boundary to whole side of model domain. 
     1132            nbdybeg  = 2 
     1133            nbdyend  = jpjglo - 1 
     1134         ENDIF 
     1135         nbdysegw = nbdysegw + 1 
     1136         npckgw(nbdysegw) = kb_bdy ! Save bdy package number 
     1137         jpiwob(nbdysegw) = nbdyind 
     1138         jpjwdt(nbdysegw) = nbdybeg 
     1139         jpjwft(nbdysegw) = nbdyend 
     1140         ! 
     1141      CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 
     1142      END SELECT 
     1143       
     1144      ! For simplicity we assume that in case of straight bdy, arrays have the same length 
     1145      ! (even if it is true that last tangential velocity points 
     1146      ! are useless). This simplifies a little bit boundary data format (and agrees with format 
     1147      ! used so far in obc package) 
     1148       
     1149      knblendta(1:jpbgrd) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy) 
     1150       
     1151   END SUBROUTINE bdy_read_seg 
     1152 
     1153    
    12521154   SUBROUTINE bdy_ctl_seg 
    12531155      !!---------------------------------------------------------------------- 
     
    12791181            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
    12801182         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
    1281          IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
    1282          IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
     1183         IF (jpindt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     1184         IF (jpinft(ib).gt.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
    12831185      END DO 
    12841186      ! 
     
    12881190            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
    12891191         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
    1290          IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
    1291          IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
     1192         IF (jpisdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     1193         IF (jpisft(ib).gt.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
    12921194      END DO 
    12931195      ! 
     
    12971199            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
    12981200         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
    1299          IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
    1300          IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
     1201         IF (jpjedt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     1202         IF (jpjeft(ib).gt.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
    13011203      END DO 
    13021204      ! 
     
    13061208            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
    13071209         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
    1308          IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
    1309          IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
     1210         IF (jpjwdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     1211         IF (jpjwft(ib).gt.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
    13101212      ENDDO 
    13111213      ! 
     
    13361238                     icorns(ib2,1) = npckgw(ib1) 
    13371239                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN 
    1338                      WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1240                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 
    13391241                        &                                     jpisft(ib2), jpjwft(ib1) 
    1340                      WRITE(ctmp2,*) ' ==========  Not allowed yet' 
    1341                      WRITE(ctmp3,*) '             Crossing problem with West segment: ',npckgw(ib1), &  
    1342                         &                                        ' and South segment: ',npckgs(ib2) 
    1343                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1242                     WRITE(ctmp2,*) ' Not allowed yet' 
     1243                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &  
     1244                        &                            ' and South segment: ',npckgs(ib2) 
     1245                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 
    13441246                  ELSE 
    1345                      WRITE(ctmp1,*) ' E R R O R : Check South and West Open boundary indices' 
    1346                      WRITE(ctmp2,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , & 
    1347                         &                                         ' and South segment: ',npckgs(ib2) 
    1348                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1247                     WRITE(ctmp1,*) ' Check South and West Open boundary indices' 
     1248                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , & 
     1249                        &                            ' and South segment: ',npckgs(ib2) 
     1250                     CALL ctl_stop( ctmp1, ctmp2 ) 
    13491251                  END IF 
    13501252               END IF 
     
    13681270                     icorns(ib2,2) = npckge(ib1) 
    13691271                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN 
    1370                      WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1272                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 
    13711273                        &                                     jpisdt(ib2), jpjeft(ib1) 
    1372                      WRITE(ctmp2,*) ' ==========  Not allowed yet' 
    1373                      WRITE(ctmp3,*) '             Crossing problem with East segment: ',npckge(ib1), & 
    1374                         &                                        ' and South segment: ',npckgs(ib2) 
    1375                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1274                     WRITE(ctmp2,*) ' Not allowed yet' 
     1275                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), & 
     1276                        &                            ' and South segment: ',npckgs(ib2) 
     1277                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 
    13761278                  ELSE 
    1377                      WRITE(ctmp1,*) ' E R R O R : Check South and East Open boundary indices' 
    1378                      WRITE(ctmp2,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
    1379                      &                                           ' and South segment: ',npckgs(ib2) 
    1380                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1279                     WRITE(ctmp1,*) ' Check South and East Open boundary indices' 
     1280                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), & 
     1281                     &                               ' and South segment: ',npckgs(ib2) 
     1282                     CALL ctl_stop( ctmp1, ctmp2 ) 
    13811283                  END IF 
    13821284               END IF 
     
    14001302                     icornn(ib2,1) = npckgw(ib1) 
    14011303                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN 
    1402                      WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1304                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 
    14031305                     &                                     jpinft(ib2), jpjwdt(ib1) 
    1404                      WRITE(ctmp2,*) ' ==========  Not allowed yet' 
    1405                      WRITE(ctmp3,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
    1406                      &                                                    ' and North segment: ',npckgn(ib2) 
    1407                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1306                     WRITE(ctmp2,*) ' Not allowed yet' 
     1307                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
     1308                     &                               ' and North segment: ',npckgn(ib2) 
     1309                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 
    14081310                  ELSE 
    1409                      WRITE(ctmp1,*) ' E R R O R : Check North and West Open boundary indices' 
    1410                      WRITE(ctmp2,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), & 
    1411                      &                                                    ' and North segment: ',npckgn(ib2) 
    1412                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1311                     WRITE(ctmp1,*) ' Check North and West Open boundary indices' 
     1312                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
     1313                     &                               ' and North segment: ',npckgn(ib2) 
     1314                     CALL ctl_stop( ctmp1, ctmp2 ) 
    14131315                  END IF 
    14141316               END IF 
     
    14321334                     icornn(ib2,2) = npckge(ib1) 
    14331335                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN 
    1434                      WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1336                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 
    14351337                     &                                     jpindt(ib2), jpjedt(ib1) 
    1436                      WRITE(ctmp2,*) ' ==========  Not allowed yet' 
    1437                      WRITE(ctmp3,*) '             Crossing problem with East segment: ',npckge(ib1), & 
    1438                      &                                           ' and North segment: ',npckgn(ib2) 
    1439                      CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 
     1338                     WRITE(ctmp2,*) ' Not allowed yet' 
     1339                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), & 
     1340                     &                               ' and North segment: ',npckgn(ib2) 
     1341                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 
    14401342                  ELSE 
    1441                      WRITE(ctmp1,*) ' E R R O R : Check North and East Open boundary indices' 
    1442                      WRITE(ctmp2,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
    1443                      &                                           ' and North segment: ',npckgn(ib2) 
    1444                      CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1343                     WRITE(ctmp1,*) ' Check North and East Open boundary indices' 
     1344                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), & 
     1345                     &                               ' and North segment: ',npckgn(ib2) 
     1346                     CALL ctl_stop( ctmp1, ctmp2 ) 
    14451347                  END IF 
    14461348               END IF 
     
    14681370         IF (ztestmask(1)==1) THEN  
    14691371            IF (icornw(ib,1)==0) THEN 
    1470                WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
    1471                WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                   
    1472                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1372               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib) 
     1373               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 
    14731374            ELSE 
    14741375               ! This is a corner 
     
    14801381         IF (ztestmask(2)==1) THEN 
    14811382            IF (icornw(ib,2)==0) THEN 
    1482                WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
    1483                WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                   
    1484                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1383               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib) 
     1384               CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' ) 
    14851385            ELSE 
    14861386               ! This is a corner 
     
    15081408         IF (ztestmask(1)==1) THEN 
    15091409            IF (icorne(ib,1)==0) THEN 
    1510                WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib) 
    1511                WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                   
    1512                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1410               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib) 
     1411               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 
    15131412            ELSE 
    15141413               ! This is a corner 
     
    15201419         IF (ztestmask(2)==1) THEN 
    15211420            IF (icorne(ib,2)==0) THEN 
    1522                WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib) 
    1523                WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                   
    1524                CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1421               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib) 
     1422               CALL ctl_stop( ctmp1, ' does not end on land or on a corner' ) 
    15251423            ELSE 
    15261424               ! This is a corner 
     
    15471445 
    15481446         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN 
    1549             WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
    1550             WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                   
    1551             CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1447            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib) 
     1448            CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 
    15521449         ENDIF 
    15531450         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN 
    1554             WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
    1555             WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                   
    1556             CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1451            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib) 
     1452            CALL ctl_stop( ctmp1, ' does not end on land or on a corner' ) 
    15571453         ENDIF 
    15581454      END DO 
     
    15731469 
    15741470         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN 
    1575             WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
    1576             WRITE(ctmp2,*) ' ==========  does not start on land'                                                   
    1577             CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1471            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib) 
     1472            CALL ctl_stop( ctmp1, ' does not start on land' ) 
    15781473         ENDIF 
    15791474         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN 
    1580             WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
    1581             WRITE(ctmp2,*) ' ==========  does not end on land'                                                   
    1582             CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1475            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib) 
     1476            CALL ctl_stop( ctmp1, ' does not end on land' ) 
    15831477         ENDIF 
    15841478      END DO 
     
    15931487   END SUBROUTINE bdy_ctl_seg 
    15941488 
    1595  
     1489    
     1490   SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta )  
     1491      !!---------------------------------------------------------------------- 
     1492      !!                 ***  ROUTINE bdy_coords_seg  *** 
     1493      !! 
     1494      !! ** Purpose :  build nbidta, nbidta, nbrdta for bdy built with segments 
     1495      !! 
     1496      !! ** Method  :   
     1497      !! 
     1498      !!---------------------------------------------------------------------- 
     1499      INTEGER, DIMENSION(:,:,:), intent(  out)  ::   nbidta, nbjdta, nbrdta   ! Index arrays: i and j indices of bdy dta 
     1500      !! 
     1501      INTEGER  ::   ii, ij, ir, iseg 
     1502      INTEGER  ::   igrd         ! grid type (t=1, u=2, v=3) 
     1503      INTEGER  ::   icount       !  
     1504      INTEGER  ::   ib_bdy       ! bdy number 
     1505      !!---------------------------------------------------------------------- 
     1506 
     1507      ! East 
     1508      !----- 
     1509      DO iseg = 1, nbdysege 
     1510         ib_bdy = npckge(iseg) 
     1511         ! 
     1512         ! ------------ T points ------------- 
     1513         igrd=1 
     1514         icount=0 
     1515         DO ir = 1, nn_rimwidth(ib_bdy) 
     1516            DO ij = jpjedt(iseg), jpjeft(iseg) 
     1517               icount = icount + 1 
     1518               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
     1519               nbjdta(icount, igrd, ib_bdy) = ij 
     1520               nbrdta(icount, igrd, ib_bdy) = ir 
     1521            ENDDO 
     1522         ENDDO 
     1523         ! 
     1524         ! ------------ U points ------------- 
     1525         igrd=2 
     1526         icount=0 
     1527         DO ir = 1, nn_rimwidth(ib_bdy) 
     1528            DO ij = jpjedt(iseg), jpjeft(iseg) 
     1529               icount = icount + 1 
     1530               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 
     1531               nbjdta(icount, igrd, ib_bdy) = ij 
     1532               nbrdta(icount, igrd, ib_bdy) = ir 
     1533            ENDDO 
     1534         ENDDO 
     1535         ! 
     1536         ! ------------ V points ------------- 
     1537         igrd=3 
     1538         icount=0 
     1539         DO ir = 1, nn_rimwidth(ib_bdy) 
     1540            !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
     1541            DO ij = jpjedt(iseg), jpjeft(iseg) 
     1542               icount = icount + 1 
     1543               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
     1544               nbjdta(icount, igrd, ib_bdy) = ij 
     1545               nbrdta(icount, igrd, ib_bdy) = ir 
     1546            ENDDO 
     1547            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1548            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1549         ENDDO 
     1550      ENDDO 
     1551      ! 
     1552      ! West 
     1553      !----- 
     1554      DO iseg = 1, nbdysegw 
     1555         ib_bdy = npckgw(iseg) 
     1556         ! 
     1557         ! ------------ T points ------------- 
     1558         igrd=1 
     1559         icount=0 
     1560         DO ir = 1, nn_rimwidth(ib_bdy) 
     1561            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     1562               icount = icount + 1 
     1563               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     1564               nbjdta(icount, igrd, ib_bdy) = ij 
     1565               nbrdta(icount, igrd, ib_bdy) = ir 
     1566            ENDDO 
     1567         ENDDO 
     1568         ! 
     1569         ! ------------ U points ------------- 
     1570         igrd=2 
     1571         icount=0 
     1572         DO ir = 1, nn_rimwidth(ib_bdy) 
     1573            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     1574               icount = icount + 1 
     1575               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     1576               nbjdta(icount, igrd, ib_bdy) = ij 
     1577               nbrdta(icount, igrd, ib_bdy) = ir 
     1578            ENDDO 
     1579         ENDDO 
     1580         ! 
     1581         ! ------------ V points ------------- 
     1582         igrd=3 
     1583         icount=0 
     1584         DO ir = 1, nn_rimwidth(ib_bdy) 
     1585            !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
     1586            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     1587               icount = icount + 1 
     1588               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     1589               nbjdta(icount, igrd, ib_bdy) = ij 
     1590               nbrdta(icount, igrd, ib_bdy) = ir 
     1591            ENDDO 
     1592            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1593            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1594         ENDDO 
     1595      ENDDO 
     1596      ! 
     1597      ! North 
     1598      !----- 
     1599      DO iseg = 1, nbdysegn 
     1600         ib_bdy = npckgn(iseg) 
     1601         ! 
     1602         ! ------------ T points ------------- 
     1603         igrd=1 
     1604         icount=0 
     1605         DO ir = 1, nn_rimwidth(ib_bdy) 
     1606            DO ii = jpindt(iseg), jpinft(iseg) 
     1607               icount = icount + 1 
     1608               nbidta(icount, igrd, ib_bdy) = ii 
     1609               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir  
     1610               nbrdta(icount, igrd, ib_bdy) = ir 
     1611            ENDDO 
     1612         ENDDO 
     1613         ! 
     1614         ! ------------ U points ------------- 
     1615         igrd=2 
     1616         icount=0 
     1617         DO ir = 1, nn_rimwidth(ib_bdy) 
     1618            !            DO ii = jpindt(iseg), jpinft(iseg) - 1 
     1619            DO ii = jpindt(iseg), jpinft(iseg) 
     1620               icount = icount + 1 
     1621               nbidta(icount, igrd, ib_bdy) = ii 
     1622               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
     1623               nbrdta(icount, igrd, ib_bdy) = ir 
     1624            ENDDO 
     1625            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1626            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1627         ENDDO 
     1628         ! 
     1629         ! ------------ V points ------------- 
     1630         igrd=3 
     1631         icount=0 
     1632         DO ir = 1, nn_rimwidth(ib_bdy) 
     1633            DO ii = jpindt(iseg), jpinft(iseg) 
     1634               icount = icount + 1 
     1635               nbidta(icount, igrd, ib_bdy) = ii 
     1636               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 
     1637               nbrdta(icount, igrd, ib_bdy) = ir 
     1638            ENDDO 
     1639         ENDDO 
     1640      ENDDO 
     1641      ! 
     1642      ! South 
     1643      !----- 
     1644      DO iseg = 1, nbdysegs 
     1645         ib_bdy = npckgs(iseg) 
     1646         ! 
     1647         ! ------------ T points ------------- 
     1648         igrd=1 
     1649         icount=0 
     1650         DO ir = 1, nn_rimwidth(ib_bdy) 
     1651            DO ii = jpisdt(iseg), jpisft(iseg) 
     1652               icount = icount + 1 
     1653               nbidta(icount, igrd, ib_bdy) = ii 
     1654               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     1655               nbrdta(icount, igrd, ib_bdy) = ir 
     1656            ENDDO 
     1657         ENDDO 
     1658         ! 
     1659         ! ------------ U points ------------- 
     1660         igrd=2 
     1661         icount=0 
     1662         DO ir = 1, nn_rimwidth(ib_bdy) 
     1663            !            DO ii = jpisdt(iseg), jpisft(iseg) - 1 
     1664            DO ii = jpisdt(iseg), jpisft(iseg) 
     1665               icount = icount + 1 
     1666               nbidta(icount, igrd, ib_bdy) = ii 
     1667               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     1668               nbrdta(icount, igrd, ib_bdy) = ir 
     1669            ENDDO 
     1670            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1671            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     1672         ENDDO 
     1673         ! 
     1674         ! ------------ V points ------------- 
     1675         igrd=3 
     1676         icount=0 
     1677         DO ir = 1, nn_rimwidth(ib_bdy) 
     1678            DO ii = jpisdt(iseg), jpisft(iseg) 
     1679               icount = icount + 1 
     1680               nbidta(icount, igrd, ib_bdy) = ii 
     1681               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     1682               nbrdta(icount, igrd, ib_bdy) = ir 
     1683            ENDDO 
     1684         ENDDO 
     1685      ENDDO 
     1686 
     1687       
     1688   END SUBROUTINE bdy_coords_seg 
     1689    
     1690    
    15961691   SUBROUTINE bdy_ctl_corn( ib1, ib2 ) 
    15971692      !!---------------------------------------------------------------------- 
     
    16191714      ! 
    16201715      IF( itest>0 ) THEN 
    1621          WRITE(ctmp1,*) ' E R R O R : Segments ', ib1, 'and ', ib2 
    1622          WRITE(ctmp2,*) ' ==========  have different open bdy schemes'                                                   
    1623          CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 
     1716         WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2 
     1717         CALL ctl_stop( ctmp1, ' have different open bdy schemes' ) 
    16241718      ENDIF 
    16251719      ! 
    16261720   END SUBROUTINE bdy_ctl_corn 
    16271721 
     1722 
     1723   SUBROUTINE bdy_meshwri() 
     1724      !!---------------------------------------------------------------------- 
     1725      !!                 ***  ROUTINE bdy_meshwri  *** 
     1726      !!          
     1727      !! ** Purpose :   write netcdf file with nbr, flagu, flagv, ntreat for T, U  
     1728      !!                and V points in 2D arrays for easier visualisation/control 
     1729      !! 
     1730      !! ** Method  :   use iom_rstput as in domwri.F 
     1731      !!----------------------------------------------------------------------       
     1732      INTEGER  ::   ib_bdy, ii, ij, igrd, ib     ! dummy loop indices 
     1733      INTEGER  ::   inum                                   !   -       - 
     1734      REAL(wp), POINTER, DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields 
     1735      REAL(wp)         , DIMENSION(jpi,jpj) ::   ztmp 
     1736      CHARACTER(LEN=1) , DIMENSION(jpbgrd)  ::   cgrid 
     1737      !!----------------------------------------------------------------------       
     1738      cgrid = (/'t','u','v'/) 
     1739      CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. ) 
     1740      DO igrd = 1, jpbgrd 
     1741         SELECT CASE( igrd ) 
     1742         CASE( 1 )   ;   zmask => tmask(:,:,1) 
     1743         CASE( 2 )   ;   zmask => umask(:,:,1) 
     1744         CASE( 3 )   ;   zmask => vmask(:,:,1) 
     1745         END SELECT 
     1746         ztmp(:,:) = zmask(:,:) 
     1747         DO ib_bdy = 1, nb_bdy 
     1748            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)      ! nbr deined for all rims 
     1749               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1750               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1751               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10. 
     1752               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 
     1753            END DO 
     1754         END DO 
     1755         CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 
     1756         ztmp(:,:) = zmask(:,:) 
     1757         DO ib_bdy = 1, nb_bdy 
     1758            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagu defined only for rims 0 and 1 
     1759               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1760               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1761               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10. 
     1762               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 
     1763            END DO 
     1764         END DO 
     1765         CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 
     1766         ztmp(:,:) = zmask(:,:) 
     1767         DO ib_bdy = 1, nb_bdy 
     1768            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagv defined only for rims 0 and 1 
     1769               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1770               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1771               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10. 
     1772               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 
     1773            END DO 
     1774         END DO 
     1775         CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 
     1776         ztmp(:,:) = zmask(:,:) 
     1777         DO ib_bdy = 1, nb_bdy 
     1778            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! ntreat defined only for rims 0 and 1 
     1779               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     1780               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     1781               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10. 
     1782               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 
     1783            END DO 
     1784         END DO 
     1785         CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 
     1786      END DO 
     1787      CALL iom_close( inum ) 
     1788 
     1789   END SUBROUTINE bdy_meshwri 
     1790    
    16281791   !!================================================================================= 
    16291792END MODULE bdyini 
Note: See TracChangeset for help on using the changeset viewer.