Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 7351)
@@ -70,4 +70,10 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   ht_i   !: Now ice  thickness climatology
       REAL(wp), POINTER, DIMENSION(:,:) ::   ht_s   !: now snow thickness
+#endif
+#if defined key_top
+      CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply
+      REAL(wp)                            :: rn_fac  !: multiplicative scaling factor
+      REAL(wp), POINTER, DIMENSION(:,:)   :: trc     !: now field of the tracer
+      LOGICAL                             :: dmp     !: obc damping term
 #endif
    END TYPE OBC_DATA
Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 7351)
@@ -45,5 +45,4 @@
    INTEGER, ALLOCATABLE, DIMENSION(:)   ::   nb_bdy_fld        ! Number of fields to update for each boundary set.
    INTEGER                              ::   nb_bdy_fld_sum    ! Total number of fields to update for all boundary sets.
-
    LOGICAL,           DIMENSION(jp_bdy) ::   ln_full_vel_array ! =T => full velocities in 3D boundary conditions
                                                                ! =F => baroclinic velocities in 3D boundary conditions
@@ -58,5 +57,4 @@
 #endif
 
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -75,26 +73,24 @@
       !!                
       !!----------------------------------------------------------------------
-      !!
-      INTEGER, INTENT( in )           ::   kt    ! ocean time-step index 
-      INTEGER, INTENT( in ), OPTIONAL ::   jit   ! subcycle time-step index (for timesplitting option)
-      INTEGER, INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit
-                                                        ! is present then units = subcycle timesteps.
-                                                        ! time_offset = 0 => get data at "now" time level
-                                                        ! time_offset = -1 => get data at "before" time level
-                                                        ! time_offset = +1 => get data at "after" time level
-                                                        ! etc.
-      !!
-      INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices
+      INTEGER, INTENT(in)           ::   kt           ! ocean time-step index 
+      INTEGER, INTENT(in), OPTIONAL ::   jit          ! subcycle time-step index (for timesplitting option)
+      INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit
+      !                                               ! is present then units = subcycle timesteps.
+      !                                               ! time_offset = 0 => get data at "now" time level
+      !                                               ! time_offset = -1 => get data at "before" time level
+      !                                               ! time_offset = +1 => get data at "after" time level
+      !                                               ! etc.
+      !
+      INTEGER ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices
       INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
       INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
       TYPE(OBC_DATA), POINTER             ::   dta              ! short cut
-      !!
       !!---------------------------------------------------------------------------
-      !!
-      IF( nn_timing == 1 ) CALL timing_start('bdy_dta')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_dta')
+      !
       ! Initialise data arrays once for all from initial conditions where required
       !---------------------------------------------------------------------------
-      IF( kt .eq. nit000 .and. .not. PRESENT(jit) ) THEN
+      IF( kt == nit000 .AND. .NOT.PRESENT(jit) ) THEN
 
          ! Calculate depth-mean currents
@@ -102,10 +98,10 @@
          
          DO ib_bdy = 1, nb_bdy
-
+            !
             nblen => idx_bdy(ib_bdy)%nblen
             nblenrim => idx_bdy(ib_bdy)%nblenrim
             dta => dta_bdy(ib_bdy)
 
-            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 
+            IF( nn_dyn2d_dta(ib_bdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_ssh ) THEN 
@@ -135,5 +131,5 @@
             ENDIF
 
-            IF( nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
+            IF( nn_dyn3d_dta(ib_bdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_u3d ) THEN 
@@ -159,5 +155,5 @@
             ENDIF
 
-            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
+            IF( nn_tra_dta(ib_bdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_tem ) THEN
@@ -184,5 +180,5 @@
 
 #if defined key_lim2
-            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
+            IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_frld ) THEN
@@ -212,5 +208,5 @@
             ENDIF
 #elif defined key_lim3
-            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
+            IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_a_i ) THEN
@@ -246,9 +242,7 @@
             ENDIF
 #endif
-
-         ENDDO ! ib_bdy
-
-
-      ENDIF ! kt .eq. nit000
+         END DO ! ib_bdy
+         !
+      ENDIF ! kt == nit000
 
       ! update external data from files
@@ -258,5 +252,5 @@
       DO ib_bdy = 1, nb_bdy   
          dta => dta_bdy(ib_bdy)
-         IF( nn_dta(ib_bdy) .eq. 1 ) THEN ! skip this bit if no external data required
+         IF( nn_dta(ib_bdy) == 1 ) THEN ! skip this bit if no external data required
       
             IF( PRESENT(jit) ) THEN
@@ -264,11 +258,11 @@
                ! jit is optional argument for fld_read and bdytide_update
                IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
-                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                     IF( dta%ll_ssh ) dta%ssh(:) = 0.0
-                     IF( dta%ll_u2d ) dta%u2d(:) = 0.0
-                     IF( dta%ll_u3d ) dta%v2d(:) = 0.0
+                  IF( nn_dyn2d_dta(ib_bdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
+                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
+                     IF( dta%ll_u3d ) dta%v2d(:) = 0._wp
                   ENDIF
                   IF (cn_tra(ib_bdy) /= 'runoff') THEN
-                     IF( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 ) THEN
+                     IF( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 ) THEN
 
                         jend = jstart + dta%nread(2) - 1
@@ -278,9 +272,9 @@
                         ! If full velocities in boundary data then extract barotropic velocities from 3D fields
                         IF( ln_full_vel_array(ib_bdy) .AND.                                             &
-                          &    ( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR.  &
-                          &      nn_dyn3d_dta(ib_bdy) .EQ. 1 ) )THEN
+                          &    ( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 .OR.  &
+                          &      nn_dyn3d_dta(ib_bdy) == 1 ) )THEN
 
                            igrd = 2                      ! zonal velocity
-                           dta%u2d(:) = 0.0
+                           dta%u2d(:) = 0._wp
                            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                               ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -288,10 +282,10 @@
                               DO ik = 1, jpkm1
                                  dta%u2d(ib) = dta%u2d(ib) &
-                       &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
+                       &                          + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
                               END DO
-                              dta%u2d(ib) =  dta%u2d(ib) * hur(ii,ij)
+                              dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
                            END DO
                            igrd = 3                      ! meridional velocity
-                           dta%v2d(:) = 0.0
+                           dta%v2d(:) = 0._wp
                            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                               ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -299,7 +293,7 @@
                               DO ik = 1, jpkm1
                                  dta%v2d(ib) = dta%v2d(ib) &
-                       &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
+                       &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
                               END DO
-                              dta%v2d(ib) =  dta%v2d(ib) * hvr(ii,ij)
+                              dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
                            END DO
                         ENDIF                    
@@ -331,8 +325,8 @@
                   END DO
                ELSE
-                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                     IF( dta%ll_ssh ) dta%ssh(:) = 0.0
-                     IF( dta%ll_u2d ) dta%u2d(:) = 0.0
-                     IF( dta%ll_v2d ) dta%v2d(:) = 0.0
+                  IF( nn_dyn2d_dta(ib_bdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
+                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
+                     IF( dta%ll_v2d ) dta%v2d(:) = 0._wp
                   ENDIF
                   IF( dta%nread(1) .gt. 0 ) THEN ! update external data
@@ -343,8 +337,8 @@
                   ! If full velocities in boundary data then split into barotropic and baroclinic data
                   IF( ln_full_vel_array(ib_bdy) .and.                                             &
-                    & ( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR. &
-                    &   nn_dyn3d_dta(ib_bdy) .EQ. 1 ) ) THEN
+                    & ( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 .OR. &
+                    &   nn_dyn3d_dta(ib_bdy) == 1 ) ) THEN
                      igrd = 2                      ! zonal velocity
-                     dta%u2d(:) = 0.0
+                     dta%u2d(:) = 0._wp
                      DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                         ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -352,7 +346,7 @@
                         DO ik = 1, jpkm1
                            dta%u2d(ib) = dta%u2d(ib) &
-                 &                       + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
+                 &                       + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
                         END DO
-                        dta%u2d(ib) =  dta%u2d(ib) * hur(ii,ij)
+                        dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
                         DO ik = 1, jpkm1
                            dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib)
@@ -360,5 +354,5 @@
                      END DO
                      igrd = 3                      ! meridional velocity
-                     dta%v2d(:) = 0.0
+                     dta%v2d(:) = 0._wp
                      DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                         ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -366,7 +360,7 @@
                         DO ik = 1, jpkm1
                            dta%v2d(ib) = dta%v2d(ib) &
-                 &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
+                 &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
                         END DO
-                        dta%v2d(ib) =  dta%v2d(ib) * hvr(ii,ij)
+                        dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
                         DO ik = 1, jpkm1
                            dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib)
@@ -413,15 +407,15 @@
                   ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
                   dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij)
-               ENDDO
-            ENDIF
-         ENDDO
+               END DO
+            ENDIF
+         END DO
       ENDIF
-
+      !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dta')
-
-      END SUBROUTINE bdy_dta
-
-
-      SUBROUTINE bdy_dta_init
+      !
+   END SUBROUTINE bdy_dta
+
+
+   SUBROUTINE bdy_dta_init
       !!----------------------------------------------------------------------
       !!                   ***  SUBROUTINE bdy_dta_init  ***
@@ -433,8 +427,6 @@
       !!                
       !!----------------------------------------------------------------------
-      !!
-      INTEGER     ::  ib_bdy, jfld, jstart, jend, ierror  ! local indices
-      INTEGER      ::   ios                               ! Local integer output status for namelist read
-      !!
+      INTEGER ::   ib_bdy, jfld, jstart, jend, ierror, ios     ! Local integers
+      !
       CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
       CHARACTER(len=100), DIMENSION(nb_bdy)  ::   cn_dir_array  ! Root directory for location of data files
@@ -469,7 +461,7 @@
       NAMELIST/nambdy_dta/ ln_full_vel
       !!---------------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_init')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_dta_init')
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries'
@@ -486,5 +478,5 @@
 #endif
                               )
-         IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1
+         IF(nn_dta(ib_bdy) > 1) nn_dta(ib_bdy) = 1
       END DO
 
@@ -494,20 +486,20 @@
       nb_bdy_fld(:) = 0
       DO ib_bdy = 1, nb_bdy         
-         IF( cn_dyn2d(ib_bdy) /= 'none' .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN
+         IF( cn_dyn2d(ib_bdy) /= 'none' .and. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
          ENDIF
-         IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN
+         IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) == 1 ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
          ENDIF
-         IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN
+         IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) == 1  ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
          ENDIF
 #if ( defined key_lim2 || defined key_lim3 )
-         IF( cn_ice_lim(ib_bdy) /= 'none' .and. nn_ice_lim_dta(ib_bdy) .eq. 1  ) THEN
+         IF( cn_ice_lim(ib_bdy) /= 'none' .and. nn_ice_lim_dta(ib_bdy) == 1  ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
          ENDIF
 #endif               
          IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy)
-      ENDDO            
+      END DO            
 
       nb_bdy_fld_sum = SUM( nb_bdy_fld )
@@ -535,11 +527,11 @@
       jfld = 0 
       DO ib_bdy = 1, nb_bdy         
-         IF( nn_dta(ib_bdy) .eq. 1 ) THEN
+         IF( nn_dta(ib_bdy) == 1 ) THEN
             READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
-901         IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
+901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
 
             READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 )
-902         IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist', lwp )
-            IF(lwm) WRITE ( numond, nambdy_dta )
+902         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist', lwp )
+            IF(lwm) WRITE( numond, nambdy_dta )
 
             cn_dir_array(ib_bdy) = cn_dir
@@ -553,5 +545,5 @@
             ! Only read in necessary fields for this set.
             ! Important that barotropic variables come first.
-            IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN 
+            IF( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) THEN 
 
                IF( dta%ll_ssh ) THEN 
@@ -592,8 +584,8 @@
             ! read 3D velocities if baroclinic velocities require OR if
             ! barotropic velocities required and ln_full_vel set to .true.
-            IF( nn_dyn3d_dta(ib_bdy) .eq. 1 .or. &
-           &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN
-
-               IF( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN
+            IF( nn_dyn3d_dta(ib_bdy) == 1 .OR. &
+           &  ( ln_full_vel_array(ib_bdy) .AND. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) ) THEN
+
+               IF( dta%ll_u3d .OR. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in u3d field'
                   jfld = jfld + 1
@@ -606,5 +598,5 @@
                ENDIF
 
-               IF( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN
+               IF( dta%ll_v3d .OR. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in v3d field'
                   jfld = jfld + 1
@@ -620,5 +612,5 @@
 
             ! temperature and salinity
-            IF( nn_tra_dta(ib_bdy) .eq. 1 ) THEN
+            IF( nn_tra_dta(ib_bdy) == 1 ) THEN
 
                IF( dta%ll_tem ) THEN
@@ -646,5 +638,5 @@
 #if defined key_lim2
             ! sea ice
-            IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN
+            IF( nn_ice_lim_dta(ib_bdy) == 1 ) THEN
 
                IF( dta%ll_frld ) THEN
@@ -678,5 +670,5 @@
 #elif defined key_lim3
             ! sea ice
-            IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN
+            IF( nn_ice_lim_dta(ib_bdy) == 1 ) THEN
                ! Test for types of ice input (lim2 or lim3) 
                ! Build file name to find dimensions 
@@ -733,5 +725,5 @@
             ! Recalculate field counts
             !-------------------------
-            IF( ib_bdy .eq. 1 ) THEN 
+            IF( ib_bdy == 1 ) THEN 
                nb_bdy_fld_sum = 0
                nb_bdy_fld(ib_bdy) = jfld
@@ -744,5 +736,5 @@
             dta%nread(1) = nb_bdy_fld(ib_bdy)
 
-         ENDIF ! nn_dta .eq. 1
+         ENDIF ! nn_dta == 1
       ENDDO ! ib_bdy
 
@@ -785,5 +777,5 @@
          endif
 
-         IF ( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN
+         IF ( nn_dyn2d_dta(ib_bdy) == 0 .or. nn_dyn2d_dta(ib_bdy) == 2 ) THEN
             if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space'
             IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) )
@@ -791,5 +783,5 @@
             IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) )
          ENDIF
-         IF ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN
+         IF ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) THEN
             IF( dta%ll_ssh ) THEN
                if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow'
@@ -819,11 +811,11 @@
          ENDIF
 
-         IF ( nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN
+         IF ( nn_dyn3d_dta(ib_bdy) == 0 ) THEN
             if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space'
             IF( dta%ll_u3d ) ALLOCATE( dta_bdy(ib_bdy)%u3d(nblen(2),jpk) )
             IF( dta%ll_v3d ) ALLOCATE( dta_bdy(ib_bdy)%v3d(nblen(3),jpk) )
          ENDIF
-         IF ( nn_dyn3d_dta(ib_bdy) .eq. 1 .or. &
-           &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN
+         IF ( nn_dyn3d_dta(ib_bdy) == 1 .or. &
+           &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) ) THEN
             IF ( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN
                if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow'
@@ -838,5 +830,5 @@
          ENDIF
 
-         IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN
+         IF( nn_tra_dta(ib_bdy) == 0 ) THEN
             if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space'
             IF( dta%ll_tem ) ALLOCATE( dta_bdy(ib_bdy)%tem(nblen(1),jpk) )
@@ -857,5 +849,5 @@
 #if defined key_lim2
          IF (cn_ice_lim(ib_bdy) /= 'none') THEN
-            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN
+            IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN
                ALLOCATE( dta_bdy(ib_bdy)%frld(nblen(1)) )
                ALLOCATE( dta_bdy(ib_bdy)%hicif(nblen(1)) )
@@ -872,5 +864,5 @@
 #elif defined key_lim3
          IF (cn_ice_lim(ib_bdy) /= 'none') THEN
-            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN
+            IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN
                ALLOCATE( dta_bdy(ib_bdy)%a_i (nblen(1),jpl) )
                ALLOCATE( dta_bdy(ib_bdy)%ht_i(nblen(1),jpl) )
@@ -892,18 +884,18 @@
                   ALLOCATE( dta_bdy(ib_bdy)%ht_i(nblen(1),jpl) )
                   ALLOCATE( dta_bdy(ib_bdy)%ht_s(nblen(1),jpl) )
-                  dta_bdy(ib_bdy)%a_i (:,:) = 0.0
-                  dta_bdy(ib_bdy)%ht_i(:,:) = 0.0
-                  dta_bdy(ib_bdy)%ht_s(:,:) = 0.0
-               ENDIF
-
-            ENDIF
-         ENDIF
-#endif
-
-      ENDDO ! ib_bdy 
-
+                  dta_bdy(ib_bdy)%a_i (:,:) = 0._wp
+                  dta_bdy(ib_bdy)%ht_i(:,:) = 0._wp
+                  dta_bdy(ib_bdy)%ht_s(:,:) = 0._wp
+               ENDIF
+
+            ENDIF
+         ENDIF
+#endif
+         !
+      END DO ! ib_bdy 
+      !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_init')
-
-      END SUBROUTINE bdy_dta_init
+      !
+   END SUBROUTINE bdy_dta_init
 
 #else
Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 7351)
@@ -36,5 +36,4 @@
    PUBLIC   bdy_dyn    ! routine called in dyn_nxt
 
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -51,27 +50,26 @@
       !!
       !!----------------------------------------------------------------------
-      !!
-      INTEGER, INTENT( in )           :: kt               ! Main time step counter
-      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities
-      !!
-      INTEGER               :: jk,ii,ij,ib_bdy,ib,igrd     ! Loop counter
-      LOGICAL               :: ll_dyn2d, ll_dyn3d, ll_orlanski
-      !!
+      INTEGER, INTENT(in)           ::   kt           ! Main time step counter
+      LOGICAL, INTENT(in), OPTIONAL ::   dyn3d_only   ! T => only update baroclinic velocities
+      !
+      INTEGER ::   jk, ii, ij, ib_bdy, ib, igrd     ! Loop counter
+      LOGICAL ::   ll_dyn2d, ll_dyn3d, ll_orlanski
       REAL(wp), POINTER, DIMENSION(:,:) :: pua2d, pva2d     ! after barotropic velocities
-
-      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn')
-
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_dyn')
+      !
       ll_dyn2d = .true.
       ll_dyn3d = .true.
-
+      !
       IF( PRESENT(dyn3d_only) ) THEN
-         IF( dyn3d_only ) ll_dyn2d = .false.
+         IF( dyn3d_only )   ll_dyn2d = .false.
       ENDIF
-
+      !
       ll_orlanski = .false.
       DO ib_bdy = 1, nb_bdy
-         IF ( cn_dyn2d(ib_bdy) == 'orlanski' .or. cn_dyn2d(ib_bdy) == 'orlanski_npo' &
-     &   .or. cn_dyn3d(ib_bdy) == 'orlanski' .or. cn_dyn3d(ib_bdy) == 'orlanski_npo') ll_orlanski = .true.
-      ENDDO
+         IF ( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' &
+     &   .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo')   ll_orlanski = .true.
+      END DO
 
       !-------------------------------------------------------
@@ -79,5 +77,5 @@
       !-------------------------------------------------------
 
-      CALL wrk_alloc(jpi,jpj,pua2d,pva2d) 
+      CALL wrk_alloc( jpi,jpj,   pua2d, pva2d ) 
 
       !-------------------------------------------------------
@@ -85,29 +83,26 @@
       !-------------------------------------------------------
 
-      ! "After" velocities: 
+      !                          ! "After" velocities: 
+      pua2d(:,:) = 0._wp
+      pva2d(:,:) = 0._wp     
+      DO jk = 1, jpkm1
+         pua2d(:,:) = pua2d(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+         pva2d(:,:) = pva2d(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+      END DO
+      pua2d(:,:) = pua2d(:,:) * r1_hu_a(:,:)
+      pva2d(:,:) = pva2d(:,:) * r1_hv_a(:,:)
 
-      pua2d(:,:) = 0.e0
-      pva2d(:,:) = 0.e0      
-      DO jk = 1, jpkm1
-         pua2d(:,:) = pua2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-         pva2d(:,:) = pva2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+      DO jk = 1 , jpkm1
+         ua(:,:,jk) = ( ua(:,:,jk) - pua2d(:,:) ) * umask(:,:,jk)
+         va(:,:,jk) = ( va(:,:,jk) - pva2d(:,:) ) * vmask(:,:,jk)
       END DO
 
-      pua2d(:,:) = pua2d(:,:) * hur_a(:,:)
-      pva2d(:,:) = pva2d(:,:) * hvr_a(:,:)
 
-      DO jk = 1 , jpkm1
-         ua(:,:,jk) = (ua(:,:,jk) - pua2d(:,:)) * umask(:,:,jk)
-         va(:,:,jk) = (va(:,:,jk) - pva2d(:,:)) * vmask(:,:,jk)
-      END DO
-
-      ! "Before" velocities (required for Orlanski condition): 
-
-      IF ( ll_orlanski ) THEN          
+      IF( ll_orlanski ) THEN     ! "Before" velocities (Orlanski condition only) 
          DO jk = 1 , jpkm1
-            ub(:,:,jk) = (ub(:,:,jk) - ub_b(:,:)) * umask(:,:,jk)
-            vb(:,:,jk) = (vb(:,:,jk) - vb_b(:,:)) * vmask(:,:,jk)
+            ub(:,:,jk) = ( ub(:,:,jk) - ub_b(:,:) ) * umask(:,:,jk)
+            vb(:,:,jk) = ( vb(:,:,jk) - vb_b(:,:) ) * vmask(:,:,jk)
          END DO
-      END IF
+      ENDIF
 
       !-------------------------------------------------------
@@ -116,17 +111,17 @@
       !-------------------------------------------------------
 
-      IF( ll_dyn2d ) CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, hur_a(:,:), hvr_a(:,:), ssha )
+      IF( ll_dyn2d )   CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, r1_hu_a(:,:), r1_hv_a(:,:), ssha )
 
-      IF( ll_dyn3d ) CALL bdy_dyn3d( kt )
+      IF( ll_dyn3d )   CALL bdy_dyn3d( kt )
 
       !-------------------------------------------------------
       ! Recombine velocities
       !-------------------------------------------------------
-
+      !
       DO jk = 1 , jpkm1
          ua(:,:,jk) = ( ua(:,:,jk) + pua2d(:,:) ) * umask(:,:,jk)
          va(:,:,jk) = ( va(:,:,jk) + pva2d(:,:) ) * vmask(:,:,jk)
       END DO
-
+      !
       IF ( ll_orlanski ) THEN
          DO jk = 1 , jpkm1
@@ -135,9 +130,9 @@
          END DO
       END IF
-
-      CALL wrk_dealloc(jpi,jpj,pua2d,pva2d) 
-
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn')
-
+      !
+      CALL wrk_dealloc( jpi,jpj,  pua2d, pva2d ) 
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_dyn')
+      !
    END SUBROUTINE bdy_dyn
 
Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 7351)
@@ -29,6 +29,4 @@
    PUBLIC   bdy_dyn3d_dmp ! routine called by step
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -45,30 +43,24 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) :: kt     ! Main time step counter
-      !!
-      INTEGER               :: ib_bdy ! loop index
-      !!
-
+      INTEGER, INTENT(in) ::   kt   ! Main time step counter
+      !
+      INTEGER ::   ib_bdy   ! loop index
+      !!----------------------------------------------------------------------
+      !
       DO ib_bdy=1, nb_bdy
-
+         !
          SELECT CASE( cn_dyn3d(ib_bdy) )
-         CASE('none')
-            CYCLE
-         CASE('frs')
-            CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE('specified')
-            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE('zero')
-            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE('orlanski')
-            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
-         CASE('orlanski_npo')
-            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
-         CASE DEFAULT
-            CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
+         CASE('none')        ;   CYCLE
+         CASE('frs' )        ;   CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+         CASE('specified')   ;   CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+         CASE('zero')        ;   CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+         CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
+         CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
+         CASE DEFAULT        ;   CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
          END SELECT
-      ENDDO
-
+      END DO
+      !
    END SUBROUTINE bdy_dyn3d
+
 
    SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy )
@@ -80,9 +72,9 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER                     ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
-      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
-      !!
+      INTEGER        , INTENT(in) ::   kt      ! time step index
+      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
+      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
+      INTEGER        , INTENT(in) ::   ib_bdy  ! BDY set index
+      !
       INTEGER  ::   jb, jk         ! dummy loop indices
       INTEGER  ::   ii, ij, igrd   ! local integers
@@ -112,10 +104,11 @@
       CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
       !
-      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
-
+      IF( kt == nit000 )   CLOSE( unit = 102 )
+      !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe')
-
+      !
    END SUBROUTINE bdy_dyn3d_spe
 
+
    SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
       !!----------------------------------------------------------------------
@@ -125,11 +118,11 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER                     ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      INTEGER        , INTENT(in) ::   kt      ! time step index
+      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
+      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
-      !!
+      !
       INTEGER  ::   ib, ik         ! dummy loop indices
-      INTEGER  ::   ii, ij, igrd, zcoef   ! local integers
+      INTEGER  ::   ii, ij, igrd   ! local integers
       REAL(wp) ::   zwgt           ! boundary weight
       !!----------------------------------------------------------------------
@@ -157,9 +150,10 @@
       CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
       !
-      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
-
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro')
-
+      IF( kt == nit000 )   CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_dyn3d_zro')
+      !
    END SUBROUTINE bdy_dyn3d_zro
+
 
    SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy )
@@ -174,9 +168,9 @@
       !!               topography. Tellus, 365-382.
       !!----------------------------------------------------------------------
-      INTEGER                     ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      INTEGER        , INTENT(in) ::   kt      ! time step index
+      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
+      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
-      !!
+      !
       INTEGER  ::   jb, jk         ! dummy loop indices
       INTEGER  ::   ii, ij, igrd   ! local integers
@@ -208,9 +202,10 @@
       CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
       !
-      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
-
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_frs')
-
+      IF( kt == nit000 )   CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_dyn3d_frs')
+      !
    END SUBROUTINE bdy_dyn3d_frs
+
 
    SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, ll_npo )
@@ -259,16 +254,14 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER                     ::   kt
-      !!
+      INTEGER, INTENT(in) ::   kt   ! time step index
+      !
       INTEGER  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ib_bdy         ! loop index
       INTEGER  ::   ii, ij, igrd   ! local integers
       REAL(wp) ::   zwgt           ! boundary weight
-      INTEGER  ::  ib_bdy          ! loop index
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp')
-      !
-      !-------------------------------------------------------
-
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_dyn3d_dmp')
+      !
       DO ib_bdy=1, nb_bdy
          IF ( ln_dyn3d_dmp(ib_bdy) .and. cn_dyn3d(ib_bdy) /= 'none' ) THEN
@@ -295,10 +288,10 @@
             END DO
          ENDIF
-      ENDDO
+      END DO
       !
       CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
       !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp')
-
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_dyn3d_dmp')
+      !
    END SUBROUTINE bdy_dyn3d_dmp
 
@@ -311,9 +304,7 @@
       WRITE(*,*) 'bdy_dyn3d: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_dyn3d
-
    SUBROUTINE bdy_dyn3d_dmp( kt )      ! Empty routine
       WRITE(*,*) 'bdy_dyn3d_dmp: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_dyn3d_dmp
-
 #endif
 
Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 7351)
@@ -12,6 +12,5 @@
    !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
    !!            3.4  !  2012     (J. Chanut) straight open boundary case update
-   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Updates for the 
-   !!                             optimization of BDY communications
+   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) optimization of BDY communications
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -19,17 +18,18 @@
    !!   'key_bdy'                     Unstructured Open Boundary Conditions
    !!----------------------------------------------------------------------
-   !!   bdy_init       : Initialization of unstructured open boundaries
+   !!   bdy_init      : Initialization of unstructured open boundaries
    !!----------------------------------------------------------------------
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain
-   USE bdy_oce         ! unstructured open boundary conditions
-   USE in_out_manager  ! I/O units
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         ! for mpp_sum  
-   USE iom             ! I/O
-   USE sbctide, ONLY: lk_tide ! Tidal forcing or not
-   USE phycst, ONLY: rday
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain
+   USE bdy_oce        ! unstructured open boundary conditions
+   USE sbctide  , ONLY: lk_tide ! Tidal forcing or not
+   USE phycst   , ONLY: rday
+   !
+   USE in_out_manager ! I/O units
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! for mpp_sum  
+   USE iom            ! I/O
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -38,15 +38,15 @@
    PUBLIC   bdy_init   ! routine called in nemo_init
 
-   INTEGER, PARAMETER          :: jp_nseg = 100
-   INTEGER, PARAMETER          :: nrimmax = 20 ! maximum rimwidth in structured
+   INTEGER, PARAMETER ::   jp_nseg = 100   ! 
+   INTEGER, PARAMETER ::   nrimmax =  20   ! maximum rimwidth in structured
                                                ! open boundary data files
    ! Straight open boundary segment parameters:
-   INTEGER  :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
-   INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge
-   INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw
-   INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn
-   INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs
+   INTEGER  ::   nbdysege, nbdysegw, nbdysegn, nbdysegs 
+   INTEGER, DIMENSION(jp_nseg) ::   jpieob, jpjedt, jpjeft, npckge   !
+   INTEGER, DIMENSION(jp_nseg) ::   jpiwob, jpjwdt, jpjwft, npckgw   !
+   INTEGER, DIMENSION(jp_nseg) ::   jpjnob, jpindt, jpinft, npckgn   !
+   INTEGER, DIMENSION(jp_nseg) ::   jpjsob, jpisdt, jpisft, npckgs   !
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -66,9 +66,4 @@
       !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
       !!----------------------------------------------------------------------      
-      ! namelist variables
-      !-------------------
-      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
-      CHARACTER(LEN=1)   ::   ctypebdy
-      INTEGER :: nbdyind, nbdybeg, nbdyend
 
       ! local variables
@@ -81,5 +76,5 @@
       INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
       INTEGER  ::   i_offset, j_offset                     !   -       -
-      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
+      INTEGER , POINTER  ::  nbi, nbj, nbr                 ! short cuts
       REAL(wp), POINTER  ::  flagu, flagv                  !    -   -
       REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields
@@ -94,5 +89,8 @@
       INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates
       REAL(wp), POINTER, DIMENSION(:,:)       ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat)
-
+      !!
+      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile     ! Namelist variables
+      CHARACTER(LEN=1)                     ::   ctypebdy   !     -        - 
+      INTEGER                              ::   nbdyind, nbdybeg, nbdyend
       !!
       NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,                 &
@@ -103,32 +101,30 @@
          &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 &
          &             ln_vol, nn_volctl, nn_rimwidth
-      !!
+         !
       NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('bdy_init')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_init')
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
       IF(lwp) WRITE(numout,*) '~~~~~~~~'
       !
-
       IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
          &                               ' and general open boundary condition are not compatible' )
 
-      cgrid= (/'t','u','v'/)
+      cgrid = (/'t','u','v'/)
       
       ! ------------------------
       ! Read namelist parameters
       ! ------------------------
-
       REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries  
       READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
+      !
       REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
       READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, nambdy )
 
@@ -137,10 +133,10 @@
       ! -----------------------------------------
       !                                   ! control prints
-      IF(lwp) WRITE(numout,*) '         nambdy'
-
-      IF( nb_bdy .eq. 0 ) THEN 
+      IF(lwp) WRITE(numout,*) '   nambdy'
+
+      IF( nb_bdy == 0 ) THEN 
         IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
       ELSE
-        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
+        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
       ENDIF
 
@@ -158,25 +154,25 @@
         IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
         SELECT CASE( cn_dyn2d(ib_bdy) )                  
-          CASE('none')         
+          CASE( 'none' )         
              IF(lwp) WRITE(numout,*) '      no open boundary condition'        
              dta_bdy(ib_bdy)%ll_ssh = .false.
              dta_bdy(ib_bdy)%ll_u2d = .false.
              dta_bdy(ib_bdy)%ll_v2d = .false.
-          CASE('frs')          
+          CASE( 'frs' )          
              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
              dta_bdy(ib_bdy)%ll_ssh = .false.
              dta_bdy(ib_bdy)%ll_u2d = .true.
              dta_bdy(ib_bdy)%ll_v2d = .true.
-          CASE('flather')      
+          CASE( 'flather' )      
              IF(lwp) WRITE(numout,*) '      Flather radiation condition'
              dta_bdy(ib_bdy)%ll_ssh = .true.
              dta_bdy(ib_bdy)%ll_u2d = .true.
              dta_bdy(ib_bdy)%ll_v2d = .true.
-          CASE('orlanski')     
+          CASE( 'orlanski' )     
              IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
              dta_bdy(ib_bdy)%ll_ssh = .false.
              dta_bdy(ib_bdy)%ll_u2d = .true.
              dta_bdy(ib_bdy)%ll_v2d = .true.
-          CASE('orlanski_npo') 
+          CASE( 'orlanski_npo' ) 
              IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
              dta_bdy(ib_bdy)%ll_ssh = .false.
@@ -392,8 +388,4 @@
       REWIND( numnam_cfg )     
 
-      !!----------------------------------------------------------------------
-
-              
-               
       nblendta(:,:) = 0
       nbdysege = 0
@@ -492,10 +484,10 @@
                nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
                jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))
-            ENDDO
+            END DO
             CALL iom_close( inum )
-
+            !
          ENDIF 
-
-      ENDDO ! ib_bdy
+         !
+      END DO ! ib_bdy
 
       IF (nb_bdy>0) THEN
@@ -514,5 +506,5 @@
       ! Now look for crossings in user (namelist) defined open boundary segments:
       !--------------------------------------------------------------------------
-      IF ( icount>0 ) CALL bdy_ctl_seg
+      IF( icount>0 )   CALL bdy_ctl_seg
 
       ! Calculate global boundary index arrays or read in from file
@@ -520,7 +512,7 @@
       ! 1. Read global index arrays from boundary coordinates file.
       DO ib_bdy = 1, nb_bdy
-
+         !
          IF( ln_coords_file(ib_bdy) ) THEN
-
+            !
             CALL iom_open( cn_coords_file(ib_bdy), inum )
             DO igrd = 1, jpbgrd
@@ -537,5 +529,5 @@
                   nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
                END DO
-
+               !
                ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
                IF(lwp) WRITE(numout,*)
@@ -546,8 +538,8 @@
             END DO
             CALL iom_close( inum )
-
+            !
          ENDIF 
-
-      ENDDO      
+         !
+      END DO      
     
       ! 2. Now fill indices corresponding to straight open boundary arrays:
@@ -792,5 +784,5 @@
 
       ! Work out dimensions of boundary data on each neighbour process
-      IF(nbondi .eq. 0) THEN
+      IF(nbondi == 0) THEN
          iw_b(1) = jpizoom + nimppt(nowe+1)
          ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3
@@ -802,10 +794,10 @@
          is_b(2) = jpjzoom + njmppt(noea+1)
          in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3
-      ELSEIF(nbondi .eq. 1) THEN
+      ELSEIF(nbondi == 1) THEN
          iw_b(1) = jpizoom + nimppt(nowe+1)
          ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3
          is_b(1) = jpjzoom + njmppt(nowe+1)
          in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3
-      ELSEIF(nbondi .eq. -1) THEN
+      ELSEIF(nbondi == -1) THEN
          iw_b(2) = jpizoom + nimppt(noea+1)
          ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3
@@ -814,5 +806,5 @@
       ENDIF
 
-      IF(nbondj .eq. 0) THEN
+      IF(nbondj == 0) THEN
          iw_b(3) = jpizoom + nimppt(noso+1)
          ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3
@@ -824,10 +816,10 @@
          is_b(4) = jpjzoom + njmppt(nono+1)
          in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3
-      ELSEIF(nbondj .eq. 1) THEN
+      ELSEIF(nbondj == 1) THEN
          iw_b(3) = jpizoom + nimppt(noso+1)
          ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3
          is_b(3) = jpjzoom + njmppt(noso+1)
          in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3
-      ELSEIF(nbondj .eq. -1) THEN
+      ELSEIF(nbondj == -1) THEN
          iw_b(4) = jpizoom + nimppt(nono+1)
          ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3
@@ -867,25 +859,25 @@
          ! Allocate index arrays for this boundary set
          !--------------------------------------------
-         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
-         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) )
+         ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbi   (ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbj   (ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbr   (ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbd   (ilen1,jpbgrd) )
          ALLOCATE( idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbw   (ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
 
          ! Dispatch mapping indices and discrete distances on each processor
          ! -----------------------------------------------------------------
 
-         com_east = 0
-         com_west = 0
+         com_east  = 0
+         com_west  = 0
          com_south = 0
          com_north = 0
 
-         com_east_b = 0
-         com_west_b = 0
+         com_east_b  = 0
+         com_west_b  = 0
          com_south_b = 0
          com_north_b = 0
@@ -912,12 +904,12 @@
                      ii = idx_bdy(ib_bdy)%nbi(icount,igrd)
                      ij = idx_bdy(ib_bdy)%nbj(icount,igrd)
-                     if((com_east .ne. 1) .and. (ii .eq. (nlci-1)) .and. (nbondi .le. 0)) then
+                     if((com_east .ne. 1) .and. (ii == (nlci-1)) .and. (nbondi .le. 0)) then
                         com_east = 1
-                     elseif((com_west .ne. 1) .and. (ii .eq. 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
+                     elseif((com_west .ne. 1) .and. (ii == 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
                         com_west = 1
                      endif 
-                     if((com_south .ne. 1) .and. (ij .eq. 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
+                     if((com_south .ne. 1) .and. (ij == 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
                         com_south = 1
-                     elseif((com_north .ne. 1) .and. (ij .eq. (nlcj-1)) .and. (nbondj .le. 0)) then
+                     elseif((com_north .ne. 1) .and. (ij == (nlcj-1)) .and. (nbondj .le. 0)) then
                         com_north = 1
                      endif 
@@ -926,14 +918,14 @@
                   ENDIF
                   ! check if point has to be received from a neighbour
-                  IF(nbondi .eq. 0) THEN
+                  IF(nbondi == 0) THEN
                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   &
                        & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   &
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
-                       if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
+                       if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
-                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
+                          if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
                             com_south = 1
-                          elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
+                          elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
                             com_north = 1
                           endif
@@ -945,9 +937,9 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
-                       if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
+                       if((com_east_b .ne. 1) .and. (ii == 2)) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
-                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
+                          if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
                             com_south = 1
-                          elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
+                          elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
                             com_north = 1
                           endif
@@ -955,14 +947,14 @@
                        endif 
                      ENDIF
-                  ELSEIF(nbondi .eq. 1) THEN
+                  ELSEIF(nbondi == 1) THEN
                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   &
                        & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   &
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
-                       if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
+                       if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
-                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
+                          if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
                             com_south = 1
-                          elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
+                          elseif((ij == nlcjt(nowe+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
                             com_north = 1
                           endif
@@ -970,14 +962,14 @@
                        endif 
                      ENDIF
-                  ELSEIF(nbondi .eq. -1) THEN
+                  ELSEIF(nbondi == -1) THEN
                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   &
                        & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   &
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
-                       if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
+                       if((com_east_b .ne. 1) .and. (ii == 2)) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
-                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
+                          if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
                             com_south = 1
-                          elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
+                          elseif((ij == nlcjt(noea+1)-1) .and. (nbondj == 0 .or. nbondj == -1)) then
                             com_north = 1
                           endif
@@ -986,5 +978,5 @@
                      ENDIF
                   ENDIF
-                  IF(nbondj .eq. 0) THEN
+                  IF(nbondj == 0) THEN
                      IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  &
                        & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
@@ -1001,5 +993,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
-                       if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
+                       if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then
                           com_south_b = 1
                        endif 
@@ -1009,9 +1001,9 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
-                       if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
+                       if((com_north_b .ne. 1) .and. (ij == 2)) then
                           com_north_b = 1
                        endif 
                      ENDIF
-                  ELSEIF(nbondj .eq. 1) THEN
+                  ELSEIF(nbondj == 1) THEN
                      IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. &
                        & nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
@@ -1023,9 +1015,9 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
-                       if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
+                       if((com_south_b .ne. 1) .and. (ij == (nlcjt(noso+1)-1))) then
                           com_south_b = 1
                        endif 
                      ENDIF
-                  ELSEIF(nbondj .eq. -1) THEN
+                  ELSEIF(nbondj == -1) THEN
                      IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  &
                        & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
@@ -1037,5 +1029,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
-                       if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
+                       if((com_north_b .ne. 1) .and. (ij == 2)) then
                           com_north_b = 1
                        endif 
@@ -1046,38 +1038,22 @@
          ENDDO 
 
-         ! definition of the i- and j- direction local boundaries arrays
-         ! used for sending the boudaries
-         IF((com_east .eq. 1) .and. (com_west .eq. 1)) THEN
-            nbondi_bdy(ib_bdy) = 0
-         ELSEIF ((com_east .eq. 1) .and. (com_west .eq. 0)) THEN
-            nbondi_bdy(ib_bdy) = -1
-         ELSEIF ((com_east .eq. 0) .and. (com_west .eq. 1)) THEN
-            nbondi_bdy(ib_bdy) = 1
+         ! definition of the i- and j- direction local boundaries arrays used for sending the boundaries
+         IF(     (com_east  == 1) .and. (com_west  == 1) ) THEN   ;   nbondi_bdy(ib_bdy) =  0
+         ELSEIF( (com_east  == 1) .and. (com_west  == 0) ) THEN   ;   nbondi_bdy(ib_bdy) = -1
+         ELSEIF( (com_east  == 0) .and. (com_west  == 1) ) THEN   ;   nbondi_bdy(ib_bdy) =  1
          ENDIF
-
-         IF((com_north .eq. 1) .and. (com_south .eq. 1)) THEN
-            nbondj_bdy(ib_bdy) = 0
-         ELSEIF ((com_north .eq. 1) .and. (com_south .eq. 0)) THEN
-            nbondj_bdy(ib_bdy) = -1
-         ELSEIF ((com_north .eq. 0) .and. (com_south .eq. 1)) THEN
-            nbondj_bdy(ib_bdy) = 1
+         IF(     (com_north == 1) .and. (com_south == 1) ) THEN   ;   nbondj_bdy(ib_bdy) =  0
+         ELSEIF( (com_north == 1) .and. (com_south == 0) ) THEN   ;   nbondj_bdy(ib_bdy) = -1
+         ELSEIF( (com_north == 0) .and. (com_south == 1) ) THEN   ;   nbondj_bdy(ib_bdy) =  1
          ENDIF
 
-         ! definition of the i- and j- direction local boundaries arrays
-         ! used for receiving the boudaries
-         IF((com_east_b .eq. 1) .and. (com_west_b .eq. 1)) THEN
-            nbondi_bdy_b(ib_bdy) = 0
-         ELSEIF ((com_east_b .eq. 1) .and. (com_west_b .eq. 0)) THEN
-            nbondi_bdy_b(ib_bdy) = -1
-         ELSEIF ((com_east_b .eq. 0) .and. (com_west_b .eq. 1)) THEN
-            nbondi_bdy_b(ib_bdy) = 1
+         ! definition of the i- and j- direction local boundaries arrays used for receiving the boundaries
+         IF(     (com_east_b  == 1) .and. (com_west_b  == 1) ) THEN   ;   nbondi_bdy_b(ib_bdy) =  0
+         ELSEIF( (com_east_b  == 1) .and. (com_west_b  == 0) ) THEN   ;   nbondi_bdy_b(ib_bdy) = -1
+         ELSEIF( (com_east_b  == 0) .and. (com_west_b  == 1) ) THEN   ;   nbondi_bdy_b(ib_bdy) =  1
          ENDIF
-
-         IF((com_north_b .eq. 1) .and. (com_south_b .eq. 1)) THEN
-            nbondj_bdy_b(ib_bdy) = 0
-         ELSEIF ((com_north_b .eq. 1) .and. (com_south_b .eq. 0)) THEN
-            nbondj_bdy_b(ib_bdy) = -1
-         ELSEIF ((com_north_b .eq. 0) .and. (com_south_b .eq. 1)) THEN
-            nbondj_bdy_b(ib_bdy) = 1
+         IF(     (com_north_b == 1) .and. (com_south_b == 1) ) THEN   ;   nbondj_bdy_b(ib_bdy) =  0
+         ELSEIF( (com_north_b == 1) .and. (com_south_b == 0) ) THEN   ;   nbondj_bdy_b(ib_bdy) = -1
+         ELSEIF( (com_north_b == 0) .and. (com_south_b == 1) ) THEN   ;   nbondj_bdy_b(ib_bdy) =  1
          ENDIF
 
@@ -1087,7 +1063,7 @@
             DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
-               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
-!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic
-!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear
+               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( nbr - 1 ) *0.5 )      ! tanh formulation
+!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic
+!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy))       ! linear
             END DO
          END DO 
@@ -1099,7 +1075,7 @@
                nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
                idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
-               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic
+               & *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
                idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 
-               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic
+               & *(REAL(nn_rimwidth(ib_bdy)+1-nbr)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
             END DO
          END DO 
@@ -1122,10 +1098,10 @@
 
          ! Derive mask on U and V grid from mask on T grid
-         bdyumask(:,:) = 0.e0
-         bdyvmask(:,:) = 0.e0
+         bdyumask(:,:) = 0._wp
+         bdyvmask(:,:) = 0._wp
          DO ij=1, jpjm1
             DO ii=1, jpim1
-               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
-               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)  
+               bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij )
+               bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1)  
             END DO
          END DO
@@ -1141,10 +1117,6 @@
                   umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
                   vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
-                  bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
                END DO      
             END DO
-         END DO
-
-         DO ik = 1, jpkm1
             DO ij = 2, jpjm1
                DO ii = 2, jpim1
@@ -1154,19 +1126,17 @@
             END DO
          END DO
-
          tmask_i (:,:) = ssmask(:,:) * tmask_i(:,:)
-
+         !
       ENDIF ! ln_mask_file=.TRUE.
       
       bdytmask(:,:) = ssmask(:,:)
-      IF( .not. ln_mask_file ) THEN
-         ! If .not. ln_mask_file then we need to derive mask on U and V grid 
-         ! from mask on T grid here.
-         bdyumask(:,:) = 0.e0
-         bdyvmask(:,:) = 0.e0
-         DO ij=1, jpjm1
-            DO ii=1, jpim1
-               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
-               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)  
+      IF( .NOT.ln_mask_file ) THEN
+         ! If .not. ln_mask_file then we need to derive mask on U and V grid from mask on T grid here.
+         bdyumask(:,:) = 0._wp
+         bdyvmask(:,:) = 0._wp
+         DO ij = 1, jpjm1
+            DO ii = 1, jpim1
+               bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij )
+               bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1)  
             END DO
          END DO
@@ -1174,30 +1144,24 @@
       ENDIF
 
-      ! bdy masks and bmask are now set to zero on boundary points:
-      igrd = 1       ! In the free surface case, bmask is at T-points
-      DO ib_bdy = 1, nb_bdy
-        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
-          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
-        ENDDO
-      ENDDO
+      ! bdy masks are now set to zero on boundary points:
       !
       igrd = 1
       DO ib_bdy = 1, nb_bdy
         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)      
-          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
-        ENDDO
-      ENDDO
+          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
+        END DO
+      END DO
       !
       igrd = 2
       DO ib_bdy = 1, nb_bdy
         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
-        ENDDO
-      ENDDO
+          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
+        END DO
+      END DO
       !
       igrd = 3
       DO ib_bdy = 1, nb_bdy
         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
+          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
         ENDDO
       ENDDO
@@ -1205,5 +1169,5 @@
       ! For the flagu/flagv calculation below we require a version of fmask without
       ! the land boundary condition (shlat) included:
-      CALL wrk_alloc(jpi,jpj,zfmask) 
+      CALL wrk_alloc(jpi,jpj,  zfmask ) 
       DO ij = 2, jpjm1
          DO ii = 2, jpim1
@@ -1220,6 +1184,6 @@
       DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
 
-         idx_bdy(ib_bdy)%flagu(:,:) = 0.e0
-         idx_bdy(ib_bdy)%flagv(:,:) = 0.e0
+         idx_bdy(ib_bdy)%flagu(:,:) = 0._wp
+         idx_bdy(ib_bdy)%flagv(:,:) = 0._wp
          icount = 0 
 
@@ -1231,13 +1195,7 @@
          DO igrd = 1,jpbgrd 
             SELECT CASE( igrd )
-               CASE( 1 )
-                  pmask => umask(:,:,1)
-                  i_offset = 0
-               CASE( 2 ) 
-                  pmask => bdytmask
-                  i_offset = 1
-               CASE( 3 ) 
-                  pmask => zfmask(:,:)
-                  i_offset = 0
+               CASE( 1 )   ;   pmask => umask   (:,:,1)   ;   i_offset = 0
+               CASE( 2 )   ;   pmask => bdytmask(:,:)     ;   i_offset = 1
+               CASE( 3 )   ;   pmask => zfmask  (:,:)     ;   i_offset = 0
             END SELECT 
             icount = 0
@@ -1270,15 +1228,9 @@
          ! flagv =  1 : v is normal to the boundary and is direction is inward
 
-         DO igrd = 1,jpbgrd 
+         DO igrd = 1, jpbgrd 
             SELECT CASE( igrd )
-               CASE( 1 )
-                  pmask => vmask(:,:,1)
-                  j_offset = 0
-               CASE( 2 )
-                  pmask => zfmask(:,:)
-                  j_offset = 0
-               CASE( 3 )
-                  pmask => bdytmask
-                  j_offset = 1
+               CASE( 1 )   ;   pmask => vmask (:,:,1)   ;   j_offset = 0
+               CASE( 2 )   ;   pmask => zfmask(:,:)     ;   j_offset = 0
+               CASE( 3 )   ;   pmask => bdytmask        ;   j_offset = 1
             END SELECT 
             icount = 0
@@ -1286,6 +1238,6 @@
                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
-               znfl = pmask(nbi,nbj+j_offset-1  )
-               zsfl = pmask(nbi,nbj+j_offset)
+               znfl = pmask(nbi,nbj+j_offset-1)
+               zsfl = pmask(nbi,nbj+j_offset  )
                ! This error check only works if you are using the bdyXmask arrays
                IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN
@@ -1305,11 +1257,11 @@
             ENDIF 
          END DO
-
+         !
       END DO
 
       ! Compute total lateral surface for volume correction:
       ! ----------------------------------------------------
-      ! JC: this must be done at each time step with key_vvl
-      bdysurftot = 0.e0 
+      ! JC: this must be done at each time step with non-linear free surface
+      bdysurftot = 0._wp 
       IF( ln_vol ) THEN  
          igrd = 2      ! Lateral surface at U-points
@@ -1319,10 +1271,10 @@
                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
                flagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
-               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
+               bdysurftot = bdysurftot + hu_n   (nbi  , nbj)                           &
                   &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
                   &                    * tmask_i(nbi  , nbj)                           &
                   &                    * tmask_i(nbi+1, nbj)                   
-            ENDDO
-         ENDDO
+            END DO
+         END DO
 
          igrd=3 ! Add lateral surface at V-points
@@ -1332,10 +1284,10 @@
                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
                flagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
-               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
+               bdysurftot = bdysurftot + hv_n   (nbi, nbj  )                           &
                   &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
                   &                    * tmask_i(nbi, nbj  )                           &
                   &                    * tmask_i(nbi, nbj+1)
-            ENDDO
-         ENDDO
+            END DO
+         END DO
          !
          IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
@@ -1344,13 +1296,12 @@
       ! Tidy up
       !--------
-      IF (nb_bdy>0) THEN
-         DEALLOCATE(nbidta, nbjdta, nbrdta)
-      ENDIF
-
-      CALL wrk_dealloc(jpi,jpj,zfmask) 
-
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
-
+      IF( nb_bdy>0 )   DEALLOCATE( nbidta, nbjdta, nbrdta )
+      !
+      CALL wrk_dealloc(jpi,jpj,   zfmask ) 
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_init')
+      !
    END SUBROUTINE bdy_init
+
 
    SUBROUTINE bdy_ctl_seg
@@ -1743,15 +1694,15 @@
       itest = 0
 
-      IF (cn_dyn2d(ib1)/=cn_dyn2d(ib2)) itest = itest + 1
-      IF (cn_dyn3d(ib1)/=cn_dyn3d(ib2)) itest = itest + 1
-      IF (cn_tra(ib1)/=cn_tra(ib2)) itest = itest + 1
-      !
-      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1
-      IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1
-      IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1
-      !
-      IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1   
-      !
-      IF ( itest>0 ) THEN
+      IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) )   itest = itest + 1
+      IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) )   itest = itest + 1
+      IF( cn_tra  (ib1) /= cn_tra  (ib2) )   itest = itest + 1
+      !
+      IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) )   itest = itest + 1
+      IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) )   itest = itest + 1
+      IF( nn_tra_dta  (ib1) /= nn_tra_dta  (ib2) )   itest = itest + 1
+      !
+      IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) )   itest = itest + 1   
+      !
+      IF( itest>0 ) THEN
          IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
          IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                  
Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 7351)
@@ -4,5 +4,5 @@
    !! Unstructured Open Boundary Cond. :  Library module of generic boundary algorithms.
    !!======================================================================
-   !! History :  3.6  !  2013     (D. Storkey) new module
+   !! History :  3.6  !  2013     (D. Storkey) original code
    !!----------------------------------------------------------------------
 #if defined key_bdy 
@@ -13,11 +13,12 @@
    !!   bdy_orlanski_3d
    !!----------------------------------------------------------------------
-   USE timing          ! Timing
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain
-   USE bdy_oce         ! ocean open boundary conditions
-   USE phycst          ! physical constants
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE in_out_manager  !
+   USE oce            ! ocean dynamics and tracers 
+   USE dom_oce        ! ocean space and time domain
+   USE bdy_oce        ! ocean open boundary conditions
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager !
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -45,11 +46,11 @@
       !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),            INTENT(in)     ::   idx      ! BDY indices
-      INTEGER,                    INTENT(in)     ::   igrd     ! grid index
-      REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   phib     ! model before 2D field
-      REAL(wp), DIMENSION(:,:),   INTENT(inout)  ::   phia     ! model after 2D field (to be updated)
-      REAL(wp), DIMENSION(:),     INTENT(in)     ::   phi_ext  ! external forcing data
-      LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version
-
+      TYPE(OBC_INDEX),          INTENT(in   ) ::   idx      ! BDY indices
+      INTEGER ,                 INTENT(in   ) ::   igrd     ! grid index
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phib     ! model before 2D field
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   phia     ! model after 2D field (to be updated)
+      REAL(wp), DIMENSION(:)  , INTENT(in   ) ::   phi_ext  ! external forcing data
+      LOGICAL ,                 INTENT(in   ) ::   ll_npo   ! switch for NPO version
+      !
       INTEGER  ::   jb                                     ! dummy loop indices
       INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
@@ -70,7 +71,7 @@
       REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives
       !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_2d')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_orlanski_2d')
+      !
       ! ----------------------------------!
       ! Orlanski boundary conditions     :!
@@ -79,25 +80,25 @@
       SELECT CASE(igrd)
          CASE(1)
-            pmask => tmask(:,:,1)
+            pmask      => tmask(:,:,1)
             pmask_xdif => umask(:,:,1)
             pmask_ydif => vmask(:,:,1)
-            pe_xdif => e1u(:,:)
-            pe_ydif => e2v(:,:)
+            pe_xdif    => e1u(:,:)
+            pe_ydif    => e2v(:,:)
             ii_offset = 0
             ij_offset = 0
          CASE(2)
-            pmask => umask(:,:,1)
+            pmask      => umask(:,:,1)
             pmask_xdif => tmask(:,:,1)
             pmask_ydif => fmask(:,:,1)
-            pe_xdif => e1t(:,:)
-            pe_ydif => e2f(:,:)
+            pe_xdif    => e1t(:,:)
+            pe_ydif    => e2f(:,:)
             ii_offset = 1
             ij_offset = 0
          CASE(3)
-            pmask => vmask(:,:,1)
+            pmask      => vmask(:,:,1)
             pmask_xdif => fmask(:,:,1)
             pmask_ydif => tmask(:,:,1)
-            pe_xdif => e1f(:,:)
-            pe_ydif => e2t(:,:)
+            pe_xdif    => e1f(:,:)
+            pe_ydif    => e2t(:,:)
             ii_offset = 0
             ij_offset = 1
@@ -188,6 +189,6 @@
       END DO
       !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_2d')
-
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_orlanski_2d')
+      !
    END SUBROUTINE bdy_orlanski_2d
 
@@ -204,11 +205,11 @@
       !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),            INTENT(in)     ::   idx      ! BDY indices
-      INTEGER,                    INTENT(in)     ::   igrd     ! grid index
-      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phib     ! model before 3D field
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated)
-      REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   phi_ext  ! external forcing data
-      LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version
-
+      TYPE(OBC_INDEX),            INTENT(in   ) ::   idx      ! BDY indices
+      INTEGER ,                   INTENT(in   ) ::   igrd     ! grid index
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phib     ! model before 3D field
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phia     ! model after 3D field (to be updated)
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   phi_ext  ! external forcing data
+      LOGICAL ,                   INTENT(in   ) ::   ll_npo   ! switch for NPO version
+      !
       INTEGER  ::   jb, jk                                 ! dummy loop indices
       INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
@@ -229,34 +230,34 @@
       REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives
       !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_3d')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_orlanski_3d')
+      !
       ! ----------------------------------!
       ! Orlanski boundary conditions     :!
       ! ----------------------------------! 
-     
+      !
       SELECT CASE(igrd)
          CASE(1)
-            pmask => tmask(:,:,:)
+            pmask      => tmask(:,:,:)
             pmask_xdif => umask(:,:,:)
             pmask_ydif => vmask(:,:,:)
-            pe_xdif => e1u(:,:)
-            pe_ydif => e2v(:,:)
+            pe_xdif    => e1u(:,:)
+            pe_ydif    => e2v(:,:)
             ii_offset = 0
             ij_offset = 0
          CASE(2)
-            pmask => umask(:,:,:)
+            pmask      => umask(:,:,:)
             pmask_xdif => tmask(:,:,:)
             pmask_ydif => fmask(:,:,:)
-            pe_xdif => e1t(:,:)
-            pe_ydif => e2f(:,:)
+            pe_xdif    => e1t(:,:)
+            pe_ydif    => e2f(:,:)
             ii_offset = 1
             ij_offset = 0
          CASE(3)
-            pmask => vmask(:,:,:)
+            pmask      => vmask(:,:,:)
             pmask_xdif => fmask(:,:,:)
             pmask_ydif => tmask(:,:,:)
-            pe_xdif => e1f(:,:)
-            pe_ydif => e2t(:,:)
+            pe_xdif    => e1f(:,:)
+            pe_ydif    => e2t(:,:)
             ii_offset = 0
             ij_offset = 1
@@ -349,7 +350,7 @@
          !
       END DO
-
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_3d')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_orlanski_3d')
+      !
    END SUBROUTINE bdy_orlanski_3d
 
Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 7351)
@@ -15,22 +15,21 @@
    !!   'key_bdy'     Open Boundary Condition
    !!----------------------------------------------------------------------
-   !!   PUBLIC
-   !!      bdytide_init     : read of namelist and initialisation of tidal harmonics data
-   !!      tide_update   : calculation of tidal forcing at each timestep
-   !!----------------------------------------------------------------------
-   USE timing          ! Timing
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain
-   USE iom
-   USE in_out_manager  ! I/O units
-   USE phycst          ! physical constants
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE bdy_par         ! Unstructured boundary parameters
-   USE bdy_oce         ! ocean open boundary conditions
-   USE daymod          ! calendar
-   USE wrk_nemo        ! Memory allocation
-   USE tideini
-!   USE tide_mod       ! Useless ??
-   USE fldread
+   !!   bdytide_init  : read of namelist and initialisation of tidal harmonics data
+   !!   tide_update   : calculation of tidal forcing at each timestep
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers 
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE bdy_par        ! Unstructured boundary parameters
+   USE bdy_oce        ! ocean open boundary conditions
+   USE tideini        ! 
+   USE daymod         ! calendar
+   !
+   USE in_out_manager ! I/O units
+   USE iom            ! xIO server
+   USE fldread        !
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE wrk_nemo       ! Memory allocation
+   USE timing         ! timing
 
    IMPLICIT NONE
@@ -42,10 +41,8 @@
 
    TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0       !: Tidal constituents : SSH0 (read in file)
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0         !: Tidal constituents : U0   (read in file)
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v0         !: Tidal constituents : V0   (read in file)
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH  (after nodal cor.)
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U    (after nodal cor.)
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0     !: Tidal constituents : SSH0   (read in file)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0, v0   !: Tidal constituents : U0, V0 (read in file)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh      !: Tidal constituents : SSH    (after nodal cor.)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u , v    !: Tidal constituents : U , V  (after nodal cor.)
    END TYPE TIDES_DATA
 
@@ -57,5 +54,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
+   !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -91,7 +88,7 @@
       NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
       !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('bdytide_init')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdytide_init')
+      !
       IF (nb_bdy>0) THEN
          IF(lwp) WRITE(numout,*)
@@ -263,5 +260,5 @@
             ENDIF ! ln_bdytide_2ddta=.true.
             !
-            IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files
+            IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files
                td%ssh0(:,:,2) = - td%ssh0(:,:,2)
                td%u0  (:,:,2) = - td%u0  (:,:,2)
@@ -274,17 +271,18 @@
             ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
             ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
-            dta_bdy_s(ib_bdy)%ssh(:) = 0.e0
-            dta_bdy_s(ib_bdy)%u2d(:) = 0.e0
-            dta_bdy_s(ib_bdy)%v2d(:) = 0.e0
+            dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
+            dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
+            dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
             !
          ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
          !
       END DO ! loop on ib_bdy
-
-      IF( nn_timing == 1 ) CALL timing_stop('bdytide_init')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdytide_init')
+      !
    END SUBROUTINE bdytide_init
 
-   SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset )
+
+   SUBROUTINE bdytide_update( kt, idx, dta, td, jit, time_offset )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdytide_update  ***
@@ -293,25 +291,25 @@
       !!                
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
-      TYPE(OBC_INDEX), INTENT( in )    ::   idx         ! OBC indices
-      TYPE(OBC_DATA),  INTENT(inout)   ::   dta         ! OBC external data
-      TYPE(TIDES_DATA),INTENT( inout ) ::   td          ! tidal harmonics data
-      INTEGER,INTENT(in),OPTIONAL      ::   jit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER,INTENT( in ), OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if jit
-                                                        ! is present then units = subcycle timesteps.
-                                                        ! time_offset = 0  => get data at "now"    time level
-                                                        ! time_offset = -1 => get data at "before" time level
-                                                        ! time_offset = +1 => get data at "after"  time level
-                                                        ! etc.
-      !!
-      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
-      INTEGER                          :: itide, igrd, ib   ! dummy loop indices
-      INTEGER                          :: time_add          ! time offset in units of timesteps
-      REAL(wp)                         :: z_arg, z_sarg, zflag, zramp      
+      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter
+      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices
+      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
+      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
+      INTEGER, OPTIONAL, INTENT(in   ) ::   jit         ! Barotropic timestep counter (for timesplitting option)
+      INTEGER, OPTIONAL, INTENT(in   ) ::   time_offset ! time offset in units of timesteps. NB. if jit
+      !                                                 ! is present then units = subcycle timesteps.
+      !                                                 ! time_offset = 0  => get data at "now"    time level
+      !                                                 ! time_offset = -1 => get data at "before" time level
+      !                                                 ! time_offset = +1 => get data at "after"  time level
+      !                                                 ! etc.
+      !
+      INTEGER  ::   itide, igrd, ib       ! dummy loop indices
+      INTEGER  ::   time_add              ! time offset in units of timesteps
+      INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays)
+      REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars    
       REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
       !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('bdytide_update')
-
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdytide_update')
+      !
       ilen0(1) =  SIZE(td%ssh(:,1,1))
       ilen0(2) =  SIZE(td%u(:,1,1))
@@ -323,7 +321,7 @@
       ENDIF
 
-      IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. zflag==1 ) THEN
+      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
         !
-        kt_tide = kt
+        kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
         !
         IF(lwp) THEN
@@ -374,5 +372,5 @@
       END DO
       !
-      IF( nn_timing == 1 ) CALL timing_stop('bdytide_update')
+      IF( nn_timing == 1 )   CALL timing_stop('bdytide_update')
       !
    END SUBROUTINE bdytide_update
@@ -385,23 +383,23 @@
       !!                
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
-      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit
-                                                        ! is present then units = subcycle timesteps.
-                                                        ! time_offset = 0  => get data at "now"    time level
-                                                        ! time_offset = -1 => get data at "before" time level
-                                                        ! time_offset = +1 => get data at "after"  time level
-                                                        ! etc.
-      !!
-      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step
-      INTEGER, DIMENSION(jpbgrd) :: ilen0 
-      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts
-      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices
-      INTEGER  :: time_add                                    ! time offset in units of timesteps
-      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist      
-      !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides')
-
+      INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
+      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
+      INTEGER, OPTIONAL, INTENT(in) ::   time_offset ! time offset in units of timesteps. NB. if kit
+      !                                              ! is present then units = subcycle timesteps.
+      !                                              ! time_offset = 0  => get data at "now"    time level
+      !                                              ! time_offset = -1 => get data at "before" time level
+      !                                              ! time_offset = +1 => get data at "after"  time level
+      !                                              ! etc.
+      !
+      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
+      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
+      INTEGER  ::   time_add                  ! time offset in units of timesteps
+      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
+      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
+      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist      
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_dta_tides')
+      !
       lk_first_btstp=.TRUE.
       IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
@@ -438,7 +436,7 @@
             ! We refresh nodal factors every day below
             ! This should be done somewhere else
-            IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. lk_first_btstp ) THEN
-               !
-               kt_tide = kt               
+            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
+               !
+               kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
                !
                IF(lwp) THEN
@@ -502,22 +500,23 @@
    END SUBROUTINE bdy_dta_tides
 
+
    SUBROUTINE tide_init_elevation( idx, td )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE tide_init_elevation  ***
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
-      TYPE(TIDES_DATA),INTENT( inout )   ::   td      ! tidal harmonics data
-      !! * Local declarations
-      INTEGER, DIMENSION(1)            ::   ilen0       !: length of boundary data (from OBC arrays)
+      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
+      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
+      !
+      INTEGER ::   itide, igrd, ib       ! dummy loop indices
+      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
-      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
-
+      !!----------------------------------------------------------------------
+      !
       igrd=1   
                               ! SSH on tracer grid.
-   
       ilen0(1) =  SIZE(td%ssh0(:,1,1))
-
-      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
-
+      !
+      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
+      !
       DO itide = 1, nb_harmo
          DO ib = 1, ilen0(igrd)
@@ -534,27 +533,29 @@
          ENDDO
       END DO
-
-      DEALLOCATE(mod_tide,phi_tide)
-
+      !
+      DEALLOCATE( mod_tide, phi_tide )
+      !
    END SUBROUTINE tide_init_elevation
 
+
    SUBROUTINE tide_init_velocities( idx, td )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE tide_init_elevation  ***
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
-      TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data
-      !! * Local declarations
-      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
+      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
+      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
+      !
+      INTEGER ::   itide, igrd, ib       ! dummy loop indices
+      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
-      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
-
+      !!----------------------------------------------------------------------
+      !
       ilen0(2) =  SIZE(td%u0(:,1,1))
       ilen0(3) =  SIZE(td%v0(:,1,1))
-
+      !
       igrd=2                                 ! U grid.
-
-      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
-
+      !
+      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
+      !
       DO itide = 1, nb_harmo
          DO ib = 1, ilen0(igrd)
@@ -571,10 +572,10 @@
          ENDDO
       END DO
-
-      DEALLOCATE(mod_tide,phi_tide)
-
+      !
+      DEALLOCATE( mod_tide , phi_tide )
+      !
       igrd=3                                 ! V grid.
-
-      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
+      !
+      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
 
       DO itide = 1, nb_harmo
@@ -592,8 +593,9 @@
          ENDDO
       END DO
-
-      DEALLOCATE(mod_tide,phi_tide)
-
+      !
+      DEALLOCATE( mod_tide, phi_tide )
+      !
   END SUBROUTINE tide_init_velocities
+
 #else
    !!----------------------------------------------------------------------
Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 7351)
@@ -16,19 +16,19 @@
    !!   bdy_tra_frs        : Apply Flow Relaxation Scheme
    !!----------------------------------------------------------------------
-   USE timing          ! Timing
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE bdy_oce         ! ocean open boundary conditions
-   USE bdylib          ! for orlanski library routines
-   USE bdydta, ONLY:   bf
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE in_out_manager  ! I/O manager
-
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables 
+   USE bdy_oce        ! ocean open boundary conditions
+   USE bdylib         ! for orlanski library routines
+   USE bdydta   , ONLY:   bf   ! 
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE timing         ! Timing
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC bdy_tra      ! routine called in tranxt.F90 
-   PUBLIC bdy_tra_dmp  ! routine called in step.F90 
+   PUBLIC   bdy_tra      ! called in tranxt.F90 
+   PUBLIC   bdy_tra_dmp  ! called in step.F90 
 
    !!----------------------------------------------------------------------
@@ -46,36 +46,29 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) :: kt     ! Main time step counter
-      !!
-      INTEGER               :: ib_bdy ! Loop index
+      INTEGER, INTENT(in) ::   kt   ! Main time step counter
+      !
+      INTEGER ::   ib_bdy   ! Loop index
+      !!----------------------------------------------------------------------
 
       DO ib_bdy=1, nb_bdy
-
+         !
          SELECT CASE( cn_tra(ib_bdy) )
-         CASE('none')
-            CYCLE
-         CASE('frs')
-            CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
-         CASE('specified')
-            CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
-         CASE('neumann')
-            CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
-         CASE('orlanski')
-            CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.false. )
-         CASE('orlanski_npo')
-            CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.true. )
-         CASE('runoff')
-            CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
-         CASE DEFAULT
-            CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
+         CASE('none'        )   ;   CYCLE
+         CASE('frs'         )   ;   CALL bdy_tra_frs     ( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE('specified'   )   ;   CALL bdy_tra_spe     ( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE('neumann'     )   ;   CALL bdy_tra_nmn     ( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE('orlanski'    )   ;   CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.false. )
+         CASE('orlanski_npo')   ;   CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.true. )
+         CASE('runoff'      )   ;   CALL bdy_tra_rnf     ( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE DEFAULT           ;   CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
          END SELECT
          ! Boundary points should be updated
          CALL lbc_bdy_lnk( tsa(:,:,:,jp_tem), 'T', 1., ib_bdy )
          CALL lbc_bdy_lnk( tsa(:,:,:,jp_sal), 'T', 1., ib_bdy )
-      ENDDO
-      !
-
+      END DO
+      !
    END SUBROUTINE bdy_tra
 
+
    SUBROUTINE bdy_tra_frs( idx, dta, kt )
       !!----------------------------------------------------------------------
@@ -86,8 +79,8 @@
       !! Reference : Engedahl H., 1995, Tellus, 365-382.
       !!----------------------------------------------------------------------
-      INTEGER,         INTENT(in) ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
-      !! 
+      INTEGER,         INTENT(in) ::   kt    !
+      TYPE(OBC_INDEX), INTENT(in) ::   idx   ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta   ! OBC external data
+      !
       REAL(wp) ::   zwgt           ! boundary weight
       INTEGER  ::   ib, ik, igrd   ! dummy loop indices
@@ -95,5 +88,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_frs')
+      IF( nn_timing == 1 )   CALL timing_start('bdy_tra_frs')
       !
       igrd = 1                       ! Everything is at T-points here
@@ -108,10 +101,11 @@
       END DO 
       !
-      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_frs')
+      IF( kt .eq. nit000 )   CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_tra_frs')
       !
    END SUBROUTINE bdy_tra_frs
-  
+
+
    SUBROUTINE bdy_tra_spe( idx, dta, kt )
       !!----------------------------------------------------------------------
@@ -121,8 +115,8 @@
       !! 
       !!----------------------------------------------------------------------
-      INTEGER,         INTENT(in) ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
-      !! 
+      INTEGER,         INTENT(in) ::   kt    !
+      TYPE(OBC_INDEX), INTENT(in) ::   idx   ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta   ! OBC external data
+      !
       REAL(wp) ::   zwgt           ! boundary weight
       INTEGER  ::   ib, ik, igrd   ! dummy loop indices
@@ -142,10 +136,11 @@
       END DO
       !
-      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_spe')
+      IF( kt == nit000 )   CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_tra_spe')
       !
    END SUBROUTINE bdy_tra_spe
 
+
    SUBROUTINE bdy_tra_nmn( idx, dta, kt )
       !!----------------------------------------------------------------------
@@ -155,8 +150,8 @@
       !! 
       !!----------------------------------------------------------------------
-      INTEGER,         INTENT(in) ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
-      !! 
+      INTEGER,         INTENT(in) ::   kt    ! 
+      TYPE(OBC_INDEX), INTENT(in) ::   idx   ! OBC indices
+      TYPE(OBC_DATA) , INTENT(in) ::   dta   ! OBC external data
+      !
       REAL(wp) ::   zwgt           ! boundary weight
       INTEGER  ::   ib, ik, igrd   ! dummy loop indices
@@ -164,5 +159,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_nmn')
+      IF( nn_timing == 1 )   CALL timing_start('bdy_tra_nmn')
       !
       igrd = 1                       ! Everything is at T-points here
@@ -196,7 +191,7 @@
       END DO
       !
-      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_nmn')
+      IF( kt == nit000 )   CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_tra_nmn')
       !
    END SUBROUTINE bdy_tra_nmn
@@ -213,11 +208,11 @@
       !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
-      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version
-
+      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
+      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
+      LOGICAL        , INTENT(in) ::   ll_npo  ! switch for NPO version
+      !
       INTEGER  ::   igrd                                    ! grid index
       !!----------------------------------------------------------------------
-
+      !
       IF( nn_timing == 1 ) CALL timing_start('bdy_tra_orlanski')
       !
@@ -230,7 +225,6 @@
       CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_sal), tsa(:,:,:,jp_sal), dta%sal, ll_npo )
       !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_orlanski')
-      !
-
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_tra_orlanski')
+      !
    END SUBROUTINE bdy_tra_orlanski
 
@@ -245,8 +239,8 @@
       !! 
       !!----------------------------------------------------------------------
-      INTEGER,         INTENT(in) ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
-      !! 
+      INTEGER        , INTENT(in) ::   kt    ! 
+      TYPE(OBC_INDEX), INTENT(in) ::   idx   ! OBC indices
+      TYPE(OBC_DATA) , INTENT(in) ::   dta   ! OBC external data
+      !
       REAL(wp) ::   zwgt           ! boundary weight
       INTEGER  ::   ib, ik, igrd   ! dummy loop indices
@@ -254,5 +248,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_rnf')
+      IF( nn_timing == 1 )   CALL timing_start('bdy_tra_rnf')
       !
       igrd = 1                       ! Everything is at T-points here
@@ -268,10 +262,11 @@
       END DO
       !
-      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_rnf')
+      IF( kt == nit000 )   CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_tra_rnf')
       !
    END SUBROUTINE bdy_tra_rnf
 
+
    SUBROUTINE bdy_tra_dmp( kt )
       !!----------------------------------------------------------------------
@@ -281,6 +276,6 @@
       !! 
       !!----------------------------------------------------------------------
-      INTEGER,         INTENT(in) ::   kt
-      !! 
+      INTEGER, INTENT(in) ::   kt   !
+      !
       REAL(wp) ::   zwgt           ! boundary weight
       REAL(wp) ::   zta, zsa, ztime
@@ -290,8 +285,8 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_dmp')
-      !
-      DO ib_bdy=1, nb_bdy
-         IF ( ln_tra_dmp(ib_bdy) ) THEN
+      IF( nn_timing == 1 )   CALL timing_start('bdy_tra_dmp')
+      !
+      DO ib_bdy = 1, nb_bdy
+         IF( ln_tra_dmp(ib_bdy) ) THEN
             igrd = 1                       ! Everything is at T-points here
             DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
@@ -307,7 +302,7 @@
             END DO
          ENDIF
-      ENDDO
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp')
+      END DO
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_tra_dmp')
       !
    END SUBROUTINE bdy_tra_dmp
@@ -325,5 +320,4 @@
       WRITE(*,*) 'bdy_tra_dmp: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_tra_dmp
-
 #endif
 
Index: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
===================================================================
--- branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 7349)
+++ branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 7351)
@@ -10,29 +10,27 @@
    !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
    !!----------------------------------------------------------------------
-#if   defined key_bdy
+#if defined key_bdy
    !!----------------------------------------------------------------------
-   !!   'key_bdy'      unstructured open boundary conditions
+   !!   'key_bdy'                     unstructured open boundary conditions
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers 
-   USE bdy_oce         ! ocean open boundary conditions
-   USE sbc_oce         ! ocean surface boundary conditions
-   USE dom_oce         ! ocean space and time domain 
-   USE phycst          ! physical constants
-   USE sbcisf          ! ice shelf
+   USE oce            ! ocean dynamics and tracers 
+   USE bdy_oce        ! ocean open boundary conditions
+   USE sbc_oce        ! ocean surface boundary conditions
+   USE dom_oce        ! ocean space and time domain 
+   USE phycst         ! physical constants
+   USE sbcisf         ! ice shelf
    !
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! for mppsum
-   USE timing          ! Timing
-   USE lib_fortran     ! Fortran routines library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! for mppsum
+   USE timing         ! Timing
+   USE lib_fortran    ! Fortran routines library
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC bdy_vol      
+   PUBLIC   bdy_vol    ! called by ???
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$ 
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -45,8 +43,8 @@
       !!
       !! ** Purpose :   This routine controls the volume of the system. 
-      !!      A correction velocity is calculated
-      !!      to correct the total transport through the unstructured OBC. 
+      !!      A correction velocity is calculated to correct the total transport 
+      !!      through the unstructured OBC. 
       !!      The total depth used is constant (H0) to be consistent with the 
-      !!      linear free surface coded in OPA 8.2
+      !!      linear free surface coded in OPA 8.2    <<<=== !!gm  ???? true ????
       !!
       !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
@@ -72,6 +70,6 @@
       !!            (set nn_volctl to 1 in tne namelist for this option)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      !!
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
       INTEGER  ::   ji, jj, jk, jb, jgrd
       INTEGER  ::   ib_bdy, ii, ij
@@ -93,5 +91,5 @@
       ! -----------------------------------------------------------------------
 !!gm replace these lines :
-      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0
+      z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0
       IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
 !!gm   by :
@@ -110,5 +108,5 @@
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
+               zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * e3u_n(ii,ij,jk)
             END DO
          END DO
@@ -118,5 +116,5 @@
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
+               zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * e3v_n(ii,ij,jk) 
             END DO
          END DO
@@ -127,6 +125,6 @@
       ! The normal velocity correction
       ! ------------------------------
-      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot 
-      ELSE                      ;   zubtpecor =   zubtpecor             / bdysurftot
+      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot 
+      ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot
       END IF
 
@@ -143,5 +141,5 @@
                ij = idx%nbj(jb,jgrd)
                ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb,jgrd) * zubtpecor * umask(ii,ij,jk)
-               ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
+               ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * e3u_n(ii,ij,jk)
             END DO
          END DO
@@ -152,5 +150,5 @@
                ij = idx%nbj(jb,jgrd)
                va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb,jgrd) * zubtpecor * vmask(ii,ij,jk)
-               ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk)
+               ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * e3v_n(ii,ij,jk)
             END DO
          END DO
@@ -161,5 +159,5 @@
       ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
       ! ------------------------------------------------------
-      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
+      IF( lwp .AND. MOD( kt, nwrite ) == 0 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
@@ -171,5 +169,5 @@
       END IF 
       !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_vol')
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_vol')
       !
       END IF ! ln_vol
