Changeset 3490


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

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

Location:
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC
Files:
10 edited

Legend:

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

    r3367 r3490  
    2828      INTEGER, POINTER, DIMENSION(:,:)   ::  nbmap 
    2929      REAL   , POINTER, DIMENSION(:,:)   ::  nbw 
     30      REAL   , POINTER, DIMENSION(:,:)   ::  nbd 
    3031      REAL   , POINTER, DIMENSION(:)     ::  flagu 
    3132      REAL   , POINTER, DIMENSION(:)     ::  flagv 
     
    7374   INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta             !: = 0 use the initial state as bdy dta ;  
    7475                                                            !: = 1 read it in a NetCDF file 
    75    LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp 
    76    LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp 
    77    REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp 
    78    LOGICAL                    ::   ln_tra_dmp_tot 
    79    LOGICAL                    ::   ln_dyn3d_dmp_tot 
     76   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping 
     77   LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping 
     78   REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
    8079 
    8180#if defined key_lim2 
     
    107106   INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions 
    108107                                                                          !: =1 => some data to be read in from data files 
    109    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays 
     108   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy) 
     109   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
    110110   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
    111111   TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process) 
  • branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r3367 r3490  
    192192            IF( PRESENT(jit) ) THEN 
    193193               ! Update barotropic boundary conditions only 
    194                ! jit is optional argument for fld_read and tide_update 
     194               ! jit is optional argument for fld_read and bdytide_update 
    195195               IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
    196196                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     
    200200                  ENDIF 
    201201                  IF (nn_tra(ib_bdy).ne.4) THEN 
    202                    IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN 
    203                      ! For the runoff case, no need to update the forcing (already done in the baroclinic part) 
    204                      jend = nb_bdy_fld(ib_bdy) 
    205                      IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2 
    206                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset ) 
    207                      IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2 
    208                      ! If full velocities in boundary data then split into barotropic and baroclinic data 
    209                      IF( ln_full_vel_array(ib_bdy) .and.                                             & 
    210                     &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 
    211                         igrd = 2                      ! zonal velocity 
    212                         dta_bdy(ib_bdy)%u2d(:) = 0.0 
    213                         DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    214                            ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    215                            ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    216                            DO ik = 1, jpkm1 
    217                               dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
    218                     &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
     202                     IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN 
     203                        ! For the runoff case, no need to update the forcing (already done in the baroclinic part) 
     204                        jend = nb_bdy_fld(ib_bdy) 
     205                        IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2 
     206                        CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset ) 
     207                        IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2 
     208                        ! If full velocities in boundary data then split into barotropic and baroclinic data 
     209                        IF( ln_full_vel_array(ib_bdy) .and.                                             & 
     210                       &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 
     211                           igrd = 2                      ! zonal velocity 
     212                           dta_bdy(ib_bdy)%u2d(:) = 0.0 
     213                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     214                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     215                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     216                              DO ik = 1, jpkm1 
     217                                 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
     218                       &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
     219                              END DO 
     220                              dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
     221                              DO ik = 1, jpkm1 
     222                                 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 
     223                              END DO 
    219224                           END DO 
    220                            dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
    221                            DO ik = 1, jpkm1 
    222                               dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 
     225                           igrd = 3                      ! meridional velocity 
     226                           dta_bdy(ib_bdy)%v2d(:) = 0.0 
     227                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     228                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     229                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     230                              DO ik = 1, jpkm1 
     231                                 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
     232                       &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
     233                              END DO 
     234                              dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
     235                              DO ik = 1, jpkm1 
     236                                 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 
     237                              END DO 
    223238                           END DO 
    224                         END DO 
    225                         igrd = 3                      ! meridional velocity 
    226                         dta_bdy(ib_bdy)%v2d(:) = 0.0 
    227                         DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    228                            ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    229                            ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    230                            DO ik = 1, jpkm1 
    231                               dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
    232                     &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
    233                            END DO 
    234                            dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
    235                            DO ik = 1, jpkm1 
    236                               dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 
    237                            END DO 
    238                         END DO 
    239                      ENDIF                     
    240                    ENDIF 
    241                    IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    242                       CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   &  
    243                       &                 jit=jit, time_offset=time_offset ) 
    244                    ENDIF 
     239                        ENDIF                     
     240                     ENDIF 
     241                     IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     242                        CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   &  
     243                          &                 jit=jit, time_offset=time_offset ) 
     244                     ENDIF 
    245245                  ENDIF 
    246246               ENDIF 
     
    306306                  ENDIF 
    307307                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    308                      CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 
     308                     CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 
    309309                  ENDIF 
    310310               ENDIF 
     
    378378      ! Set nn_dta 
    379379      DO ib_bdy = 1, nb_bdy 
    380          IF (nn_dyn2d_dta(ib_bdy).eq.1.OR.nn_dyn2d_dta(ib_bdy).eq.3) nn_dta(ib_bdy) = 1 
    381          IF (nn_dyn3d_dta(ib_bdy).ge.1) nn_dta(ib_bdy) = 1 
    382          IF (nn_tra_dta  (ib_bdy).ge.1) nn_dta(ib_bdy) = 1 
    383 #if defined key_lim2 
    384          IF (nn_ice_lim2_dta(ib_bdy).eq.1) nn_dta(ib_bdy) = 1 
    385 #endif 
    386          IF (lwp) THEN 
    387             WRITE(numout,*) 'Bdy package number ',ib_bdy 
    388             IF (nn_dta(ib_bdy).eq.1) THEN 
    389                WRITE(numout,*) 'We use external data' 
    390             ELSE 
    391                WRITE(numout,*) 'We do not use external data' 
    392             ENDIF 
    393          ENDIF 
     380         nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       & 
     381                               ,nn_dyn3d_dta(ib_bdy)       & 
     382                               ,nn_tra_dta(ib_bdy)         & 
     383#if defined key_lim2 
     384                               ,nn_ice_lim2_dta(ib_bdy)    & 
     385#endif 
     386                              ) 
     387         IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1 
    394388      END DO 
    395389 
  • branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r3367 r3490  
    1414   !!---------------------------------------------------------------------- 
    1515   USE timing          ! Timing 
     16   USE wrk_nemo        ! Memory Allocation 
    1617   USE oce             ! ocean dynamics and tracers  
    1718   USE dom_oce         ! ocean space and time domain 
     
    2728   PUBLIC   bdy_dyn3d_dmp ! routine called by step 
    2829 
     30   !! * Substitutions 
     31#  include "domzgr_substitute.h90" 
    2932   !!---------------------------------------------------------------------- 
    3033   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     
    222225      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp') 
    223226      ! 
     227      !------------------------------------------------------- 
     228      ! Remove barotropic part from before velocity 
     229      !------------------------------------------------------- 
     230      CALL wrk_alloc(jpi,jpj,pu2d,pv2d)  
     231 
     232      pu2d(:,:) = 0.e0 
     233      pv2d(:,:) = 0.e0 
     234 
     235      DO jk = 1, jpkm1 
     236#if defined key_vvl 
     237         pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk)  
     238         pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk) 
     239#else 
     240         pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk) 
     241         pv2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk) 
     242#endif 
     243      END DO 
     244 
     245      IF( lk_vvl ) THEN 
     246         pu2d(:,:) = pu2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
     247         pv2d(:,:) = pv2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
     248      ELSE 
     249         pu2d(:,:) = pv2d(:,:) * hur(:,:) 
     250         pv2d(:,:) = pu2d(:,:) * hvr(:,:) 
     251      ENDIF 
     252 
    224253      DO ib_bdy=1, nb_bdy 
    225254         IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN 
     
    228257               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    229258               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    230                zwgt = idx_bdy(ib_bdy)%nbw(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday ) 
     259               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 
    231260               DO jk = 1, jpkm1 
    232                   ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) ) ) * umask(ii,ij,jk) 
     261                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) + pu2d(ii,ij)) ) * umask(ii,ij,jk) 
    233262               END DO 
    234263            END DO 
     
    238267               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    239268               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    240                zwgt = idx_bdy(ib_bdy)%nbw(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday ) 
     269               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 
    241270               DO jk = 1, jpkm1 
    242                   va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) ) ) * vmask(ii,ij,jk) 
     271                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) + pv2d(ii,ij)) ) * vmask(ii,ij,jk) 
    243272               END DO 
    244273            END DO 
    245274         ENDIF 
    246275      ENDDO 
     276      ! 
     277      CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
    247278      ! 
    248279      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 
  • 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 
  • branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r3367 r3490  
    88   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
    99   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes 
     10   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_bdy 
     
    2728   USE bdy_oce         ! ocean open boundary conditions 
    2829   USE daymod          ! calendar 
     30   USE wrk_nemo        ! Memory allocation 
    2931   USE tideini 
    30    USE tide_mod 
     32!   USE tide_mod       ! Useless ?? 
    3133   USE fldread, ONLY: fld_map 
    3234 
     
    3436   PRIVATE 
    3537 
    36    PUBLIC   bdytide_init  ! routine called in nemo_init 
    37    PUBLIC   tide_update   ! routine called in bdydyn 
     38   PUBLIC   bdytide_init     ! routine called in bdy_init 
     39   PUBLIC   bdytide_update   ! routine called in bdy_dta 
    3840 
    3941   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
     
    4749 
    4850   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data 
    49     
     51 
    5052   !!---------------------------------------------------------------------- 
    5153   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     
    6668      !!------------------- 
    6769      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files 
    68       LOGICAL                                   ::   ln_conjug = .FALSE. !: F/T : tidal data in complex/complex conjugate 
     70      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data 
     71      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data 
    6972      !! 
    7073      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices 
     74      INTEGER                                   ::   ii, ij              !: dummy loop indices 
    7175      INTEGER                                   ::   inum, igrd 
    72       INTEGER, DIMENSION(3)                     ::   ilen0               !: length of boundary data (from OBC arrays) 
     76      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays) 
     77      INTEGER, POINTER, DIMENSION(:)            ::   nblen, nblenrim     ! short cuts 
    7378      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file  
    7479      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data 
     80      REAL(wp), POINTER, DIMENSION(:,:)         ::   ztr, zti            !:  "     "    "   "   "   "        "      "  
    7581      !! 
    7682      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut    
    7783      !! 
    78       NAMELIST/nambdy_tide/filtide, ln_conjug 
     84      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 
    7985      !!---------------------------------------------------------------------- 
    8086 
    8187      IF( nn_timing == 1 ) CALL timing_start('bdytide_init') 
    8288 
    83       IF(lwp) WRITE(numout,*) 
    84       IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 
    85       IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
     89      IF (nb_bdy>0) THEN 
     90         IF(lwp) WRITE(numout,*) 
     91         IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 
     92         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     93      ENDIF 
     94 
     95      ln_bdytide_2ddta = .FALSE. 
     96      ln_bdytide_conj  = .FALSE. 
    8697 
    8798      REWIND(numnam) 
     
    90101 
    91102            td => tides(ib_bdy) 
     103            nblen => idx_bdy(ib_bdy)%nblen 
     104            nblenrim => idx_bdy(ib_bdy)%nblenrim 
    92105 
    93106            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 
     
    97110            READ  ( numnam, nambdy_tide ) 
    98111            !                                               ! Parameter control and print 
     112            IF(lwp) WRITE(numout,*) '  ' 
    99113            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
    100             IF(lwp) WRITE(numout,*) '             tidal components specified ', nb_harmo 
    101             IF(lwp) WRITE(numout,*) '                ', Wave(ntide(1:nb_harmo))%cname_tide 
    102             IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
    103             IF(lwp) WRITE(numout,*) '                ', omega_tide(1:nb_harmo) 
     114            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta 
     115            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj 
     116            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo 
     117            IF(lwp) THEN  
     118                    WRITE(numout,*) '             Tidal cpt name    -     Phase speed (deg/hr)'             
     119               DO itide = 1, nb_harmo 
     120                  WRITE(numout,*)  '             ', Wave(ntide(itide))%cname_tide, omega_tide(itide)    
     121               END DO 
     122            ENDIF  
     123            IF(lwp) WRITE(numout,*) ' ' 
    104124 
    105125            ! Allocate space for tidal harmonics data - get size from OBC data arrays 
    106126            ! ----------------------------------------------------------------------- 
    107127 
    108             ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh )  
     128            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 
     129            ! relaxation area       
     130            IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     131               ilen0(:)=nblen(:) 
     132            ELSE 
     133               ilen0(:)=nblenrim(:) 
     134            ENDIF 
     135 
    109136            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) ) 
    110137            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) ) 
    111138 
    112             ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d )  
    113139            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) ) 
    114140            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) ) 
    115141 
    116             ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d )  
    117142            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) ) 
    118143            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 
    119144 
    120             ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 
    121  
    122             ! Open files and read in tidal forcing data 
    123             ! ----------------------------------------- 
    124  
    125             DO itide = 1, nb_harmo 
    126                !                                                              ! SSH fields 
    127                clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 
    128                CALL iom_open( clfile, inum ) 
    129                CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
    130                td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 
    131                CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
    132                td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 
     145            td%ssh0(:,:,:) = 0.e0 
     146            td%ssh(:,:,:) = 0.e0 
     147            td%u0(:,:,:) = 0.e0 
     148            td%u(:,:,:) = 0.e0 
     149            td%v0(:,:,:) = 0.e0 
     150            td%v(:,:,:) = 0.e0 
     151 
     152            IF (ln_bdytide_2ddta) THEN 
     153               ! It is assumed that each data file contains all complex harmonic amplitudes 
     154               ! given on the data domain (ie global, jpidta x jpjdta) 
     155               ! 
     156               CALL wrk_alloc( jpi, jpj, zti, ztr ) 
     157               ! 
     158               ! SSH fields 
     159               clfile = TRIM(filtide)//'_grid_T.nc' 
     160               CALL iom_open (clfile , inum )  
     161               igrd = 1                       ! Everything is at T-points here 
     162               DO itide = 1, nb_harmo 
     163                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) ) 
     164                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )  
     165                  DO ib = 1, ilen0(igrd) 
     166                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     167                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     168                     td%ssh0(ib,itide,1) = ztr(ii,ij) 
     169                     td%ssh0(ib,itide,2) = zti(ii,ij) 
     170                  END DO 
     171               END DO  
    133172               CALL iom_close( inum ) 
    134                !                                                              ! U fields 
    135                clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 
    136                CALL iom_open( clfile, inum ) 
    137                CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
    138                td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 
    139                CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
    140                td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 
     173               ! 
     174               ! U fields 
     175               clfile = TRIM(filtide)//'_grid_U.nc' 
     176               CALL iom_open (clfile , inum )  
     177               igrd = 2                       ! Everything is at U-points here 
     178               DO itide = 1, nb_harmo 
     179                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) ) 
     180                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) ) 
     181                  DO ib = 1, ilen0(igrd) 
     182                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     183                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     184                     td%u0(ib,itide,1) = ztr(ii,ij) 
     185                     td%u0(ib,itide,2) = zti(ii,ij) 
     186                  END DO 
     187               END DO 
    141188               CALL iom_close( inum ) 
    142                !                                                              ! V fields 
    143                clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 
    144                CALL iom_open( clfile, inum ) 
    145                CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
    146                td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 
    147                CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
    148                td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 
     189               ! 
     190               ! V fields 
     191               clfile = TRIM(filtide)//'_grid_V.nc' 
     192               CALL iom_open (clfile , inum )  
     193               igrd = 3                       ! Everything is at V-points here 
     194               DO itide = 1, nb_harmo 
     195                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) ) 
     196                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) ) 
     197                  DO ib = 1, ilen0(igrd) 
     198                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     199                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     200                     td%v0(ib,itide,1) = ztr(ii,ij) 
     201                     td%v0(ib,itide,2) = zti(ii,ij) 
     202                  END DO 
     203               END DO   
    149204               CALL iom_close( inum ) 
    150                IF ( ln_conjug ) THEN 
    151                   IF(lwp) WRITE(numout,*) '       The tidal input data are written in complex conjugate' 
    152                   td%ssh0(:,itide,2) = - td%ssh0(:,itide,2) 
    153                   td%u0  (:,itide,2) = - td%u0  (:,itide,2) 
    154                   td%v0  (:,itide,2) = - td%v0  (:,itide,2) 
    155                ENDIF 
    156                ! 
    157             END DO ! end loop on tidal components 
     205               ! 
     206               CALL wrk_dealloc( jpi, jpj, ztr, zti )  
     207               ! 
     208            ELSE             
     209               ! 
     210               ! Read tidal data only on bdy segments 
     211               !  
     212               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) ) 
     213 
     214               ! Open files and read in tidal forcing data 
     215               ! ----------------------------------------- 
     216 
     217               DO itide = 1, nb_harmo 
     218                  !                                                              ! SSH fields 
     219                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 
     220                  CALL iom_open( clfile, inum ) 
     221                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     222                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 
     223                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     224                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 
     225                  CALL iom_close( inum ) 
     226                  !                                                              ! U fields 
     227                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 
     228                  CALL iom_open( clfile, inum ) 
     229                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     230                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 
     231                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     232                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 
     233                  CALL iom_close( inum ) 
     234                  !                                                              ! V fields 
     235                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 
     236                  CALL iom_open( clfile, inum ) 
     237                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     238                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 
     239                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     240                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 
     241                  CALL iom_close( inum ) 
     242                  ! 
     243               END DO ! end loop on tidal components 
     244               ! 
     245               DEALLOCATE( dta_read ) 
     246            ENDIF ! ln_bdytide_2ddta=.true. 
    158247            ! 
    159             DEALLOCATE( dta_read ) 
     248            IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files 
     249               td%ssh0(:,:,2) = - td%ssh0(:,:,2) 
     250               td%u0  (:,:,2) = - td%u0  (:,:,2) 
     251               td%v0  (:,:,2) = - td%v0  (:,:,2) 
     252            ENDIF 
    160253            ! 
    161254         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
     
    167260   END SUBROUTINE bdytide_init 
    168261 
    169    SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) 
    170       !!---------------------------------------------------------------------- 
    171       !!                 ***  SUBROUTINE tide_update  *** 
     262   SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset ) 
     263      !!---------------------------------------------------------------------- 
     264      !!                 ***  SUBROUTINE bdytide_update  *** 
    172265      !!                 
    173266      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
     
    186279                                                        ! etc. 
    187280      !! 
     281      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays) 
    188282      INTEGER                          :: itide, igrd, ib   ! dummy loop indices 
    189283      INTEGER                          :: time_add          ! time offset in units of timesteps 
    190       REAL(wp)                         :: z_arg, z_sarg, zflag       
     284      REAL(wp)                         :: z_arg, z_sarg, zflag, zramp       
    191285      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost 
    192286      !!---------------------------------------------------------------------- 
    193287 
    194       IF( nn_timing == 1 ) CALL timing_start('tide_update') 
     288      IF( nn_timing == 1 ) CALL timing_start('bdytide_update') 
     289 
     290      ilen0(1) =  SIZE(td%ssh(:,1,1)) 
     291      ilen0(2) =  SIZE(td%u(:,1,1)) 
     292      ilen0(3) =  SIZE(td%v(:,1,1)) 
    195293 
    196294      zflag=1 
     
    205303        IF(lwp) THEN 
    206304           WRITE(numout,*) 
    207            WRITE(numout,*) 'bdy_tide : (re)Initialization of the tidal bdy forcing at kt=',kt 
    208            WRITE(numout,*) '~~~~~~~ ' 
     305           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt 
     306           WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
    209307        ENDIF 
    210308        ! 
     
    225323      ENDIF 
    226324 
     325      ! Linear ramp on tidal component at open boundaries  
     326      zramp = 1. 
     327      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
     328 
    227329      DO itide = 1, nb_harmo 
    228330         z_sarg = z_arg * omega_tide(itide) 
     
    233335      DO itide = 1, nb_harmo 
    234336         igrd=1                              ! SSH on tracer grid 
    235          DO ib = 1, idx%nblenrim(igrd) 
    236             dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) 
     337         DO ib = 1, ilen0(igrd) 
     338            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)) 
    237339         END DO 
    238340         igrd=2                              ! U grid 
    239          DO ib=1, idx%nblenrim(igrd) 
    240             dta%u2d(ib) = dta%u2d(ib) + td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide) 
     341         DO ib = 1, ilen0(igrd) 
     342            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide)) 
    241343         END DO 
    242344         igrd=3                              ! V grid 
    243          DO ib=1, idx%nblenrim(igrd) 
    244             dta%v2d(ib) = dta%v2d(ib) + td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide) 
     345         DO ib = 1, ilen0(igrd)  
     346            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide)) 
    245347         END DO 
    246348      END DO 
    247349      ! 
    248       IF( nn_timing == 1 ) CALL timing_stop('tide_update') 
     350      IF( nn_timing == 1 ) CALL timing_stop('bdytide_update') 
    249351      ! 
    250    END SUBROUTINE tide_update 
     352   END SUBROUTINE bdytide_update 
    251353 
    252354   SUBROUTINE tide_init_elevation( idx, td ) 
     
    257359      TYPE(TIDES_DATA),INTENT( inout )   ::   td      ! tidal harmonics data 
    258360      !! * Local declarations 
     361      INTEGER, DIMENSION(1)            ::   ilen0       !: length of boundary data (from OBC arrays) 
    259362      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide 
    260363      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices 
    261364 
    262       igrd=1                                 ! SSH on tracer grid. 
    263  
    264       ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd))) 
     365      igrd=1    
     366                              ! SSH on tracer grid. 
     367    
     368      ilen0(1) =  SIZE(td%ssh0(:,1,1)) 
     369 
     370      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 
    265371 
    266372      DO itide = 1, nb_harmo 
    267          DO ib = 1, idx%nblenrim(igrd) 
     373         DO ib = 1, ilen0(igrd) 
    268374            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.) 
    269375            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1)) 
    270376         END DO 
    271          DO ib = 1, idx%nblenrim(igrd) 
     377         DO ib = 1 , ilen0(igrd) 
    272378            mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    273379            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
    274380         ENDDO 
    275          DO ib = 1, idx%nblenrim(igrd) 
     381         DO ib = 1 , ilen0(igrd) 
    276382            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
    277383            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     
    290396      TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data 
    291397      !! * Local declarations 
     398      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays) 
    292399      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide 
    293400      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices 
    294401 
     402      ilen0(2) =  SIZE(td%u0(:,1,1)) 
     403      ilen0(3) =  SIZE(td%v0(:,1,1)) 
     404 
    295405      igrd=2                                 ! U grid. 
    296406 
    297       ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd))) 
     407      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 
    298408 
    299409      DO itide = 1, nb_harmo 
    300          DO ib = 1, idx%nblenrim(igrd) 
     410         DO ib = 1, ilen0(igrd) 
    301411            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.) 
    302412            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1)) 
    303413         END DO 
    304          DO ib = 1, idx%nblenrim(igrd) 
     414         DO ib = 1, ilen0(igrd) 
    305415            mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    306416            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
    307417         ENDDO 
    308          DO ib = 1, idx%nblenrim(igrd) 
     418         DO ib = 1, ilen0(igrd) 
    309419            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
    310420            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     
    314424      DEALLOCATE(mod_tide,phi_tide) 
    315425 
    316       igrd=3                                 ! U grid. 
    317  
    318       ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd))) 
     426      igrd=3                                 ! V grid. 
     427 
     428      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 
    319429 
    320430      DO itide = 1, nb_harmo 
    321          DO ib = 1, idx%nblenrim(igrd) 
     431         DO ib = 1, ilen0(igrd) 
    322432            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.) 
    323433            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1)) 
    324434         END DO 
    325          DO ib = 1, idx%nblenrim(igrd) 
     435         DO ib = 1, ilen0(igrd) 
    326436            mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    327437            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
    328438         ENDDO 
    329          DO ib = 1, idx%nblenrim(igrd) 
     439         DO ib = 1, ilen0(igrd) 
    330440            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 
    331441            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 
     
    341451   !!---------------------------------------------------------------------- 
    342452CONTAINS 
    343    SUBROUTINE bdytide_init                ! Empty routine 
     453   SUBROUTINE bdytide_init             ! Empty routine 
     454      WRITE(*,*) 'bdytide_init: You should not have seen this print! error?' 
    344455   END SUBROUTINE bdytide_init 
    345    SUBROUTINE tide_data                ! Empty routine 
    346    END SUBROUTINE tide_data 
    347    SUBROUTINE tide_update( kt, jit )   ! Empty routine 
    348       WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit 
    349    END SUBROUTINE tide_update 
     456   SUBROUTINE bdytide_update( kt, jit )   ! Empty routine 
     457      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 
     458   END SUBROUTINE bdytide_update 
    350459#endif 
    351460 
  • branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r3367 r3490  
    2222   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2323   USE in_out_manager  ! I/O manager 
    24    USE phycst 
     24 
    2525 
    2626   IMPLICIT NONE 
     
    2828 
    2929   PUBLIC bdy_tra      ! routine called in tranxt.F90  
    30    PUBLIC bdy_tra_dmp  ! routine called in tranxt.F90  
     30   PUBLIC bdy_tra_dmp  ! routine called in step.F90  
    3131 
    3232   !!---------------------------------------------------------------------- 
     
    6565         END SELECT 
    6666      ENDDO 
     67      ! 
     68      ! Boundary points should be updated 
     69      IF (nb_bdy>0) CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) 
     70      IF (nb_bdy>0) CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 
    6771 
    6872   END SUBROUTINE bdy_tra 
     
    98102      END DO  
    99103      ! 
    100       CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated 
    101       ! 
    102104      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
    103105      ! 
     
    133135         END DO 
    134136      END DO 
    135       ! 
    136       CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated 
    137137      ! 
    138138      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     
    190190      END DO 
    191191      ! 
    192       CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated 
    193       ! 
    194192      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
    195193      ! 
     
    229227         END DO 
    230228      END DO 
    231       ! 
    232       CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated 
    233229      ! 
    234230      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     
    262258               ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    263259               ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    264                zwgt = idx_bdy(ib_bdy)%nbw(ib,igrd) / ( rn_time_dmp(ib_bdy) * rday ) 
     260               zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd) 
    265261               DO ik = 1, jpkm1 
    266                   tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) + zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) ) * tmask(ii,ij,ik) 
    267                   tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) + zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) ) * tmask(ii,ij,ik) 
     262                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) * tmask(ii,ij,ik) 
     263                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) * tmask(ii,ij,ik) 
     264                  tsa(ii,ij,ik,jp_tem) = tsa(ii,ij,ik,jp_tem) + zta 
     265                  tsa(ii,ij,ik,jp_sal) = tsa(ii,ij,ik,jp_sal) + zsa 
    268266               END DO 
    269267            END DO 
    270268         ENDIF 
    271269      ENDDO 
    272       ! 
    273       CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )  ;  CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated 
    274270      ! 
    275271      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp') 
  • branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90

    r3367 r3490  
    2626       ftide 
    2727 
    28   LOGICAL, PUBLIC :: ln_tide_pot = .false. 
     28  LOGICAL, PUBLIC :: ln_tide_pot = .false., ln_tide_ramp = .false. 
     29  REAL(wp), PUBLIC :: rdttideramp  
    2930  INTEGER, PUBLIC :: nb_harmo 
    3031  INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide 
     
    4647    CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 
    4748    ! 
    48     NAMELIST/nam_tide/ln_tide_pot,clname 
     49    NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname 
    4950    !!---------------------------------------------------------------------- 
    5051 
     
    5455          WRITE(numout,*) 
    5556          WRITE(numout,*) 'tide_init : Initialization of the tidal components' 
    56           WRITE(numout,*) '~~~~~~~ ' 
     57          WRITE(numout,*) '~~~~~~~~~ ' 
    5758       ENDIF 
    5859       ! 
     
    7778          WRITE(numout,*) '        Namelist nam_tide' 
    7879          WRITE(numout,*) '        nb_harmo    = ', nb_harmo 
     80          WRITE(numout,*) '        ln_tide_ramp = ', ln_tide_ramp  
     81          WRITE(numout,*) '        rdttideramp = ', rdttideramp 
     82          IF (ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp)) & 
     83          & CALL ctl_stop('rdttideramp must be lower than run duration') 
     84          IF (ln_tide_ramp.AND.(rdttideramp<0.)) & 
     85          & CALL ctl_stop('rdttideramp must be positive') 
    7986          CALL flush(numout) 
    8087       ENDIF 
  • branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90

    r3294 r3490  
    1313  USE sbctide 
    1414  USE dynspg_oce 
     15  USE tideini, ONLY: ln_tide_ramp, rdttideramp 
    1516 
    1617  IMPLICIT NONE 
     
    3334    INTEGER, INTENT( in ) ::   kt,kit      ! ocean time-step index 
    3435    INTEGER  :: ji,jj,jk 
     36    REAL (wp) :: zramp 
    3537    REAL (wp), DIMENSION(nb_harmo) :: zwt  
    36     !............................................................................... 
    37     ! Potentiel astronomique 
    3838    !............................................................................... 
    3939 
    4040    pot_astro(:,:)=0.e0 
     41    zramp = 1.e0 
    4142 
    4243    IF (lk_dynspg_ts) THEN 
    4344       zwt(:) = omega_tide(:)* ((kt-kt_tide)*rdt + kit*(rdt/REAL(nn_baro,wp))) 
     45       IF (ln_tide_ramp) THEN 
     46          zramp = MIN(MAX( ((kt-nit000)*rdt + kit*(rdt/REAL(nn_baro,wp)))/(rdttideramp*rday),0.),1.) 
     47       ENDIF 
    4448    ELSE 
    4549       zwt(:) = omega_tide(:)*(kt-kt_tide)*rdt 
     50       IF (ln_tide_ramp) THEN 
     51          zramp = MIN(MAX( ((kt-nit000)*rdt)/(rdttideramp*rday),0.),1.)  
     52       ENDIF   
    4653    ENDIF 
    4754 
     
    4956       do ji=1,jpi 
    5057          do jj=1,jpj 
    51              pot_astro(ji,jj)=pot_astro(ji,jj) + (amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))       
     58             pot_astro(ji,jj)=pot_astro(ji,jj) + zramp*(amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))       
    5259          enddo 
    5360       enddo 
  • branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r3367 r3490  
    7676   USE mod_ioclient 
    7777#endif 
     78   USE sbctide, ONLY: lk_tide 
     79 
    7880 
    7981   IMPLICIT NONE 
     
    318320                            CALL  istate_init   ! ocean initial state (Dynamics and tracers) 
    319321 
    320                             CALL tide_init( nit000 )    ! Initialisation of the tidal harmonics 
     322      IF( lk_tide       )   CALL tide_init( nit000 )    ! Initialisation of the tidal harmonics 
    321323 
    322324      IF( lk_bdy        )   CALL     bdy_init       ! Open boundaries initialisation 
  • branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step.F90

    r3367 r3490  
    9595      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    9696                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
    97       IF( kstp /= nit000 )   CALL tide_init ( kstp ) 
     97      IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp ) 
    9898      IF( lk_tide    )   CALL sbc_tide( kstp ) 
    9999      IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries 
     
    188188      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme 
    189189      IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends 
    190       IF( lk_bdy .AND.     & 
    191          & ln_tra_dmp_tot)   CALL bdy_tra_dmp( kstp )       ! bdy damping trends 
     190      IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends 
    192191                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection 
    193192      IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes 
     
    222221         & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment 
    223222      IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified) 
    224       IF( lk_bdy .AND.     & 
    225          & ln_dyn3d_dmp_tot)   CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends 
     223      IF( lk_bdy           )   CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends 
    226224                               CALL dyn_adv( kstp )         ! advection (vector or flux form) 
    227225                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis 
Note: See TracChangeset for help on using the changeset viewer.