New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3583 for branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 – NEMO

Ignore:
Timestamp:
2012-11-16T17:18:17+01:00 (11 years ago)
Author:
cbricaud
Message:

add modification from dev_r3327_MERCATOR1_BDY branch in dev_MERCATOR_2012_rev3555 branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r3424 r3583  
    1111   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions 
    1212   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     13   !!            3.4  !  2012     (J. Chanut) straight open boundary case update 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_bdy 
     
    2627   USE lib_mpp         ! for mpp_sum   
    2728   USE iom             ! I/O 
     29   USE sbctide, ONLY: lk_tide ! Tidal forcing or not 
     30   USE phycst, ONLY: rday 
    2831 
    2932   IMPLICIT NONE 
     
    3235   PUBLIC   bdy_init   ! routine called in nemo_init 
    3336 
     37   INTEGER, PARAMETER          :: jp_nseg = 100 
     38   INTEGER, PARAMETER          :: nrimmax = 20 ! maximum rimwidth in structured 
     39                                               ! open boundary data files 
     40   ! Straight open boundary segment parameters: 
     41   INTEGER  :: nbdysege, nbdysegw, nbdysegn, nbdysegs  
     42   INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge 
     43   INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw 
     44   INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn 
     45   INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs 
    3446   !!---------------------------------------------------------------------- 
    3547   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     
    5365      ! namelist variables 
    5466      !------------------- 
    55       INTEGER, PARAMETER          :: jp_nseg = 100 
    56       INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs  
    57       INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft 
    58       INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft 
    59       INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft 
    60       INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft 
     67      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile 
     68      CHARACTER(LEN=1)   ::   ctypebdy 
     69      INTEGER :: nbdyind, nbdybeg, nbdyend 
    6170 
    6271      ! local variables 
     
    6675      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       - 
    6776      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       - 
     77      INTEGER  ::   jpbdtau, jpbdtas                       !   -       - 
     78      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       - 
    6879      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts 
    6980      REAL   , POINTER  ::  flagu, flagv                   !    -   - 
    7081      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    71       INTEGER, DIMENSION (2)                ::   kdimsz 
     82      INTEGER, DIMENSION (2)                  ::   kdimsz 
    7283      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays  
    7384      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
    7485      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points 
    75       REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask 
    76       CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile 
    77       CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid 
     86      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid 
    7887      !! 
    7988      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             & 
    8089         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, & 
    8190         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &   
     91         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,              & 
    8292#if defined key_lim2 
    8393         &             nn_ice_lim2, nn_ice_lim2_dta,                       & 
     
    8595         &             ln_vol, nn_volctl, nn_rimwidth 
    8696      !! 
    87       NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             & 
    88                              nbdysegw, jpiwob, jpjwdt, jpjwft,             & 
    89                              nbdysegn, jpjnob, jpindt, jpinft,             & 
    90                              nbdysegs, jpjsob, jpisdt, jpisft 
     97      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
    9198 
    9299      !!---------------------------------------------------------------------- 
     
    105112 
    106113      cgrid= (/'t','u','v'/) 
    107  
     114       
    108115      ! ----------------------------------------- 
    109116      ! Initialise and read namelist parameters 
     
    121128      nn_tra(:)         = 0 
    122129      nn_tra_dta(:)     = -1  ! uninitialised flag 
     130      ln_tra_dmp(:)     = .false. 
     131      ln_dyn3d_dmp(:)   = .false. 
     132      rn_time_dmp(:)    = 1. 
    123133#if defined key_lim2 
    124134      nn_ice_lim2(:)    = 0 
     
    135145      ! Check and write out namelist parameters 
    136146      ! ----------------------------------------- 
    137  
    138147      !                                   ! control prints 
    139148      IF(lwp) WRITE(numout,*) '         nambdy' 
     
    158167        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
    159168        SELECT CASE( nn_dyn2d(ib_bdy) )                   
    160           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    161           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    162           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
     169          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     170          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     171          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    163172          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
    164173        END SELECT 
     
    171180              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
    172181           END SELECT 
     182           IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.lk_tide)) THEN 
     183             CALL ctl_stop( 'You must activate key_tide to add tidal forcing at open boundaries' ) 
     184           ENDIF 
    173185        ENDIF 
    174186        IF(lwp) WRITE(numout,*) 
     
    176188        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
    177189        SELECT CASE( nn_dyn3d(ib_bdy) )                   
    178           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    179           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     190          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     191          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     192          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
     193          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
    180194          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 
    181195        END SELECT 
     
    187201           END SELECT 
    188202        ENDIF 
     203 
     204        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
     205           IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN 
     206              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
     207              ln_dyn3d_dmp(ib_bdy)=.false. 
     208           ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN 
     209              CALL ctl_stop( 'Use FRS OR relaxation' ) 
     210           ELSE 
     211              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
     212              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     213              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     214           ENDIF 
     215        ELSE 
     216           IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities' 
     217        ENDIF 
    189218        IF(lwp) WRITE(numout,*) 
    190219 
    191220        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
    192221        SELECT CASE( nn_tra(ib_bdy) )                   
    193           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    194           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     222          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     223          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     224          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
     225          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions' 
     226          CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
    195227          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
    196228        END SELECT 
     
    201233              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
    202234           END SELECT 
     235        ENDIF 
     236 
     237        IF ( ln_tra_dmp(ib_bdy) ) THEN 
     238           IF ( nn_tra(ib_bdy).EQ.0 ) THEN 
     239              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
     240              ln_tra_dmp(ib_bdy)=.false. 
     241           ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN 
     242              CALL ctl_stop( 'Use FRS OR relaxation' ) 
     243           ELSE 
     244              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
     245              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     246              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
     247           ENDIF 
     248        ELSE 
     249           IF(lwp) WRITE(numout,*) '      NO T/S relaxation' 
    203250        ENDIF 
    204251        IF(lwp) WRITE(numout,*) 
     
    221268#endif 
    222269 
    223         IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy) 
     270        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
    224271        IF(lwp) WRITE(numout,*) 
    225272 
    226273      ENDDO 
    227274 
    228      IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
    229        IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
    230        IF(lwp) WRITE(numout,*) 
    231        SELECT CASE ( nn_volctl ) 
    232          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    233          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    234          CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
    235        END SELECT 
    236        IF(lwp) WRITE(numout,*) 
    237      ELSE 
    238        IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
    239        IF(lwp) WRITE(numout,*) 
     275     IF (nb_bdy .gt. 0) THEN 
     276        IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
     277          IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
     278          IF(lwp) WRITE(numout,*) 
     279          SELECT CASE ( nn_volctl ) 
     280            CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
     281            CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
     282            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
     283          END SELECT 
     284          IF(lwp) WRITE(numout,*) 
     285        ELSE 
     286          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
     287          IF(lwp) WRITE(numout,*) 
     288        ENDIF 
    240289     ENDIF 
    241290 
     
    247296      ! --------------------------------------------- 
    248297      REWIND( numnam )                     
     298                
     299      nblendta(:,:) = 0 
     300      nbdysege = 0 
     301      nbdysegw = 0 
     302      nbdysegn = 0 
     303      nbdysegs = 0 
     304      icount   = 0 ! count user defined segments 
     305      ! Dimensions below are used to allocate arrays to read external data 
     306      jpbdtas = 1 ! Maximum size of boundary data (structured case) 
     307      jpbdtau = 1 ! Maximum size of boundary data (unstructured case) 
     308 
    249309      DO ib_bdy = 1, nb_bdy 
    250310 
    251          jpbdta = 1 
    252311         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
    253312  
     313            icount = icount + 1 
    254314            ! No REWIND here because may need to read more than one nambdy_index namelist. 
    255315            READ  ( numnam, nambdy_index ) 
    256316 
    257             ! Automatic boundary definition: if nbdysegX = -1 
    258             ! set boundary to whole side of model domain. 
    259             IF( nbdysege == -1 ) THEN 
    260                nbdysege = 1 
    261                jpieob(1) = jpiglo - 1 
    262                jpjedt(1) = 2 
    263                jpjeft(1) = jpjglo - 1 
    264             ENDIF 
    265             IF( nbdysegw == -1 ) THEN 
    266                nbdysegw = 1 
    267                jpiwob(1) = 2 
    268                jpjwdt(1) = 2 
    269                jpjwft(1) = jpjglo - 1 
    270             ENDIF 
    271             IF( nbdysegn == -1 ) THEN 
    272                nbdysegn = 1 
    273                jpjnob(1) = jpjglo - 1 
    274                jpindt(1) = 2 
    275                jpinft(1) = jpiglo - 1 
    276             ENDIF 
    277             IF( nbdysegs == -1 ) THEN 
    278                nbdysegs = 1 
    279                jpjsob(1) = 2 
    280                jpisdt(1) = 2 
    281                jpisft(1) = jpiglo - 1 
    282             ENDIF 
    283  
    284             nblendta(:,ib_bdy) = 0 
    285             DO iseg = 1, nbdysege 
    286                igrd = 1 
    287                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
    288                igrd = 2 
    289                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
    290                igrd = 3 
    291                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)                
    292             ENDDO 
    293             DO iseg = 1, nbdysegw 
    294                igrd = 1 
    295                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
    296                igrd = 2 
    297                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
    298                igrd = 3 
    299                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)                
    300             ENDDO 
    301             DO iseg = 1, nbdysegn 
    302                igrd = 1 
    303                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1                
    304                igrd = 2 
    305                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)                
    306                igrd = 3 
    307                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1 
    308             ENDDO 
    309             DO iseg = 1, nbdysegs 
    310                igrd = 1 
    311                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
    312                igrd = 2 
    313                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) 
    314                igrd = 3 
    315                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
    316             ENDDO 
    317  
    318             nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 
    319             jpbdta = MAXVAL(nblendta(:,ib_bdy))                
    320  
     317            SELECT CASE ( TRIM(ctypebdy) ) 
     318              CASE( 'N' ) 
     319                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     320                    nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain. 
     321                    nbdybeg  = 2 
     322                    nbdyend  = jpiglo - 1 
     323                 ENDIF 
     324                 nbdysegn = nbdysegn + 1 
     325                 npckgn(nbdysegn) = ib_bdy ! Save bdy package number 
     326                 jpjnob(nbdysegn) = nbdyind 
     327                 jpindt(nbdysegn) = nbdybeg 
     328                 jpinft(nbdysegn) = nbdyend 
     329                 ! 
     330              CASE( 'S' ) 
     331                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     332                    nbdyind  = 2           ! set boundary to whole side of model domain. 
     333                    nbdybeg  = 2 
     334                    nbdyend  = jpiglo - 1 
     335                 ENDIF 
     336                 nbdysegs = nbdysegs + 1 
     337                 npckgs(nbdysegs) = ib_bdy ! Save bdy package number 
     338                 jpjsob(nbdysegs) = nbdyind 
     339                 jpisdt(nbdysegs) = nbdybeg 
     340                 jpisft(nbdysegs) = nbdyend 
     341                 ! 
     342              CASE( 'E' ) 
     343                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     344                    nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain. 
     345                    nbdybeg  = 2 
     346                    nbdyend  = jpjglo - 1 
     347                 ENDIF 
     348                 nbdysege = nbdysege + 1  
     349                 npckge(nbdysege) = ib_bdy ! Save bdy package number 
     350                 jpieob(nbdysege) = nbdyind 
     351                 jpjedt(nbdysege) = nbdybeg 
     352                 jpjeft(nbdysege) = nbdyend 
     353                 ! 
     354              CASE( 'W' ) 
     355                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
     356                    nbdyind  = 2           ! set boundary to whole side of model domain. 
     357                    nbdybeg  = 2 
     358                    nbdyend  = jpjglo - 1 
     359                 ENDIF 
     360                 nbdysegw = nbdysegw + 1 
     361                 npckgw(nbdysegw) = ib_bdy ! Save bdy package number 
     362                 jpiwob(nbdysegw) = nbdyind 
     363                 jpjwdt(nbdysegw) = nbdybeg 
     364                 jpjwft(nbdysegw) = nbdyend 
     365                 ! 
     366              CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 
     367            END SELECT 
     368 
     369            ! For simplicity we assume that in case of straight bdy, arrays have the same length 
     370            ! (even if it is true that last tangential velocity points 
     371            ! are useless). This simplifies a little bit boundary data format (and agrees with format 
     372            ! used so far in obc package) 
     373 
     374            nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy) 
     375            jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1)) 
     376            IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) & 
     377            & CALL ctl_stop( 'rimwidth must be lower than nrimmax' ) 
    321378 
    322379         ELSE            ! Read size of arrays in boundary coordinates file. 
    323  
    324  
    325380            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    326             jpbdta = 1 
    327381            DO igrd = 1, jpbgrd 
    328382               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    329383               nblendta(igrd,ib_bdy) = kdimsz(1) 
    330                jpbdta = MAX(jpbdta, kdimsz(1)) 
    331             ENDDO 
     384               jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     385            ENDDO 
     386            CALL iom_close( inum ) 
    332387 
    333388         ENDIF  
     
    335390      ENDDO ! ib_bdy 
    336391 
    337       ! Allocate arrays 
    338       !--------------- 
    339       ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
    340          &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
    341  
    342       ALLOCATE( dta_global(jpbdta, 1, jpk) ) 
     392      IF (nb_bdy>0) THEN 
     393         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
     394 
     395         ! Allocate arrays 
     396         !--------------- 
     397         ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
     398            &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
     399 
     400         ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 
     401         IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 
     402         !  
     403      ENDIF 
     404 
     405      ! Now look for crossings in user (namelist) defined open boundary segments: 
     406      !-------------------------------------------------------------------------- 
     407      IF ( icount>0 ) CALL bdy_ctl_seg 
    343408 
    344409      ! Calculate global boundary index arrays or read in from file 
    345       !------------------------------------------------------------ 
    346       REWIND( numnam )                     
     410      !------------------------------------------------------------                
     411      ! 1. Read global index arrays from boundary coordinates file. 
    347412      DO ib_bdy = 1, nb_bdy 
    348413 
    349          IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters 
    350  
    351             ! No REWIND here because may need to read more than one nambdy_index namelist. 
    352             READ  ( numnam, nambdy_index ) 
    353  
    354             ! Automatic boundary definition: if nbdysegX = -1 
    355             ! set boundary to whole side of model domain. 
    356             IF( nbdysege == -1 ) THEN 
    357                nbdysege = 1 
    358                jpieob(1) = jpiglo - 1 
    359                jpjedt(1) = 2 
    360                jpjeft(1) = jpjglo - 1 
    361             ENDIF 
    362             IF( nbdysegw == -1 ) THEN 
    363                nbdysegw = 1 
    364                jpiwob(1) = 2 
    365                jpjwdt(1) = 2 
    366                jpjwft(1) = jpjglo - 1 
    367             ENDIF 
    368             IF( nbdysegn == -1 ) THEN 
    369                nbdysegn = 1 
    370                jpjnob(1) = jpjglo - 1 
    371                jpindt(1) = 2 
    372                jpinft(1) = jpiglo - 1 
    373             ENDIF 
    374             IF( nbdysegs == -1 ) THEN 
    375                nbdysegs = 1 
    376                jpjsob(1) = 2 
    377                jpisdt(1) = 2 
    378                jpisft(1) = jpiglo - 1 
    379             ENDIF 
    380  
    381             ! ------------ T points ------------- 
    382             igrd = 1   
    383             icount = 0 
    384             DO ir = 1, nn_rimwidth(ib_bdy) 
    385                ! east 
    386                DO iseg = 1, nbdysege 
    387                   DO ij = jpjedt(iseg), jpjeft(iseg) 
    388                      icount = icount + 1 
    389                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
    390                      nbjdta(icount, igrd, ib_bdy) = ij 
    391                      nbrdta(icount, igrd, ib_bdy) = ir 
    392                   ENDDO 
    393                ENDDO 
    394                ! west 
    395                DO iseg = 1, nbdysegw 
    396                   DO ij = jpjwdt(iseg), jpjwft(iseg) 
    397                      icount = icount + 1 
    398                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    399                      nbjdta(icount, igrd, ib_bdy) = ij 
    400                      nbrdta(icount, igrd, ib_bdy) = ir 
    401                   ENDDO 
    402                ENDDO 
    403                ! north 
    404                DO iseg = 1, nbdysegn 
    405                   DO ii = jpindt(iseg), jpinft(iseg) 
    406                      icount = icount + 1 
    407                      nbidta(icount, igrd, ib_bdy) = ii 
    408                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
    409                      nbrdta(icount, igrd, ib_bdy) = ir 
    410                   ENDDO 
    411                ENDDO 
    412                ! south 
    413                DO iseg = 1, nbdysegs 
    414                   DO ii = jpisdt(iseg), jpisft(iseg) 
    415                      icount = icount + 1 
    416                      nbidta(icount, igrd, ib_bdy) = ii 
    417                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    418                      nbrdta(icount, igrd, ib_bdy) = ir 
    419                   ENDDO 
    420                ENDDO 
    421             ENDDO 
    422  
    423             ! ------------ U points ------------- 
    424             igrd = 2   
    425             icount = 0 
    426             DO ir = 1, nn_rimwidth(ib_bdy) 
    427                ! east 
    428                DO iseg = 1, nbdysege 
    429                   DO ij = jpjedt(iseg), jpjeft(iseg) 
    430                      icount = icount + 1 
    431                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir 
    432                      nbjdta(icount, igrd, ib_bdy) = ij 
    433                      nbrdta(icount, igrd, ib_bdy) = ir 
    434                   ENDDO 
    435                ENDDO 
    436                ! west 
    437                DO iseg = 1, nbdysegw 
    438                   DO ij = jpjwdt(iseg), jpjwft(iseg) 
    439                      icount = icount + 1 
    440                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    441                      nbjdta(icount, igrd, ib_bdy) = ij 
    442                      nbrdta(icount, igrd, ib_bdy) = ir 
    443                   ENDDO 
    444                ENDDO 
    445                ! north 
    446                DO iseg = 1, nbdysegn 
    447                   DO ii = jpindt(iseg), jpinft(iseg) - 1 
    448                      icount = icount + 1 
    449                      nbidta(icount, igrd, ib_bdy) = ii 
    450                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
    451                      nbrdta(icount, igrd, ib_bdy) = ir 
    452                   ENDDO 
    453                ENDDO 
    454                ! south 
    455                DO iseg = 1, nbdysegs 
    456                   DO ii = jpisdt(iseg), jpisft(iseg) - 1 
    457                      icount = icount + 1 
    458                      nbidta(icount, igrd, ib_bdy) = ii 
    459                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    460                      nbrdta(icount, igrd, ib_bdy) = ir 
    461                   ENDDO 
    462                ENDDO 
    463             ENDDO 
    464  
    465             ! ------------ V points ------------- 
    466             igrd = 3   
    467             icount = 0 
    468             DO ir = 1, nn_rimwidth(ib_bdy) 
    469                ! east 
    470                DO iseg = 1, nbdysege 
    471                   DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
    472                      icount = icount + 1 
    473                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
    474                      nbjdta(icount, igrd, ib_bdy) = ij 
    475                      nbrdta(icount, igrd, ib_bdy) = ir 
    476                   ENDDO 
    477                ENDDO 
    478                ! west 
    479                DO iseg = 1, nbdysegw 
    480                   DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
    481                      icount = icount + 1 
    482                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    483                      nbjdta(icount, igrd, ib_bdy) = ij 
    484                      nbrdta(icount, igrd, ib_bdy) = ir 
    485                   ENDDO 
    486                ENDDO 
    487                ! north 
    488                DO iseg = 1, nbdysegn 
    489                   DO ii = jpindt(iseg), jpinft(iseg) 
    490                      icount = icount + 1 
    491                      nbidta(icount, igrd, ib_bdy) = ii 
    492                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir 
    493                      nbrdta(icount, igrd, ib_bdy) = ir 
    494                   ENDDO 
    495                ENDDO 
    496                ! south 
    497                DO iseg = 1, nbdysegs 
    498                   DO ii = jpisdt(iseg), jpisft(iseg) 
    499                      icount = icount + 1 
    500                      nbidta(icount, igrd, ib_bdy) = ii 
    501                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    502                      nbrdta(icount, igrd, ib_bdy) = ir 
    503                   ENDDO 
    504                ENDDO 
    505             ENDDO 
    506  
    507          ELSE            ! Read global index arrays from boundary coordinates file. 
    508  
     414         IF( ln_coords_file(ib_bdy) ) THEN 
     415 
     416            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    509417            DO igrd = 1, jpbgrd 
    510418               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     
    527435               IF (ibr_max < nn_rimwidth(ib_bdy))   & 
    528436                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
    529  
    530437            END DO 
    531438            CALL iom_close( inum ) 
     
    533440         ENDIF  
    534441 
    535       ENDDO  
     442      ENDDO       
     443     
     444      ! 2. Now fill indices corresponding to straight open boundary arrays: 
     445      ! East 
     446      !----- 
     447      DO iseg = 1, nbdysege 
     448         ib_bdy = npckge(iseg) 
     449         ! 
     450         ! ------------ T points ------------- 
     451         igrd=1 
     452         icount=0 
     453         DO ir = 1, nn_rimwidth(ib_bdy) 
     454            DO ij = jpjedt(iseg), jpjeft(iseg) 
     455               icount = icount + 1 
     456               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
     457               nbjdta(icount, igrd, ib_bdy) = ij 
     458               nbrdta(icount, igrd, ib_bdy) = ir 
     459            ENDDO 
     460         ENDDO 
     461         ! 
     462         ! ------------ U points ------------- 
     463         igrd=2 
     464         icount=0 
     465         DO ir = 1, nn_rimwidth(ib_bdy) 
     466            DO ij = jpjedt(iseg), jpjeft(iseg) 
     467               icount = icount + 1 
     468               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 
     469               nbjdta(icount, igrd, ib_bdy) = ij 
     470               nbrdta(icount, igrd, ib_bdy) = ir 
     471            ENDDO 
     472         ENDDO 
     473         ! 
     474         ! ------------ V points ------------- 
     475         igrd=3 
     476         icount=0 
     477         DO ir = 1, nn_rimwidth(ib_bdy) 
     478!            DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
     479            DO ij = jpjedt(iseg), jpjeft(iseg) 
     480               icount = icount + 1 
     481               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
     482               nbjdta(icount, igrd, ib_bdy) = ij 
     483               nbrdta(icount, igrd, ib_bdy) = ir 
     484            ENDDO 
     485            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     486            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     487         ENDDO 
     488      ENDDO 
     489      ! 
     490      ! West 
     491      !----- 
     492      DO iseg = 1, nbdysegw 
     493         ib_bdy = npckgw(iseg) 
     494         ! 
     495         ! ------------ T points ------------- 
     496         igrd=1 
     497         icount=0 
     498         DO ir = 1, nn_rimwidth(ib_bdy) 
     499            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     500               icount = icount + 1 
     501               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     502               nbjdta(icount, igrd, ib_bdy) = ij 
     503               nbrdta(icount, igrd, ib_bdy) = ir 
     504            ENDDO 
     505         ENDDO 
     506         ! 
     507         ! ------------ U points ------------- 
     508         igrd=2 
     509         icount=0 
     510         DO ir = 1, nn_rimwidth(ib_bdy) 
     511            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     512               icount = icount + 1 
     513               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     514               nbjdta(icount, igrd, ib_bdy) = ij 
     515               nbrdta(icount, igrd, ib_bdy) = ir 
     516            ENDDO 
     517         ENDDO 
     518         ! 
     519         ! ------------ V points ------------- 
     520         igrd=3 
     521         icount=0 
     522         DO ir = 1, nn_rimwidth(ib_bdy) 
     523!            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
     524            DO ij = jpjwdt(iseg), jpjwft(iseg) 
     525               icount = icount + 1 
     526               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     527               nbjdta(icount, igrd, ib_bdy) = ij 
     528               nbrdta(icount, igrd, ib_bdy) = ir 
     529            ENDDO 
     530            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     531            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     532         ENDDO 
     533      ENDDO 
     534      ! 
     535      ! North 
     536      !----- 
     537      DO iseg = 1, nbdysegn 
     538         ib_bdy = npckgn(iseg) 
     539         ! 
     540         ! ------------ T points ------------- 
     541         igrd=1 
     542         icount=0 
     543         DO ir = 1, nn_rimwidth(ib_bdy) 
     544            DO ii = jpindt(iseg), jpinft(iseg) 
     545               icount = icount + 1 
     546               nbidta(icount, igrd, ib_bdy) = ii 
     547               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir  
     548               nbrdta(icount, igrd, ib_bdy) = ir 
     549            ENDDO 
     550         ENDDO 
     551         ! 
     552         ! ------------ U points ------------- 
     553         igrd=2 
     554         icount=0 
     555         DO ir = 1, nn_rimwidth(ib_bdy) 
     556!            DO ii = jpindt(iseg), jpinft(iseg) - 1 
     557            DO ii = jpindt(iseg), jpinft(iseg) 
     558               icount = icount + 1 
     559               nbidta(icount, igrd, ib_bdy) = ii 
     560               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
     561               nbrdta(icount, igrd, ib_bdy) = ir 
     562            ENDDO 
     563            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     564            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     565         ENDDO 
     566         ! 
     567         ! ------------ V points ------------- 
     568         igrd=3 
     569         icount=0 
     570         DO ir = 1, nn_rimwidth(ib_bdy) 
     571            DO ii = jpindt(iseg), jpinft(iseg) 
     572               icount = icount + 1 
     573               nbidta(icount, igrd, ib_bdy) = ii 
     574               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 
     575               nbrdta(icount, igrd, ib_bdy) = ir 
     576            ENDDO 
     577         ENDDO 
     578      ENDDO 
     579      ! 
     580      ! South 
     581      !----- 
     582      DO iseg = 1, nbdysegs 
     583         ib_bdy = npckgs(iseg) 
     584         ! 
     585         ! ------------ T points ------------- 
     586         igrd=1 
     587         icount=0 
     588         DO ir = 1, nn_rimwidth(ib_bdy) 
     589            DO ii = jpisdt(iseg), jpisft(iseg) 
     590               icount = icount + 1 
     591               nbidta(icount, igrd, ib_bdy) = ii 
     592               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     593               nbrdta(icount, igrd, ib_bdy) = ir 
     594            ENDDO 
     595         ENDDO 
     596         ! 
     597         ! ------------ U points ------------- 
     598         igrd=2 
     599         icount=0 
     600         DO ir = 1, nn_rimwidth(ib_bdy) 
     601!            DO ii = jpisdt(iseg), jpisft(iseg) - 1 
     602            DO ii = jpisdt(iseg), jpisft(iseg) 
     603               icount = icount + 1 
     604               nbidta(icount, igrd, ib_bdy) = ii 
     605               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     606               nbrdta(icount, igrd, ib_bdy) = ir 
     607            ENDDO 
     608            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     609            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
     610         ENDDO 
     611         ! 
     612         ! ------------ V points ------------- 
     613         igrd=3 
     614         icount=0 
     615         DO ir = 1, nn_rimwidth(ib_bdy) 
     616            DO ii = jpisdt(iseg), jpisft(iseg) 
     617               icount = icount + 1 
     618               nbidta(icount, igrd, ib_bdy) = ii 
     619               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     620               nbrdta(icount, igrd, ib_bdy) = ir 
     621            ENDDO 
     622         ENDDO 
     623      ENDDO 
     624 
     625      !  Deal with duplicated points 
     626      !----------------------------- 
     627      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated) 
     628      ! if their distance to the bdy is greater than the other 
     629      ! If their distance are the same, just keep only one to avoid updating a point twice 
     630      DO igrd = 1, jpbgrd 
     631         DO ib_bdy1 = 1, nb_bdy 
     632            DO ib_bdy2 = 1, nb_bdy 
     633               IF (ib_bdy1/=ib_bdy2) THEN 
     634                  DO ib1 = 1, nblendta(igrd,ib_bdy1) 
     635                     DO ib2 = 1, nblendta(igrd,ib_bdy2) 
     636                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. & 
     637                        &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 
     638!                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &  
     639!                                                       &              nbidta(ib1, igrd, ib_bdy1),      &  
     640!                                                       &              nbjdta(ib2, igrd, ib_bdy2) 
     641                           ! keep only points with the lowest distance to boundary: 
     642                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN 
     643                             nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 
     644                             nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 
     645                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN 
     646                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     647                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     648                           ! Arbitrary choice if distances are the same: 
     649                           ELSE 
     650                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     651                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
     652                           ENDIF 
     653                        END IF 
     654                     END DO 
     655                  END DO 
     656               ENDIF 
     657            END DO 
     658         END DO 
     659      END DO 
    536660 
    537661      ! Work out dimensions of boundary data on each processor 
    538662      ! ------------------------------------------------------ 
    539       
    540       iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
    541       ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
    542       is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
    543       in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1 
     663 
     664      ! Rather assume that boundary data indices are given on global domain 
     665      ! TO BE DISCUSSED ? 
     666!      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
     667!      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
     668!      is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
     669!      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1       
     670      iw = mig(1) - jpizoom + 2         ! if monotasking and no zoom, iw=2 
     671      ie = mig(1) + nlci - jpizoom - 1  ! if monotasking and no zoom, ie=jpim1 
     672      is = mjg(1) - jpjzoom + 2         ! if monotasking and no zoom, is=2 
     673      in = mjg(1) + nlcj - jpjzoom - 1  ! if monotasking and no zoom, in=jpjm1 
    544674 
    545675      DO ib_bdy = 1, nb_bdy 
     
    577707         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
    578708         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 
     709         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) ) 
    579710         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 
    580711         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
     
    596727                     ! 
    597728                     icount = icount  + 1 
    598                      idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
    599                      idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
     729 
     730                     ! Rather assume that boundary data indices are given on global domain 
     731                     ! TO BE DISCUSSED ? 
     732!                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
     733!                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
     734                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom 
     735                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom 
    600736                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy) 
    601737                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 
     
    611747               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    612748               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
    613 !              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic 
    614 !              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear 
     749!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
     750!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear 
     751            END DO 
     752         END DO  
     753 
     754         ! Compute damping coefficients 
     755         ! ---------------------------- 
     756         DO igrd = 1, jpbgrd 
     757            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     758               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
     759               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &  
     760               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
    615761            END DO 
    616762         END DO  
     
    627773      !          = 0  elsewhere    
    628774  
    629       IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution 
    630          zmask(         :                ,:) = 0.e0 
    631          zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0           
    632       ELSE IF( ln_mask_file ) THEN 
     775      IF( ln_mask_file ) THEN 
    633776         CALL iom_open( cn_mask_file, inum ) 
    634          CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 
     777         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) ) 
    635778         CALL iom_close( inum ) 
    636       ELSE 
    637          zmask(:,:) = 1.e0 
    638       ENDIF 
    639  
    640       DO ij = 1, nlcj      ! Save mask over local domain       
    641          DO ii = 1, nlci 
    642             bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) ) 
     779 
     780         ! Derive mask on U and V grid from mask on T grid 
     781         bdyumask(:,:) = 0.e0 
     782         bdyvmask(:,:) = 0.e0 
     783         DO ij=1, jpjm1 
     784            DO ii=1, jpim1 
     785               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 
     786               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)   
     787            END DO 
    643788         END DO 
    644       END DO 
    645  
    646       ! Derive mask on U and V grid from mask on T grid 
    647       bdyumask(:,:) = 0.e0 
    648       bdyvmask(:,:) = 0.e0 
    649       DO ij=1, jpjm1 
    650          DO ii=1, jpim1 
    651             bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 
    652             bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)   
     789         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond. 
     790 
     791 
     792         ! Mask corrections 
     793         ! ---------------- 
     794         DO ik = 1, jpkm1 
     795            DO ij = 1, jpj 
     796               DO ii = 1, jpi 
     797                  tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 
     798                  umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 
     799                  vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 
     800                  bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij) 
     801               END DO       
     802            END DO 
    653803         END DO 
    654       END DO 
    655       CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond. 
    656  
    657  
    658       ! Mask corrections 
    659       ! ---------------- 
    660       DO ik = 1, jpkm1 
    661          DO ij = 1, jpj 
    662             DO ii = 1, jpi 
    663                tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 
    664                umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 
    665                vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 
    666                bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij) 
    667             END DO       
     804 
     805         DO ik = 1, jpkm1 
     806            DO ij = 2, jpjm1 
     807               DO ii = 2, jpim1 
     808                  fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   & 
     809                     &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 
     810               END DO       
     811            END DO 
    668812         END DO 
    669       END DO 
    670  
    671       DO ik = 1, jpkm1 
    672          DO ij = 2, jpjm1 
    673             DO ii = 2, jpim1 
    674                fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   & 
    675                   &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 
    676             END DO       
    677          END DO 
    678       END DO 
    679  
    680       tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)              
     813 
     814         tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 
     815 
     816      ENDIF ! ln_mask_file=.TRUE. 
     817       
    681818      bdytmask(:,:) = tmask(:,:,1) 
    682819 
     
    753890            END IF 
    754891         END DO 
    755   
     892 
    756893         IF( icount /= 0 ) THEN 
    757894            IF(lwp) WRITE(numout,*) 
     
    767904      ! Compute total lateral surface for volume correction: 
    768905      ! ---------------------------------------------------- 
     906      ! JC: this must be done at each time step with key_vvl 
    769907      bdysurftot = 0.e0  
    770908      IF( ln_vol ) THEN   
     
    800938      ! Tidy up 
    801939      !-------- 
    802       DEALLOCATE(nbidta, nbjdta, nbrdta) 
     940      IF (nb_bdy>0) THEN 
     941         DEALLOCATE(nbidta, nbjdta, nbrdta) 
     942      ENDIF 
    803943 
    804944      IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 
    805945 
    806946   END SUBROUTINE bdy_init 
     947 
     948   SUBROUTINE bdy_ctl_seg 
     949      !!---------------------------------------------------------------------- 
     950      !!                 ***  ROUTINE bdy_ctl_seg  *** 
     951      !! 
     952      !! ** Purpose :   Check straight open boundary segments location 
     953      !! 
     954      !! ** Method  :   - Look for open boundary corners 
     955      !!                - Check that segments start or end on land  
     956      !!---------------------------------------------------------------------- 
     957      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest   
     958      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns   
     959      REAL(wp), DIMENSION(2) ::   ztestmask 
     960      !!---------------------------------------------------------------------- 
     961      ! 
     962      IF (lwp) WRITE(numout,*) ' ' 
     963      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments' 
     964      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     965      ! 
     966      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege 
     967      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw 
     968      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn 
     969      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs 
     970      ! 1. Check bounds 
     971      !---------------- 
     972      DO ib = 1, nbdysegn 
     973         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib) 
     974         IF ((jpjnob(ib).ge.jpjglo-1).or.&  
     975            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
     976         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
     977         IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     978         IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
     979      END DO 
     980      ! 
     981      DO ib = 1, nbdysegs 
     982         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib) 
     983         IF ((jpjsob(ib).ge.jpjglo-1).or.&  
     984            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
     985         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
     986         IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     987         IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
     988      END DO 
     989      ! 
     990      DO ib = 1, nbdysege 
     991         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib) 
     992         IF ((jpieob(ib).ge.jpiglo-1).or.&  
     993            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
     994         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
     995         IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     996         IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
     997      END DO 
     998      ! 
     999      DO ib = 1, nbdysegw 
     1000         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib) 
     1001         IF ((jpiwob(ib).ge.jpiglo-1).or.&  
     1002            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
     1003         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
     1004         IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
     1005         IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
     1006      ENDDO 
     1007      ! 
     1008      !       
     1009      ! 2. Look for segment crossings 
     1010      !------------------------------  
     1011      IF (lwp) WRITE(numout,*) '**Look for segments corners  :' 
     1012      ! 
     1013      itest = 0 ! corner number 
     1014      ! 
     1015      ! flag to detect if start or end of open boundary belongs to a corner 
     1016      ! if not (=0), it must be on land. 
     1017      ! if a corner is detected, save bdy package number for further tests 
     1018      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0. 
     1019      ! South/West crossings 
     1020      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN 
     1021         DO ib1 = 1, nbdysegw         
     1022            DO ib2 = 1, nbdysegs 
     1023               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. & 
     1024                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. & 
     1025                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. & 
     1026                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN 
     1027                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN  
     1028                     ! We have a possible South-West corner                       
     1029!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)  
     1030!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2) 
     1031                     icornw(ib1,1) = npckgs(ib2) 
     1032                     icorns(ib2,1) = npckgw(ib1) 
     1033                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN 
     1034                     IF(lwp) WRITE(numout,*) 
     1035                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1036                     &                                     jpisft(ib2), jpjwft(ib1) 
     1037                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
     1038                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), &  
     1039                     &                                                    ' and South segment: ',npckgs(ib2) 
     1040                     IF(lwp) WRITE(numout,*) 
     1041                     nstop = nstop + 1 
     1042                  ELSE 
     1043                     IF(lwp) WRITE(numout,*) 
     1044                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices' 
     1045                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , & 
     1046                     &                                                    ' and South segment: ',npckgs(ib2) 
     1047                     IF(lwp) WRITE(numout,*) 
     1048                     nstop = nstop+1 
     1049                  END IF 
     1050               END IF 
     1051            END DO 
     1052         END DO 
     1053      END IF 
     1054      ! 
     1055      ! South/East crossings 
     1056      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN 
     1057         DO ib1 = 1, nbdysege 
     1058            DO ib2 = 1, nbdysegs 
     1059               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. & 
     1060                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. & 
     1061                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. & 
     1062                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN 
     1063                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN 
     1064                     ! We have a possible South-East corner  
     1065!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)  
     1066!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2) 
     1067                     icorne(ib1,1) = npckgs(ib2) 
     1068                     icorns(ib2,2) = npckge(ib1) 
     1069                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN 
     1070                     IF(lwp) WRITE(numout,*) 
     1071                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1072                     &                                     jpisdt(ib2), jpjeft(ib1) 
     1073                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
     1074                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), & 
     1075                     &                                                    ' and South segment: ',npckgs(ib2) 
     1076                     IF(lwp) WRITE(numout,*) 
     1077                     nstop = nstop + 1 
     1078                  ELSE 
     1079                     IF(lwp) WRITE(numout,*) 
     1080                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices' 
     1081                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
     1082                     &                                                    ' and South segment: ',npckgs(ib2) 
     1083                     IF(lwp) WRITE(numout,*) 
     1084                     nstop = nstop + 1 
     1085                  END IF 
     1086               END IF 
     1087            END DO 
     1088         END DO 
     1089      END IF 
     1090      ! 
     1091      ! North/West crossings 
     1092      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN 
     1093         DO ib1 = 1, nbdysegw         
     1094            DO ib2 = 1, nbdysegn 
     1095               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. & 
     1096                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. & 
     1097                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. & 
     1098                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN 
     1099                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN 
     1100                     ! We have a possible North-West corner  
     1101!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)  
     1102!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2) 
     1103                     icornw(ib1,2) = npckgn(ib2) 
     1104                     icornn(ib2,1) = npckgw(ib1) 
     1105                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN 
     1106                     IF(lwp) WRITE(numout,*) 
     1107                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1108                     &                                     jpinft(ib2), jpjwdt(ib1) 
     1109                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
     1110                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
     1111                     &                                                    ' and North segment: ',npckgn(ib2) 
     1112                     IF(lwp) WRITE(numout,*) 
     1113                     nstop = nstop + 1 
     1114                  ELSE 
     1115                     IF(lwp) WRITE(numout,*) 
     1116                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices' 
     1117                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), & 
     1118                     &                                                    ' and North segment: ',npckgn(ib2) 
     1119                     IF(lwp) WRITE(numout,*) 
     1120                     nstop = nstop + 1 
     1121                  END IF 
     1122               END IF 
     1123            END DO 
     1124         END DO 
     1125      END IF 
     1126      ! 
     1127      ! North/East crossings 
     1128      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN 
     1129         DO ib1 = 1, nbdysege         
     1130            DO ib2 = 1, nbdysegn 
     1131               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. & 
     1132                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. & 
     1133                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. & 
     1134                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN 
     1135                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN 
     1136                     ! We have a possible North-East corner  
     1137!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1) 
     1138!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2) 
     1139                     icorne(ib1,2) = npckgn(ib2) 
     1140                     icornn(ib2,2) = npckge(ib1) 
     1141                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN 
     1142                     IF(lwp) WRITE(numout,*) 
     1143                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
     1144                     &                                     jpindt(ib2), jpjedt(ib1) 
     1145                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
     1146                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), & 
     1147                     &                                                    ' and North segment: ',npckgn(ib2) 
     1148                     IF(lwp) WRITE(numout,*) 
     1149                     nstop = nstop + 1 
     1150                  ELSE 
     1151                     IF(lwp) WRITE(numout,*) 
     1152                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices' 
     1153                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
     1154                     &                                                    ' and North segment: ',npckgn(ib2) 
     1155                     IF(lwp) WRITE(numout,*) 
     1156                     nstop = nstop + 1 
     1157                  END IF 
     1158               END IF 
     1159            END DO 
     1160         END DO 
     1161      END IF 
     1162      ! 
     1163      ! 3. Check if segment extremities are on land 
     1164      !--------------------------------------------  
     1165      ! 
     1166      ! West segments 
     1167      DO ib = 1, nbdysegw 
     1168         ! get mask at boundary extremities: 
     1169         ztestmask(1:2)=0. 
     1170         DO ji = 1, jpi 
     1171            DO jj = 1, jpj              
     1172              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &  
     1173               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1) 
     1174              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &  
     1175               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)   
     1176            END DO 
     1177         END DO 
     1178         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
     1179 
     1180         IF (ztestmask(1)==1) THEN  
     1181            IF (icornw(ib,1)==0) THEN 
     1182               IF(lwp) WRITE(numout,*) 
     1183               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
     1184               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
     1185               IF(lwp) WRITE(numout,*) 
     1186               nstop = nstop + 1 
     1187            ELSE 
     1188               ! This is a corner 
     1189               WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib) 
     1190               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1)) 
     1191               itest=itest+1 
     1192            ENDIF 
     1193         ENDIF 
     1194         IF (ztestmask(2)==1) THEN 
     1195            IF (icornw(ib,2)==0) THEN 
     1196               IF(lwp) WRITE(numout,*) 
     1197               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
     1198               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
     1199               IF(lwp) WRITE(numout,*) 
     1200               nstop = nstop + 1 
     1201            ELSE 
     1202               ! This is a corner 
     1203               WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib) 
     1204               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2)) 
     1205               itest=itest+1 
     1206            ENDIF 
     1207         ENDIF 
     1208      END DO 
     1209      ! 
     1210      ! East segments 
     1211      DO ib = 1, nbdysege 
     1212         ! get mask at boundary extremities: 
     1213         ztestmask(1:2)=0. 
     1214         DO ji = 1, jpi 
     1215            DO jj = 1, jpj              
     1216              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &  
     1217               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1) 
     1218              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &  
     1219               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)   
     1220            END DO 
     1221         END DO 
     1222         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
     1223 
     1224         IF (ztestmask(1)==1) THEN 
     1225            IF (icorne(ib,1)==0) THEN 
     1226               IF(lwp) WRITE(numout,*) 
     1227               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib) 
     1228               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
     1229               IF(lwp) WRITE(numout,*) 
     1230               nstop = nstop + 1  
     1231            ELSE 
     1232               ! This is a corner 
     1233               WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib) 
     1234               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1)) 
     1235               itest=itest+1 
     1236            ENDIF 
     1237         ENDIF 
     1238         IF (ztestmask(2)==1) THEN 
     1239            IF (icorne(ib,2)==0) THEN 
     1240               IF(lwp) WRITE(numout,*) 
     1241               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib) 
     1242               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
     1243               IF(lwp) WRITE(numout,*) 
     1244               nstop = nstop + 1 
     1245            ELSE 
     1246               ! This is a corner 
     1247               WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib) 
     1248               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2)) 
     1249               itest=itest+1 
     1250            ENDIF 
     1251         ENDIF 
     1252      END DO 
     1253      ! 
     1254      ! South segments 
     1255      DO ib = 1, nbdysegs 
     1256         ! get mask at boundary extremities: 
     1257         ztestmask(1:2)=0. 
     1258         DO ji = 1, jpi 
     1259            DO jj = 1, jpj              
     1260              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &  
     1261               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1) 
     1262              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &  
     1263               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)   
     1264            END DO 
     1265         END DO 
     1266         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
     1267 
     1268         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN 
     1269            IF(lwp) WRITE(numout,*) 
     1270            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
     1271            IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
     1272            IF(lwp) WRITE(numout,*) 
     1273            nstop = nstop + 1 
     1274         ENDIF 
     1275         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN 
     1276            IF(lwp) WRITE(numout,*) 
     1277            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
     1278            IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
     1279            IF(lwp) WRITE(numout,*) 
     1280            nstop = nstop + 1 
     1281         ENDIF 
     1282      END DO 
     1283      ! 
     1284      ! North segments 
     1285      DO ib = 1, nbdysegn 
     1286         ! get mask at boundary extremities: 
     1287         ztestmask(1:2)=0. 
     1288         DO ji = 1, jpi 
     1289            DO jj = 1, jpj              
     1290              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &  
     1291               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1) 
     1292              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &  
     1293               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)   
     1294            END DO 
     1295         END DO 
     1296         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
     1297 
     1298         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN 
     1299            IF(lwp) WRITE(numout,*) 
     1300            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
     1301            IF(lwp) WRITE(numout,*) ' ==========  does not start on land'                                                   
     1302            IF(lwp) WRITE(numout,*) 
     1303            nstop = nstop + 1 
     1304         ENDIF 
     1305         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN 
     1306            IF(lwp) WRITE(numout,*) 
     1307            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
     1308            IF(lwp) WRITE(numout,*) ' ==========  does not end on land'                                                   
     1309            IF(lwp) WRITE(numout,*) 
     1310            nstop = nstop + 1 
     1311         ENDIF 
     1312      END DO 
     1313      ! 
     1314      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found' 
     1315      ! 
     1316      ! Other tests TBD:  
     1317      ! segments completly on land 
     1318      ! optimized open boundary array length according to landmask 
     1319      ! Nudging layers that overlap with interior domain 
     1320      ! 
     1321   END SUBROUTINE bdy_ctl_seg 
     1322 
     1323   SUBROUTINE bdy_ctl_corn( ib1, ib2 ) 
     1324      !!---------------------------------------------------------------------- 
     1325      !!                 ***  ROUTINE bdy_ctl_corn  *** 
     1326      !! 
     1327      !! ** Purpose :   Check numerical schemes consistency between 
     1328      !!                segments having a common corner 
     1329      !! 
     1330      !! ** Method  :    
     1331      !!---------------------------------------------------------------------- 
     1332      INTEGER, INTENT(in)  ::   ib1, ib2 
     1333      INTEGER :: itest 
     1334      !!---------------------------------------------------------------------- 
     1335      itest = 0 
     1336 
     1337      IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1 
     1338      IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1 
     1339      IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1 
     1340      ! 
     1341      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1 
     1342      IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1 
     1343      IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1 
     1344      ! 
     1345      IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1    
     1346      ! 
     1347      IF ( itest>0 ) THEN 
     1348         IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2 
     1349         IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                   
     1350         IF(lwp) WRITE(numout,*) 
     1351         nstop = nstop + 1 
     1352      ENDIF 
     1353      ! 
     1354   END SUBROUTINE bdy_ctl_corn 
    8071355 
    8081356#else 
Note: See TracChangeset for help on using the changeset viewer.