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 3490 for branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 – NEMO

Ignore:
Timestamp:
2012-10-08T16:27:20+02:00 (12 years ago)
Author:
cbricaud
Message:

add Jerome Chanut 's modications for BDY, Mercator_1 2012 task

File:
1 edited

Legend:

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

    r3367 r3490  
    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,             & 
     
    8695         &             ln_vol, nn_volctl, nn_rimwidth 
    8796      !! 
    88       NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             & 
    89                              nbdysegw, jpiwob, jpjwdt, jpjwft,             & 
    90                              nbdysegn, jpjnob, jpindt, jpinft,             & 
    91                              nbdysegs, jpjsob, jpisdt, jpisft 
     97      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
    9298 
    9399      !!---------------------------------------------------------------------- 
     
    106112 
    107113      cgrid= (/'t','u','v'/) 
    108  
     114       
    109115      ! ----------------------------------------- 
    110116      ! Initialise and read namelist parameters 
     
    139145      ! Check and write out namelist parameters 
    140146      ! ----------------------------------------- 
    141   
    142       ln_tra_dmp_tot=.false. 
    143       ln_dyn3d_dmp_tot=.false. 
    144147      !                                   ! control prints 
    145148      IF(lwp) WRITE(numout,*) '         nambdy' 
     
    164167        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
    165168        SELECT CASE( nn_dyn2d(ib_bdy) )                   
    166           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    167           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    168           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' 
    169172          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
    170173        END SELECT 
     
    177180              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
    178181           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 
    179185        ENDIF 
    180186        IF(lwp) WRITE(numout,*) 
     
    182188        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
    183189        SELECT CASE( nn_dyn3d(ib_bdy) )                   
    184           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    185           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' 
    186192          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    187193          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     
    195201           END SELECT 
    196202        ENDIF 
    197         IF(lwp) WRITE(numout,*) 
    198203 
    199204        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
    200205           IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN 
    201               IF(lwp) WRITE(numout,*) 'No open boundary condition for the baroclinic velocitues : we force ln_dyn3d_dmp=.false.' 
     206              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
    202207              ln_dyn3d_dmp(ib_bdy)=.false. 
     208           ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN 
     209              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    203210           ELSE 
    204               IF(lwp) WRITE(numout,*) '      Damping of the baroclinic velocities along the boundaries' 
    205               IF(lwp) WRITE(numout,*) '         Time damping :',rn_time_dmp(ib_bdy),' days' 
    206               ln_dyn3d_dmp_tot=.true. 
     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' ) 
    207214           ENDIF 
     215        ELSE 
     216           IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities' 
    208217        ENDIF 
    209218        IF(lwp) WRITE(numout,*) 
     
    211220        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
    212221        SELECT CASE( nn_tra(ib_bdy) )                   
    213           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    214           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' 
    215224          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    216225          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions' 
     
    225234           END SELECT 
    226235        ENDIF 
    227         IF(lwp) WRITE(numout,*) 
    228236 
    229237        IF ( ln_tra_dmp(ib_bdy) ) THEN 
    230238           IF ( nn_tra(ib_bdy).EQ.0 ) THEN 
    231               IF(lwp) WRITE(numout,*) 'No open boundary condition for the tracer : we force ln_tra_dmp=.false.' 
     239              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
    232240              ln_tra_dmp(ib_bdy)=.false. 
     241           ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN 
     242              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    233243           ELSE 
    234               IF(lwp) WRITE(numout,*) '      Damping of the scalars along the boundaries' 
    235               IF(lwp) WRITE(numout,*) '         Time dumping :',rn_time_dmp(ib_bdy),' days' 
    236               ln_tra_dmp_tot=.true. 
     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' ) 
    237247           ENDIF 
     248        ELSE 
     249           IF(lwp) WRITE(numout,*) '      NO T/S relaxation' 
    238250        ENDIF 
    239251        IF(lwp) WRITE(numout,*) 
     
    256268#endif 
    257269 
    258         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) 
    259271        IF(lwp) WRITE(numout,*) 
    260272 
    261273      ENDDO 
    262274 
    263      IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
    264        IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
    265        IF(lwp) WRITE(numout,*) 
    266        SELECT CASE ( nn_volctl ) 
    267          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    268          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    269          CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
    270        END SELECT 
    271        IF(lwp) WRITE(numout,*) 
    272      ELSE 
    273        IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
    274        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 
    275289     ENDIF 
    276290 
     
    282296      ! --------------------------------------------- 
    283297      REWIND( numnam )                     
    284       jpbdta = 1  
     298               
    285299      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 
    286309      DO ib_bdy = 1, nb_bdy 
    287310 
    288311         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
    289312  
     313            icount = icount + 1 
    290314            ! No REWIND here because may need to read more than one nambdy_index namelist. 
    291315            READ  ( numnam, nambdy_index ) 
    292316 
    293             ! Automatic boundary definition: if nbdysegX = -1 
    294             ! set boundary to whole side of model domain. 
    295             IF( nbdysege == -1 ) THEN 
    296                nbdysege = 1 
    297                jpieob(1) = jpiglo - 1 
    298                jpjedt(1) = 2 
    299                jpjeft(1) = jpjglo - 1 
    300             ENDIF 
    301             IF( nbdysegw == -1 ) THEN 
    302                nbdysegw = 1 
    303                jpiwob(1) = 2 
    304                jpjwdt(1) = 2 
    305                jpjwft(1) = jpjglo - 1 
    306             ENDIF 
    307             IF( nbdysegn == -1 ) THEN 
    308                nbdysegn = 1 
    309                jpjnob(1) = jpjglo - 1 
    310                jpindt(1) = 2 
    311                jpinft(1) = jpiglo - 1 
    312             ENDIF 
    313             IF( nbdysegs == -1 ) THEN 
    314                nbdysegs = 1 
    315                jpjsob(1) = 2 
    316                jpisdt(1) = 2 
    317                jpisft(1) = jpiglo - 1 
    318             ENDIF 
    319  
    320             DO iseg = 1, nbdysege 
    321                igrd = 1 
    322                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
    323                igrd = 2 
    324                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
    325                igrd = 3 
    326                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)                
    327             ENDDO 
    328             DO iseg = 1, nbdysegw 
    329                igrd = 1 
    330                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
    331                igrd = 2 
    332                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
    333                igrd = 3 
    334                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)                
    335             ENDDO 
    336             DO iseg = 1, nbdysegn 
    337                igrd = 1 
    338                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1                
    339                igrd = 2 
    340                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)                
    341                igrd = 3 
    342                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1 
    343             ENDDO 
    344             DO iseg = 1, nbdysegs 
    345                igrd = 1 
    346                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
    347                igrd = 2 
    348                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) 
    349                igrd = 3 
    350                nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
    351             ENDDO 
    352  
    353             nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 
     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' ) 
    354378 
    355379         ELSE            ! Read size of arrays in boundary coordinates file. 
    356  
    357  
    358380            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    359381            DO igrd = 1, jpbgrd 
    360382               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    361383               nblendta(igrd,ib_bdy) = kdimsz(1) 
    362             ENDDO 
     384               jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     385            ENDDO 
     386            CALL iom_close( inum ) 
    363387 
    364388         ENDIF  
     
    366390      ENDDO ! ib_bdy 
    367391 
    368       jpbdta = MAXVAL(nblendta(:,:)) 
    369  
    370       ! Allocate arrays 
    371       !--------------- 
    372       ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
    373          &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
    374  
    375       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 
    376408 
    377409      ! Calculate global boundary index arrays or read in from file 
    378       !------------------------------------------------------------ 
    379       REWIND( numnam )                     
     410      !------------------------------------------------------------                
     411      ! 1. Read global index arrays from boundary coordinates file. 
    380412      DO ib_bdy = 1, nb_bdy 
    381413 
    382          IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters 
    383  
    384             ! No REWIND here because may need to read more than one nambdy_index namelist. 
    385             READ  ( numnam, nambdy_index ) 
    386  
    387             ! Automatic boundary definition: if nbdysegX = -1 
    388             ! set boundary to whole side of model domain. 
    389             IF( nbdysege == -1 ) THEN 
    390                nbdysege = 1 
    391                jpieob(1) = jpiglo - 1 
    392                jpjedt(1) = 2 
    393                jpjeft(1) = jpjglo - 1 
    394             ENDIF 
    395             IF( nbdysegw == -1 ) THEN 
    396                nbdysegw = 1 
    397                jpiwob(1) = 2 
    398                jpjwdt(1) = 2 
    399                jpjwft(1) = jpjglo - 1 
    400             ENDIF 
    401             IF( nbdysegn == -1 ) THEN 
    402                nbdysegn = 1 
    403                jpjnob(1) = jpjglo - 1 
    404                jpindt(1) = 2 
    405                jpinft(1) = jpiglo - 1 
    406             ENDIF 
    407             IF( nbdysegs == -1 ) THEN 
    408                nbdysegs = 1 
    409                jpjsob(1) = 2 
    410                jpisdt(1) = 2 
    411                jpisft(1) = jpiglo - 1 
    412             ENDIF 
    413  
    414             ! ------------ T points ------------- 
    415             igrd = 1   
    416             icount = 0 
    417             DO ir = 1, nn_rimwidth(ib_bdy) 
    418                ! east 
    419                DO iseg = 1, nbdysege 
    420                   DO ij = jpjedt(iseg), jpjeft(iseg) 
    421                      icount = icount + 1 
    422                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
    423                      nbjdta(icount, igrd, ib_bdy) = ij 
    424                      nbrdta(icount, igrd, ib_bdy) = ir 
    425                   ENDDO 
    426                ENDDO 
    427                ! west 
    428                DO iseg = 1, nbdysegw 
    429                   DO ij = jpjwdt(iseg), jpjwft(iseg) 
    430                      icount = icount + 1 
    431                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    432                      nbjdta(icount, igrd, ib_bdy) = ij 
    433                      nbrdta(icount, igrd, ib_bdy) = ir 
    434                   ENDDO 
    435                ENDDO 
    436                ! north 
    437                DO iseg = 1, nbdysegn 
    438                   DO ii = jpindt(iseg), jpinft(iseg) 
    439                      icount = icount + 1 
    440                      nbidta(icount, igrd, ib_bdy) = ii 
    441                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
    442                      nbrdta(icount, igrd, ib_bdy) = ir 
    443                   ENDDO 
    444                ENDDO 
    445                ! south 
    446                DO iseg = 1, nbdysegs 
    447                   DO ii = jpisdt(iseg), jpisft(iseg) 
    448                      icount = icount + 1 
    449                      nbidta(icount, igrd, ib_bdy) = ii 
    450                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    451                      nbrdta(icount, igrd, ib_bdy) = ir 
    452                   ENDDO 
    453                ENDDO 
    454             ENDDO 
    455  
    456             ! ------------ U points ------------- 
    457             igrd = 2   
    458             icount = 0 
    459             DO ir = 1, nn_rimwidth(ib_bdy) 
    460                ! east 
    461                DO iseg = 1, nbdysege 
    462                   DO ij = jpjedt(iseg), jpjeft(iseg) 
    463                      icount = icount + 1 
    464                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir 
    465                      nbjdta(icount, igrd, ib_bdy) = ij 
    466                      nbrdta(icount, igrd, ib_bdy) = ir 
    467                   ENDDO 
    468                ENDDO 
    469                ! west 
    470                DO iseg = 1, nbdysegw 
    471                   DO ij = jpjwdt(iseg), jpjwft(iseg) 
    472                      icount = icount + 1 
    473                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    474                      nbjdta(icount, igrd, ib_bdy) = ij 
    475                      nbrdta(icount, igrd, ib_bdy) = ir 
    476                   ENDDO 
    477                ENDDO 
    478                ! north 
    479                DO iseg = 1, nbdysegn 
    480                   DO ii = jpindt(iseg), jpinft(iseg) - 1 
    481                      icount = icount + 1 
    482                      nbidta(icount, igrd, ib_bdy) = ii 
    483                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
    484                      nbrdta(icount, igrd, ib_bdy) = ir 
    485                   ENDDO 
    486                ENDDO 
    487                ! south 
    488                DO iseg = 1, nbdysegs 
    489                   DO ii = jpisdt(iseg), jpisft(iseg) - 1 
    490                      icount = icount + 1 
    491                      nbidta(icount, igrd, ib_bdy) = ii 
    492                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    493                      nbrdta(icount, igrd, ib_bdy) = ir 
    494                   ENDDO 
    495                ENDDO 
    496             ENDDO 
    497  
    498             ! ------------ V points ------------- 
    499             igrd = 3   
    500             icount = 0 
    501             DO ir = 1, nn_rimwidth(ib_bdy) 
    502                ! east 
    503                DO iseg = 1, nbdysege 
    504                   DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
    505                      icount = icount + 1 
    506                      nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
    507                      nbjdta(icount, igrd, ib_bdy) = ij 
    508                      nbrdta(icount, igrd, ib_bdy) = ir 
    509                   ENDDO 
    510                ENDDO 
    511                ! west 
    512                DO iseg = 1, nbdysegw 
    513                   DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
    514                      icount = icount + 1 
    515                      nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    516                      nbjdta(icount, igrd, ib_bdy) = ij 
    517                      nbrdta(icount, igrd, ib_bdy) = ir 
    518                   ENDDO 
    519                ENDDO 
    520                ! north 
    521                DO iseg = 1, nbdysegn 
    522                   DO ii = jpindt(iseg), jpinft(iseg) 
    523                      icount = icount + 1 
    524                      nbidta(icount, igrd, ib_bdy) = ii 
    525                      nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir 
    526                      nbrdta(icount, igrd, ib_bdy) = ir 
    527                   ENDDO 
    528                ENDDO 
    529                ! south 
    530                DO iseg = 1, nbdysegs 
    531                   DO ii = jpisdt(iseg), jpisft(iseg) 
    532                      icount = icount + 1 
    533                      nbidta(icount, igrd, ib_bdy) = ii 
    534                      nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    535                      nbrdta(icount, igrd, ib_bdy) = ir 
    536                   ENDDO 
    537                ENDDO 
    538             ENDDO 
    539  
    540          ELSE            ! Read global index arrays from boundary coordinates file. 
    541  
     414         IF( ln_coords_file(ib_bdy) ) THEN 
     415 
     416            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    542417            DO igrd = 1, jpbgrd 
    543418               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     
    560435               IF (ibr_max < nn_rimwidth(ib_bdy))   & 
    561436                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
    562  
    563437            END DO 
    564438            CALL iom_close( inum ) 
     
    566440         ENDIF  
    567441 
    568       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 
    569660 
    570661      ! Work out dimensions of boundary data on each processor 
    571662      ! ------------------------------------------------------ 
    572       
    573       iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
    574       ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
    575       is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
    576       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 
    577674 
    578675      DO ib_bdy = 1, nb_bdy 
     
    610707         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
    611708         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 
     709         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) ) 
    612710         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 
    613711         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
     
    629727                     ! 
    630728                     icount = icount  + 1 
    631                      idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
    632                      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 
    633736                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy) 
    634737                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 
     
    643746            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    644747               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    645 !               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
    646                idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.      ! quadratic 
    647 !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear 
     748               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
     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 
    648761            END DO 
    649762         END DO  
     
    660773      !          = 0  elsewhere    
    661774  
    662       IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution 
    663          zmask(         :                ,:) = 0.e0 
    664          zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0           
    665       ELSE IF( ln_mask_file ) THEN 
     775      IF( ln_mask_file ) THEN 
    666776         CALL iom_open( cn_mask_file, inum ) 
    667          CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) ) 
     777         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) ) 
    668778         CALL iom_close( inum ) 
    669       ELSE 
    670          zmask(:,:) = 1.e0 
    671       ENDIF 
    672  
    673       DO ij = 1, nlcj      ! Save mask over local domain       
    674          DO ii = 1, nlci 
    675             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 
    676788         END DO 
    677       END DO 
    678  
    679       ! Derive mask on U and V grid from mask on T grid 
    680       bdyumask(:,:) = 0.e0 
    681       bdyvmask(:,:) = 0.e0 
    682       DO ij=1, jpjm1 
    683          DO ii=1, jpim1 
    684             bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 
    685             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 
    686803         END DO 
    687       END DO 
    688       CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond. 
    689  
    690  
    691       ! Mask corrections 
    692       ! ---------------- 
    693       DO ik = 1, jpkm1 
    694          DO ij = 1, jpj 
    695             DO ii = 1, jpi 
    696                tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij) 
    697                umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij) 
    698                vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij) 
    699                bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij) 
    700             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 
    701812         END DO 
    702       END DO 
    703  
    704       DO ik = 1, jpkm1 
    705          DO ij = 2, jpjm1 
    706             DO ii = 2, jpim1 
    707                fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   & 
    708                   &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1) 
    709             END DO       
    710          END DO 
    711       END DO 
    712  
    713       tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)              
     813 
     814         tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 
     815 
     816      ENDIF ! ln_mask_file=.TRUE. 
     817       
    714818      bdytmask(:,:) = tmask(:,:,1) 
    715819 
     
    800904      ! Compute total lateral surface for volume correction: 
    801905      ! ---------------------------------------------------- 
     906      ! JC: this must be done at each time step with key_vvl 
    802907      bdysurftot = 0.e0  
    803908      IF( ln_vol ) THEN   
     
    806911            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    807912               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    808                nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     913               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    809914               flagu => idx_bdy(ib_bdy)%flagu(ib) 
    810915               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           & 
     
    819924            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    820925               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    821                nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     926               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    822927               flagv => idx_bdy(ib_bdy)%flagv(ib) 
    823928               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           & 
     
    833938      ! Tidy up 
    834939      !-------- 
    835       DEALLOCATE(nbidta, nbjdta, nbrdta) 
     940      IF (nb_bdy>0) THEN 
     941         DEALLOCATE(nbidta, nbjdta, nbrdta) 
     942      ENDIF 
    836943 
    837944      IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 
    838945 
    839946   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 
    8401355 
    8411356#else 
Note: See TracChangeset for help on using the changeset viewer.