Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 3991)
@@ -29,9 +29,23 @@
       REAL   , POINTER, DIMENSION(:,:)   ::  nbw
       REAL   , POINTER, DIMENSION(:,:)   ::  nbd
-      REAL   , POINTER, DIMENSION(:)     ::  flagu
-      REAL   , POINTER, DIMENSION(:)     ::  flagv
+      REAL   , POINTER, DIMENSION(:,:)   ::  nbdout
+      REAL   , POINTER, DIMENSION(:,:)   ::  flagu
+      REAL   , POINTER, DIMENSION(:,:)   ::  flagv
    END TYPE OBC_INDEX
 
+   !! Logicals in OBC_DATA structure are true if the chosen algorithm requires this
+   !! field as external data. If true the data can come from external files
+   !! or model initial conditions. If false then no "external" data array
+   !! is required for this field. 
+
    TYPE, PUBLIC ::   OBC_DATA     !: Storage for external data
+      INTEGER,       DIMENSION(2)     ::  nread
+      LOGICAL                         ::  ll_ssh
+      LOGICAL                         ::  ll_u2d
+      LOGICAL                         ::  ll_v2d
+      LOGICAL                         ::  ll_u3d
+      LOGICAL                         ::  ll_v3d
+      LOGICAL                         ::  ll_tem
+      LOGICAL                         ::  ll_sal
       REAL, POINTER, DIMENSION(:)     ::  ssh
       REAL, POINTER, DIMENSION(:)     ::  u2d
@@ -42,4 +56,7 @@
       REAL, POINTER, DIMENSION(:,:)   ::  sal
 #if defined key_lim2
+      LOGICAL                         ::  ll_frld
+      LOGICAL                         ::  ll_hicif
+      LOGICAL                         ::  ll_hsnif
       REAL, POINTER, DIMENSION(:)     ::  frld
       REAL, POINTER, DIMENSION(:)     ::  hicif
@@ -63,23 +80,24 @@
    INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P 
    !                                                        !  = 1 the volume will be constant during all the integration.
-   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d                 ! Choice of boundary condition for barotropic variables (U,V,SSH)
-   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d_dta      	    !: = 0 use the initial state as bdy dta ; 
+   CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_dyn2d       ! Choice of boundary condition for barotropic variables (U,V,SSH)
+   INTEGER, DIMENSION(jp_bdy)           ::   nn_dyn2d_dta   !: = 0 use the initial state as bdy dta ; 
                                                             !: = 1 read it in a NetCDF file
                                                             !: = 2 read tidal harmonic forcing from a NetCDF file
                                                             !: = 3 read external data AND tidal harmonic forcing from NetCDF files
-   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn3d                 ! Choice of boundary condition for baroclinic velocities 
-   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn3d_dta      	    !: = 0 use the initial state as bdy dta ; 
+   CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_dyn3d       ! Choice of boundary condition for baroclinic velocities 
+   INTEGER, DIMENSION(jp_bdy)           ::   nn_dyn3d_dta   !: = 0 use the initial state as bdy dta ; 
                                                             !: = 1 read it in a NetCDF file
-   INTEGER, DIMENSION(jp_bdy) ::   nn_tra                   ! Choice of boundary condition for active tracers (T and S)
-   INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta      	    !: = 0 use the initial state as bdy dta ; 
+   CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_tra         ! Choice of boundary condition for active tracers (T and S)
+   INTEGER, DIMENSION(jp_bdy)           ::   nn_tra_dta     !: = 0 use the initial state as bdy dta ; 
                                                             !: = 1 read it in a NetCDF file
    LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping
    LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping
    REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days
+   REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp_out          !: Damping time scale in days at radiation outflow points
 
 #if defined key_lim2
-   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables 
-   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2_dta          !: = 0 use the initial state as bdy dta ; 
-                                                            !: = 1 read it in a NetCDF file
+   CHARACTER(len=20), DIMENSION(jp_bdy) ::   nn_ice_lim2      ! Choice of boundary condition for sea ice variables 
+   INTEGER, DIMENSION(jp_bdy)           ::   nn_ice_lim2_dta  !: = 0 use the initial state as bdy dta ; 
+                                                              !: = 1 read it in a NetCDF file
 #endif
    !
@@ -88,15 +106,16 @@
    !! Global variables
    !!----------------------------------------------------------------------
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdytmask   !: Mask defining computational domain at T-points
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyumask   !: Mask defining computational domain at U-points
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyvmask   !: Mask defining computational domain at V-points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdytmask   !: Mask defining computational domain at T-points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyumask   !: Mask defining computational domain at U-points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyvmask   !: Mask defining computational domain at V-points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyfmask   !: Mask defining computational domain at F-points
 
    REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary
 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh       !: 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   phur       !: 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr       !: Pointers for barotropic fields 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   pu2d       !: 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   pv2d       !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh                  !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   phur                  !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr                  !: Pointers for barotropic fields 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pub2d, pun2d, pua2d   !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pvb2d, pvn2d, pva2d   !: 
 
    !!----------------------------------------------------------------------
@@ -109,5 +128,5 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy)
    TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process)
-   TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process)
+   TYPE(OBC_DATA) , DIMENSION(jp_bdy), TARGET      ::   dta_bdy           !: bdy external data (local process)
 
    !!----------------------------------------------------------------------
@@ -125,5 +144,5 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),                    &  
+      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj), bdyfmask(jpi,jpj),     &  
          &      STAT=bdy_oce_alloc )
          !
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90	(revision 3991)
@@ -23,11 +23,6 @@
 # endif
    INTEGER, PUBLIC, PARAMETER ::   jp_bdy  = 10       !: Maximum number of bdy sets
-   INTEGER, PUBLIC, PARAMETER ::   jpbtime = 1000     !: Max number of time dumps per file
    INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 3	      !: Number of horizontal grid types used  (T, U, V)
 
-   !! Flags for choice of schemes
-   INTEGER, PUBLIC, PARAMETER ::   jp_none         = 0        !: Flag for no open boundary condition
-   INTEGER, PUBLIC, PARAMETER ::   jp_frs          = 1        !: Flag for Flow Relaxation Scheme
-   INTEGER, PUBLIC, PARAMETER ::   jp_flather      = 2        !: Flag for Flather
 #else
    !!----------------------------------------------------------------------
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3991)
@@ -80,4 +80,5 @@
       INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
       INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
+      TYPE(OBC_DATA), POINTER             ::   dta              ! short cut
       !!
       !!---------------------------------------------------------------------------
@@ -91,15 +92,14 @@
          ! Calculate depth-mean currents
          !-----------------------------
-         CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 
-
-         pu2d(:,:) = 0.e0
-         pv2d(:,:) = 0.e0
-
+         CALL wrk_alloc(jpi,jpj,pun2d,pvn2d) 
+
+         pun2d(:,:) = 0.e0
+         pvn2d(:,:) = 0.e0
          DO ik = 1, jpkm1   !! Vertically integrated momentum trends
-             pu2d(:,:) = pu2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik)
-             pv2d(:,:) = pv2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik)
+             pun2d(:,:) = pun2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik)
+             pvn2d(:,:) = pvn2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik)
          END DO
-         pu2d(:,:) = pu2d(:,:) * hur(:,:)
-         pv2d(:,:) = pv2d(:,:) * hvr(:,:)
+         pun2d(:,:) = pun2d(:,:) * hur(:,:)
+         pvn2d(:,:) = pvn2d(:,:) * hvr(:,:)
          
          DO ib_bdy = 1, nb_bdy
@@ -107,71 +107,109 @@
             nblen => idx_bdy(ib_bdy)%nblen
             nblenrim => idx_bdy(ib_bdy)%nblenrim
-
-            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 
+            dta => dta_bdy(ib_bdy)
+
+            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 
                ilen1(:) = nblen(:)
-               igrd = 1
-               DO ib = 1, ilen1(igrd)
-                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                  dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
-               END DO 
-               igrd = 2
-               DO ib = 1, ilen1(igrd)
-                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                  dta_bdy(ib_bdy)%u2d(ib) = pu2d(ii,ij) * umask(ii,ij,1)         
-               END DO 
-               igrd = 3
-               DO ib = 1, ilen1(igrd)
-                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                  dta_bdy(ib_bdy)%v2d(ib) = pv2d(ii,ij) * vmask(ii,ij,1)         
-               END DO 
-            ENDIF
-
-            IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
-               ilen1(:) = nblen(:)
-               igrd = 2 
-               DO ib = 1, ilen1(igrd)
-                  DO ik = 1, jpkm1
+               IF( dta%ll_ssh ) THEN 
+                  igrd = 1
+                  DO ib = 1, ilen1(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - pu2d(ii,ij) ) * umask(ii,ij,ik)         
-                  END DO
-               END DO 
-               igrd = 3 
-               DO ib = 1, ilen1(igrd)
-                  DO ik = 1, jpkm1
+                     dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
+                  END DO 
+               END IF
+               IF( dta%ll_u2d ) THEN 
+                  igrd = 2
+                  DO ib = 1, ilen1(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pv2d(ii,ij) ) * vmask(ii,ij,ik)         
-                     END DO
-               END DO 
-            ENDIF
-
-            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
-               ilen1(:) = nblen(:)
-               igrd = 1                       ! Everything is at T-points here
-               DO ib = 1, ilen1(igrd)
-                  DO ik = 1, jpkm1
+                     dta_bdy(ib_bdy)%u2d(ib) = pun2d(ii,ij) * umask(ii,ij,1)         
+                  END DO 
+               END IF
+               IF( dta%ll_v2d ) THEN 
+                  igrd = 3
+                  DO ib = 1, ilen1(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
-                     dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
-                  END DO
-               END DO 
-            ENDIF
-
-#if defined key_lim2
-            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
+                     dta_bdy(ib_bdy)%v2d(ib) = pvn2d(ii,ij) * vmask(ii,ij,1)         
+                  END DO 
+               END IF
+            ENDIF
+
+            IF( nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
                ilen1(:) = nblen(:)
-               igrd = 1                       ! Everything is at T-points here
-               DO ib = 1, ilen1(igrd)
-                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                  dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)         
-                  dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)         
-                  dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)         
-               END DO 
+               IF( dta%ll_u3d ) THEN 
+                  igrd = 2 
+                  DO ib = 1, ilen1(igrd)
+                     DO ik = 1, jpkm1
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - pun2d(ii,ij) ) * umask(ii,ij,ik)         
+                     END DO
+                  END DO 
+               END IF
+               IF( dta%ll_v3d ) THEN 
+                  igrd = 3 
+                  DO ib = 1, ilen1(igrd)
+                     DO ik = 1, jpkm1
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pvn2d(ii,ij) ) * vmask(ii,ij,ik)         
+                        END DO
+                  END DO 
+               END IF
+            ENDIF
+
+            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
+               ilen1(:) = nblen(:)
+               IF( dta%ll_tem ) THEN
+                  igrd = 1 
+                  DO ib = 1, ilen1(igrd)
+                     DO ik = 1, jpkm1
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
+                     END DO
+                  END DO 
+               END IF
+               IF( dta%ll_sal ) THEN
+                  igrd = 1 
+                  DO ib = 1, ilen1(igrd)
+                     DO ik = 1, jpkm1
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
+                     END DO
+                  END DO 
+               END IF
+            ENDIF
+
+#if defined key_lim2
+            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
+               ilen1(:) = nblen(:)
+               IF( dta%ll_frld ) THEN
+                  igrd = 1 
+                  DO ib = 1, ilen1(igrd)
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)         
+                  END DO 
+               END IF
+               IF( dta%ll_hicif ) THEN
+                  igrd = 1 
+                  DO ib = 1, ilen1(igrd)
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)         
+                  END DO 
+               END IF
+               IF( dta%ll_hsnif ) THEN
+                  igrd = 1 
+                  DO ib = 1, ilen1(igrd)
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)         
+                  END DO 
+               END IF
             ENDIF
 #endif
@@ -179,5 +217,5 @@
          ENDDO ! ib_bdy
 
-         CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 
+         CALL wrk_dealloc(jpi,jpj,pun2d,pvn2d) 
 
       ENDIF ! kt .eq. nit000
@@ -188,4 +226,5 @@
       jstart = 1
       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
       
@@ -193,22 +232,18 @@
                ! Update barotropic boundary conditions only
                ! jit is optional argument for fld_read and bdytide_update
-               IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
+               IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
                   IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                     dta_bdy(ib_bdy)%ssh(:) = 0.0
-                     dta_bdy(ib_bdy)%u2d(:) = 0.0
-                     dta_bdy(ib_bdy)%v2d(:) = 0.0
+                     IF( dta%ll_ssh ) dta%ssh(:) = 0.0
+                     IF( dta%ll_u2d ) dta%u2d(:) = 0.0
+                     IF( dta%ll_u3d ) dta%v2d(:) = 0.0
                   ENDIF
-                  IF (nn_tra(ib_bdy).ne.4) THEN
-                     IF( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR.  &
-                       & (ln_full_vel_array(ib_bdy) .AND. nn_dyn3d_dta(ib_bdy).eq.1) )THEN
-
-                        ! For the runoff case, no need to update the forcing (already done in the baroclinic part)
-                        jend = nb_bdy_fld(ib_bdy)
-                        IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2
+                  IF (cn_tra(ib_bdy) /= 'runoff') THEN
+                     IF( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 ) THEN
+
+                        jend = jstart + dta%nread(2) - 1
                         CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
                                      & kit=jit, kt_offset=time_offset )
-                        IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2
-
-                        ! If full velocities in boundary data then split into barotropic and baroclinic data
+
+                        ! 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.  &
@@ -216,35 +251,29 @@
 
                            igrd = 2                      ! zonal velocity
-                           dta_bdy(ib_bdy)%u2d(:) = 0.0
+                           dta%u2d(:) = 0.0
                            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                               ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
                               ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
                               DO ik = 1, jpkm1
-                                 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &
-                       &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik)
+                                 dta%u2d(ib) = dta%u2d(ib) &
+                       &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
                               END DO
-                              dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)
-                              DO ik = 1, jpkm1
-                                 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)
-                              END DO
+                              dta%u2d(ib) =  dta%u2d(ib) * hur(ii,ij)
                            END DO
                            igrd = 3                      ! meridional velocity
-                           dta_bdy(ib_bdy)%v2d(:) = 0.0
+                           dta%v2d(:) = 0.0
                            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                               ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
                               ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
                               DO ik = 1, jpkm1
-                                 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &
-                       &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik)
+                                 dta%v2d(ib) = dta%v2d(ib) &
+                       &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
                               END DO
-                              dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)
-                              DO ik = 1, jpkm1
-                                 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)
-                              END DO
+                              dta%v2d(ib) =  dta%v2d(ib) * hvr(ii,ij)
                            END DO
                         ENDIF                    
                      ENDIF
                      IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                        CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   & 
+                        CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta, td=tides(ib_bdy),   & 
                           &                 jit=jit, time_offset=time_offset )
                      ENDIF
@@ -252,5 +281,5 @@
                ENDIF
             ELSE
-               IF (nn_tra(ib_bdy).eq.4) then      ! runoff condition
+               IF (cn_tra(ib_bdy) == 'runoff') then      ! runoff condition
                   jend = nb_bdy_fld(ib_bdy)
                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  &
@@ -261,5 +290,5 @@
                      ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
                      ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
+                     dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
                   END DO
                   !
@@ -268,14 +297,14 @@
                      ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
                      ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
+                     dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
                   END DO
                ELSE
-                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                     dta_bdy(ib_bdy)%ssh(:) = 0.0
-                     dta_bdy(ib_bdy)%u2d(:) = 0.0
-                     dta_bdy(ib_bdy)%v2d(:) = 0.0
+                  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
                   ENDIF
-                  IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data
-                     jend = nb_bdy_fld(ib_bdy)
+                  IF( dta%nread(1) .gt. 0 ) THEN ! update external data
+                     jend = jstart + dta%nread(1) - 1
                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
                                   & map=nbmap_ptr(jstart:jend), kt_offset=time_offset )
@@ -286,39 +315,39 @@
                     &   nn_dyn3d_dta(ib_bdy) .EQ. 1 ) ) THEN
                      igrd = 2                      ! zonal velocity
-                     dta_bdy(ib_bdy)%u2d(:) = 0.0
+                     dta%u2d(:) = 0.0
                      DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                         ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
                         ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
                         DO ik = 1, jpkm1
-                           dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &
-                 &                       + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik)
+                           dta%u2d(ib) = dta%u2d(ib) &
+                 &                       + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
                         END DO
-                        dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)
+                        dta%u2d(ib) =  dta%u2d(ib) * hur(ii,ij)
                         DO ik = 1, jpkm1
-                           dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)
+                           dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib)
                         END DO
                      END DO
                      igrd = 3                      ! meridional velocity
-                     dta_bdy(ib_bdy)%v2d(:) = 0.0
+                     dta%v2d(:) = 0.0
                      DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                         ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
                         ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
                         DO ik = 1, jpkm1
-                           dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &
-                 &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik)
+                           dta%v2d(ib) = dta%v2d(ib) &
+                 &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
                         END DO
-                        dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)
+                        dta%v2d(ib) =  dta%v2d(ib) * hvr(ii,ij)
                         DO ik = 1, jpkm1
-                           dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)
+                           dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib)
                         END DO
                      END DO
                   ENDIF
-                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                     CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  &
+                  IF( cn_dyn2d(ib_bdy) /= 'none' .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                     CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta,  &
                                         & td=tides(ib_bdy), time_offset=time_offset )
                   ENDIF
                ENDIF
             ENDIF
-            jstart = jend+1
+            jstart = jstart + dta%nread(1)
          END IF ! nn_dta(ib_bdy) = 1
       END DO  ! ib_bdy
@@ -326,5 +355,5 @@
       IF ( ln_apr_obc ) THEN
          DO ib_bdy = 1, nb_bdy
-            IF (nn_tra(ib_bdy).NE.4)THEN
+            IF (cn_tra(ib_bdy) /= 'runoff')THEN
                igrd = 1                      ! meridional velocity
                DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
@@ -349,5 +378,5 @@
       !!                for open boundary conditions
       !!
-      !! ** Method  :   Use fldread.F90
+      !! ** Method  :   
       !!                
       !!----------------------------------------------------------------------
@@ -361,9 +390,9 @@
                                                                 ! =F => baroclinic velocities in 3D boundary data
       INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays
-      INTEGER,              DIMENSION(jpbgrd) ::  ilen0         ! size of local arrays
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ibdy           ! bdy set for a particular jfld
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V)
       INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts
+      TYPE(OBC_DATA), POINTER                ::   dta           ! short cut
       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures
       TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! 
@@ -403,15 +432,15 @@
       nb_bdy_fld(:) = 0
       DO ib_bdy = 1, nb_bdy         
-         IF( nn_dyn2d(ib_bdy) .gt. 0 .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) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
          ENDIF
-         IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN
+         IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
          ENDIF
-         IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN
+         IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
          ENDIF
 #if defined key_lim2
-         IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN
+         IF( cn_ice_lim2(ib_bdy) /= 'none' .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN
             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
          ENDIF
@@ -471,10 +500,13 @@
             nblen => idx_bdy(ib_bdy)%nblen
             nblenrim => idx_bdy(ib_bdy)%nblenrim
+            dta => dta_bdy(ib_bdy)
+            dta%nread(2) = 0
 
             ! Only read in necessary fields for this set.
             ! Important that barotropic variables come first.
-            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
-
-               IF( nn_dyn2d(ib_bdy) .ne. jp_frs .and. nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading
+            IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN 
+
+               IF( dta%ll_ssh ) THEN 
+                  if(lwp) write(numout,*) '++++++ reading in ssh field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
@@ -483,7 +515,9 @@
                   ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
-               ENDIF
-
-               IF( .not. ln_full_vel_array(ib_bdy) ) THEN
+                  dta%nread(2) = dta%nread(2) + 1
+               ENDIF
+
+               IF( dta%ll_u2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN
+                  if(lwp) write(numout,*) '++++++ reading in u2d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_u2d
@@ -492,5 +526,9 @@
                   ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
-
+                  dta%nread(2) = dta%nread(2) + 1
+               ENDIF
+
+               IF( dta%ll_v2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN
+                  if(lwp) write(numout,*) '++++++ reading in v2d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_v2d
@@ -499,45 +537,60 @@
                   ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
-               ENDIF
-
-            ENDIF
-
-            ! baroclinic velocities
-            IF( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) .or. &
-           &      ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.  &
-           &        ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_u3d
-               ibdy(jfld) = ib_bdy
-               igrid(jfld) = 2
-               ilen1(jfld) = nblen(igrid(jfld))
-               ilen3(jfld) = jpk
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_v3d
-               ibdy(jfld) = ib_bdy
-               igrid(jfld) = 3
-               ilen1(jfld) = nblen(igrid(jfld))
-               ilen3(jfld) = jpk
+                  dta%nread(2) = dta%nread(2) + 1
+               ENDIF
+
+            ENDIF
+
+            ! 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(lwp) write(numout,*) '++++++ reading in u3d field'
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_u3d
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 2
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  ilen3(jfld) = jpk
+                  IF( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1
+               ENDIF
+
+               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
+                  blf_i(jfld) = bn_v3d
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 3
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  ilen3(jfld) = jpk
+                  IF( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1
+               ENDIF
 
             ENDIF
 
             ! temperature and salinity
-            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1 ) THEN
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_tem
-               ibdy(jfld) = ib_bdy
-               igrid(jfld) = 1
-               ilen1(jfld) = nblen(igrid(jfld))
-               ilen3(jfld) = jpk
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_sal
-               ibdy(jfld) = ib_bdy
-               igrid(jfld) = 1
-               ilen1(jfld) = nblen(igrid(jfld))
-               ilen3(jfld) = jpk
+            IF( nn_tra_dta(ib_bdy) .eq. 1 ) THEN
+
+               IF( dta%ll_tem ) THEN
+                  if(lwp) write(numout,*) '++++++ reading in tem field'
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_tem
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  ilen3(jfld) = jpk
+               ENDIF
+
+               IF( dta%ll_sal ) THEN
+                  if(lwp) write(numout,*) '++++++ reading in sal field'
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_sal
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  ilen3(jfld) = jpk
+               ENDIF
 
             ENDIF
@@ -545,26 +598,32 @@
 #if defined key_lim2
             ! sea ice
-            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_frld
-               ibdy(jfld) = ib_bdy
-               igrid(jfld) = 1
-               ilen1(jfld) = nblen(igrid(jfld))
-               ilen3(jfld) = 1
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_hicif
-               ibdy(jfld) = ib_bdy
-               igrid(jfld) = 1
-               ilen1(jfld) = nblen(igrid(jfld))
-               ilen3(jfld) = 1
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_hsnif
-               ibdy(jfld) = ib_bdy
-               igrid(jfld) = 1
-               ilen1(jfld) = nblen(igrid(jfld))
-               ilen3(jfld) = 1
+            IF( nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN
+
+               IF( dta%ll_frld ) THEN
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_frld
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  ilen3(jfld) = 1
+               ENDIF
+
+               IF( dta%ll_hicif ) THEN
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_hicif
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  ilen3(jfld) = 1
+               ENDIF
+
+               IF( dta%ll_hsnif ) THEN
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_hsnif
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblen(igrid(jfld))
+                  ilen3(jfld) = 1
+               ENDIF
 
             ENDIF
@@ -581,4 +640,6 @@
             ENDIF
 
+            dta%nread(1) = nb_bdy_fld(ib_bdy)
+
          ENDIF ! nn_dta .eq. 1
       ENDDO ! ib_bdy
@@ -609,51 +670,81 @@
 
          nblen => idx_bdy(ib_bdy)%nblen
-         nblenrim => idx_bdy(ib_bdy)%nblenrim
-
-         IF (nn_dyn2d(ib_bdy) .gt. 0) THEN
-            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN
-               ilen0(1:3) = nblen(1:3)
-               ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) )
-               ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) )
-               IF ( nn_dyn2d(ib_bdy) .ne. jp_frs .and. (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) ) THEN
-                  jfld = jfld + 1
-                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1)
+         dta => dta_bdy(ib_bdy)
+
+         if(lwp) then
+            write(numout,*) '++++++ dta%ll_ssh = ',dta%ll_ssh
+            write(numout,*) '++++++ dta%ll_u2d = ',dta%ll_u2d
+            write(numout,*) '++++++ dta%ll_v2d = ',dta%ll_v2d
+            write(numout,*) '++++++ dta%ll_u3d = ',dta%ll_u3d
+            write(numout,*) '++++++ dta%ll_v3d = ',dta%ll_v3d
+            write(numout,*) '++++++ dta%ll_tem = ',dta%ll_tem
+            write(numout,*) '++++++ dta%ll_sal = ',dta%ll_sal
+         endif
+
+         IF ( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN
+            if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space'
+            IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) )
+            IF( dta%ll_u2d ) ALLOCATE( dta%u2d(nblen(2)) )
+            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( dta%ll_ssh ) THEN
+               if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow'
+               jfld = jfld + 1
+               dta%ssh => bf(jfld)%fnow(:,1,1)
+            ENDIF
+            IF ( dta%ll_u2d ) THEN
+               IF ( ln_full_vel_array(ib_bdy) ) THEN
+                  if(lwp) write(numout,*) '++++++ dta%u2d allocated space'
+                  ALLOCATE( dta%u2d(nblen(2)) )
                ELSE
-                  ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) )
-               ENDIF
-            ELSE
-               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN
-                  jfld = jfld + 1
-                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1)
-               ENDIF
+                  if(lwp) write(numout,*) '++++++ dta%u2d pointing to fnow'
+                  jfld = jfld + 1
+                  dta%u2d => bf(jfld)%fnow(:,1,1)
+               ENDIF
+            ENDIF
+            IF ( dta%ll_v2d ) THEN
+               IF ( ln_full_vel_array(ib_bdy) ) THEN
+                  if(lwp) write(numout,*) '++++++ dta%v2d allocated space'
+                  ALLOCATE( dta%v2d(nblen(3)) )
+               ELSE
+                  if(lwp) write(numout,*) '++++++ dta%v2d pointing to fnow'
+                  jfld = jfld + 1
+                  dta%v2d => bf(jfld)%fnow(:,1,1)
+               ENDIF
+            ENDIF
+         ENDIF
+
+         IF ( nn_dyn3d_dta(ib_bdy) .eq. 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 ( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN
+               if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%u2d => bf(jfld)%fnow(:,1,1)
+               dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:)
+            ENDIF
+            IF ( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN
+               if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%v2d => bf(jfld)%fnow(:,1,1)
-            ENDIF
-         ENDIF
-
-         IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN
-            ilen0(1:3) = nblen(1:3)
-            ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) )
-            ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) )
-         ENDIF
-         IF ( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ).or. &
-           &  ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.   &
-           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN
-            jfld = jfld + 1
-            dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:)
-            jfld = jfld + 1
-            dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:)
-         ENDIF
-
-         IF (nn_tra(ib_bdy) .gt. 0) THEN
-            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN
-               ilen0(1:3) = nblen(1:3)
-               ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) )
-               ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) )
-            ELSE
+               dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:)
+            ENDIF
+         ENDIF
+
+         IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN
+            if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space'
+            IF( dta%ll_tem ) ALLOCATE( dta_bdy(ib_bdy)%tem(nblen(1),jpk) )
+            IF( dta%ll_sal ) ALLOCATE( dta_bdy(ib_bdy)%sal(nblen(1),jpk) )
+         ELSE
+            IF( dta%ll_tem ) THEN
+               if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow'
                jfld = jfld + 1
                dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:)
+            ENDIF
+            IF( dta%ll_sal ) THEN 
+               if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow'
                jfld = jfld + 1
                dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:)
@@ -664,8 +755,7 @@
          IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN
             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN
-               ilen0(1:3) = nblen(1:3)
-               ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) )
-               ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) )
-               ALLOCATE( dta_bdy(ib_bdy)%hsnif(ilen0(1)) )
+               ALLOCATE( dta_bdy(ib_bdy)%frld(nblen(1)) )
+               ALLOCATE( dta_bdy(ib_bdy)%hicif(nblen(1)) )
+               ALLOCATE( dta_bdy(ib_bdy)%hsnif(nblen(1)) )
             ELSE
                jfld = jfld + 1
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 3991)
@@ -30,5 +30,4 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  !
-   USE domvvl          ! variable volume
 
    IMPLICIT NONE
@@ -57,6 +56,6 @@
       LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities
       !!
-      INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter
-      LOGICAL               :: ll_dyn2d, ll_dyn3d  
+      INTEGER               :: jk,ii,ij,ib_bdy,ib,igrd     ! Loop counter
+      LOGICAL               :: ll_dyn2d, ll_dyn3d, ll_orlanski
       !!
 
@@ -70,4 +69,10 @@
       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
+
       !-------------------------------------------------------
       ! Set pointers
@@ -77,5 +82,6 @@
       phur => hur
       phvr => hvr
-      CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 
+      CALL wrk_alloc(jpi,jpj,pua2d,pva2d) 
+      IF ( ll_orlanski ) CALL wrk_alloc(jpi,jpj,pub2d,pvb2d) 
 
       !-------------------------------------------------------
@@ -83,25 +89,35 @@
       !-------------------------------------------------------
 
-      pu2d(:,:) = 0.e0
-      pv2d(:,:) = 0.e0
-      IF (lk_vvl) THEN
+      ! "After" velocities: 
+
+      pua2d(:,:) = 0.e0
+      pva2d(:,:) = 0.e0
+      DO jk = 1, jpkm1   !! Vertically integrated momentum trends
+          pua2d(:,:) = pua2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+          pva2d(:,:) = pva2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+      END DO
+      pua2d(:,:) = pua2d(:,:) * phur(:,:)
+      pva2d(:,:) = pva2d(:,:) * phvr(:,:)
+      DO jk = 1 , jpkm1
+         ua(:,:,jk) = ua(:,:,jk) - pua2d(:,:)
+         va(:,:,jk) = va(:,:,jk) - pva2d(:,:)
+      END DO
+
+      ! "Before" velocities (required for Orlanski condition): 
+
+      IF ( ll_orlanski ) THEN
+         pub2d(:,:) = 0.e0
+         pvb2d(:,:) = 0.e0
          DO jk = 1, jpkm1   !! Vertically integrated momentum trends
-            pu2d(:,:) = pu2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-            pv2d(:,:) = pv2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+             pub2d(:,:) = pub2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ub(:,:,jk)
+             pvb2d(:,:) = pvb2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * vb(:,:,jk)
          END DO
-         pu2d(:,:) = pu2d(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) )
-         pv2d(:,:) = pv2d(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) )
-      ELSE
-         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
-            pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-            pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+         pub2d(:,:) = pub2d(:,:) * phur(:,:)
+         pvb2d(:,:) = pvb2d(:,:) * phvr(:,:)
+         DO jk = 1 , jpkm1
+            ub(:,:,jk) = ub(:,:,jk) - pub2d(:,:)
+            vb(:,:,jk) = vb(:,:,jk) - pvb2d(:,:)
          END DO
-         pu2d(:,:) = pu2d(:,:) * phur(:,:)
-         pv2d(:,:) = pv2d(:,:) * phvr(:,:)
-      ENDIF
-      DO jk = 1 , jpkm1
-         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) * umask(:,:,jk)
-         va(:,:,jk) = va(:,:,jk) - pv2d(:,:) * vmask(:,:,jk)
-      END DO
+      END IF
 
       !-------------------------------------------------------
@@ -119,9 +135,17 @@
 
       DO jk = 1 , jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) + pu2d(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) + pv2d(:,:) ) * vmask(:,:,jk)
+         ua(:,:,jk) = ( ua(:,:,jk) + pua2d(:,:) ) * umask(:,:,jk)
+         va(:,:,jk) = ( va(:,:,jk) + pva2d(:,:) ) * vmask(:,:,jk)
       END DO
 
-      CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 
+      IF ( ll_orlanski ) THEN
+         DO jk = 1 , jpkm1
+            ub(:,:,jk) = ( ub(:,:,jk) + pub2d(:,:) ) * umask(:,:,jk)
+            vb(:,:,jk) = ( vb(:,:,jk) + pvb2d(:,:) ) * vmask(:,:,jk)
+         END DO
+      END IF
+
+      CALL wrk_dealloc(jpi,jpj,pua2d,pva2d) 
+      IF ( ll_orlanski ) CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d) 
 
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn')
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 3991)
@@ -18,4 +18,5 @@
    USE dom_oce         ! ocean space and time domain
    USE bdy_oce         ! ocean open boundary conditions
+   USE bdylib          ! BDY library routines
    USE dynspg_oce      ! for barotropic variables
    USE phycst          ! physical constants
@@ -26,5 +27,5 @@
    PRIVATE
 
-   PUBLIC   bdy_dyn2d     ! routine called in dynspg_ts and bdy_dyn
+   PUBLIC   bdy_dyn2d   ! routine called in dynspg_ts and bdy_dyn
 
    !!----------------------------------------------------------------------
@@ -48,11 +49,15 @@
       DO ib_bdy=1, nb_bdy
 
-         SELECT CASE( nn_dyn2d(ib_bdy) )
-         CASE(jp_none)
+         SELECT CASE( cn_dyn2d(ib_bdy) )
+         CASE('none')
             CYCLE
-         CASE(jp_frs)
+         CASE('frs')
             CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
-         CASE(jp_flather)
+         CASE('flather')
             CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
+         CASE('orlanski')
+            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
+         CASE('orlanski_npo')
+            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
          CASE DEFAULT
             CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
@@ -89,5 +94,5 @@
          ij   = idx%nbj(jb,igrd)
          zwgt = idx%nbw(jb,igrd)
-         pu2d(ii,ij) = ( pu2d(ii,ij) + zwgt * ( dta%u2d(jb) - pu2d(ii,ij) ) ) * umask(ii,ij,1)
+         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1)
       END DO
       !
@@ -97,8 +102,8 @@
          ij   = idx%nbj(jb,igrd)
          zwgt = idx%nbw(jb,igrd)
-         pv2d(ii,ij) = ( pv2d(ii,ij) + zwgt * ( dta%v2d(jb) - pv2d(ii,ij) ) ) * vmask(ii,ij,1)
+         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1)
       END DO 
-      CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy ) 
-      CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
+      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy ) 
+      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated
       !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs')
@@ -133,4 +138,5 @@
       INTEGER  ::   jb, igrd                         ! dummy loop indices
       INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
+      REAL(wp), POINTER :: flagu, flagv              ! short cuts
       REAL(wp) ::   zcorr                            ! Flather correction
       REAL(wp) ::   zforc                            ! temporary scalar
@@ -160,10 +166,11 @@
          ii  = idx%nbi(jb,igrd)
          ij  = idx%nbj(jb,igrd) 
-         iim1 = ii + MAX( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice inside the boundary
-         iip1 = ii - MIN( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice outside the boundary 
+         flagu => idx%flagu(jb,igrd)
+         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary
+         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary 
          !
-         zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
+         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
          zforc = dta%u2d(jb)
-         pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
+         pua2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
       END DO
       !
@@ -173,17 +180,59 @@
          ii  = idx%nbi(jb,igrd)
          ij  = idx%nbj(jb,igrd) 
-         ijm1 = ij + MAX( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice inside the boundary
-         ijp1 = ij - MIN( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice outside the boundary 
+         flagv => idx%flagv(jb,igrd)
+         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary
+         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary 
          !
-         zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
+         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
          zforc = dta%v2d(jb)
-         pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
-      END DO
-      CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
-      CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy )   !
+         pva2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
+      END DO
+      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
+      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
       !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla')
       !
    END SUBROUTINE bdy_dyn2d_fla
+
+
+   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, ll_npo )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
+      !!             
+      !!              - Apply Orlanski radiation condition adaptively:
+      !!                  - radiation plus weak nudging at outflow points
+      !!                  - no radiation and strong nudging at inflow points
+      !! 
+      !!
+      !! 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
+      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
+      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
+
+      INTEGER  ::   ib, igrd                               ! dummy loop indices
+      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski')
+      !
+      igrd = 2      ! Orlanski bc on u-velocity; 
+      !            
+      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo )
+
+      igrd = 3      ! Orlanski bc on v-velocity
+      !  
+      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
+      !
+      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
+      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   !
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
+      !
+   END SUBROUTINE bdy_dyn2d_orlanski
+
 #else
    !!----------------------------------------------------------------------
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 3991)
@@ -19,4 +19,5 @@
    USE dom_oce         ! ocean space and time domain
    USE bdy_oce         ! ocean open boundary conditions
+   USE bdylib          ! for orlanski library routines
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  !
@@ -52,17 +53,17 @@
       DO ib_bdy=1, nb_bdy
 
-!!$         IF ( using Orlanski radiation conditions ) THEN 
-!!$            CALL bdy_rad( kt,  bdyidx(ib_bdy) )
-!!$         ENDIF
-
-         SELECT CASE( nn_dyn3d(ib_bdy) )
-         CASE(jp_none)
+         SELECT CASE( cn_dyn3d(ib_bdy) )
+         CASE('none')
             CYCLE
-         CASE(jp_frs)
+         CASE('frs')
             CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE(2)
+         CASE('specified')
             CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE(3)
+         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' )
@@ -109,5 +110,6 @@
          END DO
       END DO
-      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
+      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated  
+      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
       !
       IF( kt .eq. nit000 ) CLOSE( unit = 102 )
@@ -204,5 +206,6 @@
          END DO
       END DO 
-      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
+      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
+      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
       !
       IF( kt .eq. nit000 ) CLOSE( unit = 102 )
@@ -211,4 +214,42 @@
 
    END SUBROUTINE bdy_dyn3d_frs
+
+   SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, ll_npo )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  ***
+      !!             
+      !!              - Apply Orlanski radiation to baroclinic velocities. 
+      !!              - Wrapper routine for bdy_orlanski_3d
+      !! 
+      !!
+      !! 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
+      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
+      LOGICAL,                      INTENT(in) ::   ll_npo  ! switch for NPO version
+
+      INTEGER  ::   jb, igrd                               ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_orlanski')
+      !
+      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 
+      !
+      igrd = 2      ! Orlanski bc on u-velocity; 
+      !            
+      CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, ll_npo )
+
+      igrd = 3      ! Orlanski bc on v-velocity
+      !  
+      CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo )
+      !
+      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
+      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_orlanski')
+      !
+   END SUBROUTINE bdy_dyn3d_orlanski
+
 
    SUBROUTINE bdy_dyn3d_dmp( kt )
@@ -232,29 +273,29 @@
       ! Remove barotropic part from before velocity
       !-------------------------------------------------------
-      CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 
-
-      pu2d(:,:) = 0.e0
-      pv2d(:,:) = 0.e0
+      CALL wrk_alloc(jpi,jpj,pub2d,pvb2d) 
+
+      pub2d(:,:) = 0.e0
+      pvb2d(:,:) = 0.e0
 
       DO jk = 1, jpkm1
 #if defined key_vvl
-         pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk) 
-         pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk)
+         pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk) 
+         pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk)
 #else
-         pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk)
-         pv2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk)
+         pub2d(:,:) = pub2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk)
+         pvb2d(:,:) = pvb2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk)
 #endif
       END DO
 
       IF( lk_vvl ) THEN
-         pu2d(:,:) = pu2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
-         pv2d(:,:) = pv2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
+         pub2d(:,:) = pub2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
+         pvb2d(:,:) = pvb2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
       ELSE
-         pu2d(:,:) = pv2d(:,:) * hur(:,:)
-         pv2d(:,:) = pu2d(:,:) * hvr(:,:)
+         pub2d(:,:) = pvb2d(:,:) * hur(:,:)
+         pvb2d(:,:) = pub2d(:,:) * hvr(:,:)
       ENDIF
 
       DO ib_bdy=1, nb_bdy
-         IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN
+         IF ( ln_dyn3d_dmp(ib_bdy) .and. cn_dyn3d(ib_bdy) /= 'none' ) THEN
             igrd = 2                      ! Relaxation of zonal velocity
             DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
@@ -264,5 +305,5 @@
                DO jk = 1, jpkm1
                   ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
-                                   ub(ii,ij,jk) + pu2d(ii,ij)) ) * umask(ii,ij,jk)
+                                   ub(ii,ij,jk) + pub2d(ii,ij)) ) * umask(ii,ij,jk)
                END DO
             END DO
@@ -275,5 +316,5 @@
                DO jk = 1, jpkm1
                   va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
-                                   vb(ii,ij,jk) + pv2d(ii,ij)) ) * vmask(ii,ij,jk)
+                                   vb(ii,ij,jk) + pvb2d(ii,ij)) ) * vmask(ii,ij,jk)
                END DO
             END DO
@@ -281,5 +322,5 @@
       ENDDO
       !
-      CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 
+      CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d) 
       !
       CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim2.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim2.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim2.F90	(revision 3991)
@@ -50,8 +50,8 @@
       DO ib_bdy=1, nb_bdy
 
-         SELECT CASE( nn_ice_lim2(ib_bdy) )
-         CASE(jp_none)
+         SELECT CASE( cn_ice_lim2(ib_bdy) )
+         CASE('none')
             CYCLE
-         CASE(jp_frs)
+         CASE('frs')
             CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
          CASE DEFAULT
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3991)
@@ -79,6 +79,8 @@
       INTEGER  ::   jpbdtau, jpbdtas                       !   -       -
       INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
+      INTEGER  ::   i_offset, j_offset                     !   -       -
       INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
-      REAL   , POINTER  ::  flagu, flagv                   !    -   -
+      REAL(wp), POINTER  ::  flagu, flagv                  !    -   -
+      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields
       REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
       INTEGER, DIMENSION (2)                  ::   kdimsz
@@ -92,10 +94,10 @@
 
       !!
-      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             &
-         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
-         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &  
-         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,              &
+      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,                 &
+         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
+         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &  
+         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
 #if defined key_lim2
-         &             nn_ice_lim2, nn_ice_lim2_dta,                       &
+         &             cn_ice_lim2, nn_ice_lim2_dta,                           &
 #endif
          &             ln_vol, nn_volctl, nn_rimwidth
@@ -128,9 +130,9 @@
       ln_mask_file      = .false.
       cn_mask_file(:)   = ''
-      nn_dyn2d(:)       = 0
+      cn_dyn2d(:)       = ''
       nn_dyn2d_dta(:)   = -1  ! uninitialised flag
-      nn_dyn3d(:)       = 0
+      cn_dyn3d(:)       = ''
       nn_dyn3d_dta(:)   = -1  ! uninitialised flag
-      nn_tra(:)         = 0
+      cn_tra(:)         = ''
       nn_tra_dta(:)     = -1  ! uninitialised flag
       ln_tra_dmp(:)     = .false.
@@ -138,5 +140,5 @@
       rn_time_dmp(:)    = 1.
 #if defined key_lim2
-      nn_ice_lim2(:)    = 0
+      cn_ice_lim2(:)    = ''
       nn_ice_lim2_dta(:)= -1  ! uninitialised flag
 #endif
@@ -172,11 +174,33 @@
 
         IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
-        SELECT CASE( nn_dyn2d(ib_bdy) )                  
-          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
+        SELECT CASE( cn_dyn2d(ib_bdy) )                  
+          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')          
+             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')      
+             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')     
+             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') 
+             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) 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 DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
         END SELECT
-        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
+        IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
            SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   ! 
               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
@@ -193,12 +217,32 @@
 
         IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
-        SELECT CASE( nn_dyn3d(ib_bdy) )                  
-          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
-          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
+        SELECT CASE( cn_dyn3d(ib_bdy) )                  
+          CASE('none')
+             IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+             dta_bdy(ib_bdy)%ll_u3d = .false.
+             dta_bdy(ib_bdy)%ll_v3d = .false.
+          CASE('frs')       
+             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+             dta_bdy(ib_bdy)%ll_u3d = .true.
+             dta_bdy(ib_bdy)%ll_v3d = .true.
+          CASE('specified')
+             IF(lwp) WRITE(numout,*) '      Specified value'
+             dta_bdy(ib_bdy)%ll_u3d = .true.
+             dta_bdy(ib_bdy)%ll_v3d = .true.
+          CASE('zero')
+             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
+             dta_bdy(ib_bdy)%ll_u3d = .false.
+             dta_bdy(ib_bdy)%ll_v3d = .false.
+          CASE('orlanski')
+             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
+             dta_bdy(ib_bdy)%ll_u3d = .true.
+             dta_bdy(ib_bdy)%ll_v3d = .true.
+          CASE('orlanski_npo')
+             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
+             dta_bdy(ib_bdy)%ll_u3d = .true.
+             dta_bdy(ib_bdy)%ll_v3d = .true.
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
         END SELECT
-        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN
+        IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
            SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   ! 
               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
@@ -209,8 +253,8 @@
 
         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
-           IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN
+           IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
               ln_dyn3d_dmp(ib_bdy)=.false.
-           ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN
+           ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
               CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
@@ -218,4 +262,6 @@
               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
               IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
+              dta_bdy(ib_bdy)%ll_u3d = .true.
+              dta_bdy(ib_bdy)%ll_v3d = .true.
            ENDIF
         ELSE
@@ -225,13 +271,36 @@
 
         IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
-        SELECT CASE( nn_tra(ib_bdy) )                  
-          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
-          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions'
-          CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
+        SELECT CASE( cn_tra(ib_bdy) )                  
+          CASE('none')
+             IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+             dta_bdy(ib_bdy)%ll_tem = .false.
+             dta_bdy(ib_bdy)%ll_sal = .false.
+          CASE('frs')
+             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+             dta_bdy(ib_bdy)%ll_tem = .true.
+             dta_bdy(ib_bdy)%ll_sal = .true.
+          CASE('specified')
+             IF(lwp) WRITE(numout,*) '      Specified value'
+             dta_bdy(ib_bdy)%ll_tem = .true.
+             dta_bdy(ib_bdy)%ll_sal = .true.
+          CASE('neumann')
+             IF(lwp) WRITE(numout,*) '      Neumann conditions'
+             dta_bdy(ib_bdy)%ll_tem = .false.
+             dta_bdy(ib_bdy)%ll_sal = .false.
+          CASE('runoff')
+             IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
+             dta_bdy(ib_bdy)%ll_tem = .false.
+             dta_bdy(ib_bdy)%ll_sal = .false.
+          CASE('orlanski')
+             IF(lwp) WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
+             dta_bdy(ib_bdy)%ll_tem = .true.
+             dta_bdy(ib_bdy)%ll_sal = .true.
+          CASE('orlanski_npo')
+             IF(lwp) WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
+             dta_bdy(ib_bdy)%ll_tem = .true.
+             dta_bdy(ib_bdy)%ll_sal = .true.
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_tra' )
         END SELECT
-        IF( nn_tra(ib_bdy) .gt. 0 ) THEN
+        IF( cn_tra(ib_bdy) /= 'none' ) THEN
            SELECT CASE( nn_tra_dta(ib_bdy) )                   ! 
               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
@@ -242,13 +311,16 @@
 
         IF ( ln_tra_dmp(ib_bdy) ) THEN
-           IF ( nn_tra(ib_bdy).EQ.0 ) THEN
+           IF ( cn_tra(ib_bdy) == 'none' ) THEN
               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
               ln_tra_dmp(ib_bdy)=.false.
-           ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN
+           ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
               CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
+              IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
               IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
+              dta_bdy(ib_bdy)%ll_tem = .true.
+              dta_bdy(ib_bdy)%ll_sal = .true.
            ENDIF
         ELSE
@@ -259,10 +331,18 @@
 #if defined key_lim2
         IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
-        SELECT CASE( nn_ice_lim2(ib_bdy) )                  
-          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
+        SELECT CASE( cn_ice_lim2(ib_bdy) )                  
+          CASE('none')
+             IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+             dta_bdy(ib_bdy)%ll_frld  = .false.
+             dta_bdy(ib_bdy)%ll_hicif = .false.
+             dta_bdy(ib_bdy)%ll_hsnif = .false.
+          CASE('frs')
+             IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+             dta_bdy(ib_bdy)%ll_frld  = .true.
+             dta_bdy(ib_bdy)%ll_hicif = .true.
+             dta_bdy(ib_bdy)%ll_hsnif = .true.
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice_lim2' )
         END SELECT
-        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN 
+        IF( cn_ice_lim2(ib_bdy) /= 'none' ) THEN 
            SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   ! 
               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
@@ -744,8 +824,6 @@
                IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
                   IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
-                     CALL ctl_stop('bdy_init : ERROR : boundary data in file  &
-                                    must be defined in order of distance from edge nbr.', &
-                                   'A utility for re-ordering boundary coordinates and data &
-                                    files exists in the TOOLS/OBC directory')
+                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
+                                   'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory')
                   ENDIF    
                ENDIF
@@ -770,8 +848,9 @@
          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) )
-         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
+         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
@@ -941,4 +1020,5 @@
             ENDDO
          ENDDO 
+
          ! definition of the i- and j- direction local boundaries arrays
          ! used for sending the boudaries
@@ -994,4 +1074,6 @@
                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
+               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
             END DO
@@ -1096,4 +1178,13 @@
       ENDDO
 
+      ! bdyfmask required for flagu, flagv calculations below even though F-points 
+      ! not defined for BDY grid. 
+      DO ij = 2, jpjm1
+         DO ii = 2, jpim1
+            bdyfmask(ii,ij) = bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
+           &                * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
+         END DO      
+      END DO
+
       ! Lateral boundary conditions
       CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
@@ -1102,52 +1193,95 @@
       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.e0
+         idx_bdy(ib_bdy)%flagv(:,:) = 0.e0
          icount = 0 
 
-         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
-         !flagu =  0 : u is tangential
-         !flagu =  1 : u is normal to the boundary and is direction is inward
+         ! Calculate relationship of U direction to the local orientation of the boundary
+         ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
+         ! flagu =  0 : u is tangential
+         ! flagu =  1 : u is normal to the boundary and is direction is inward
   
-         igrd = 2      ! u-component 
-         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
-            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
-            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
-            zefl = bdytmask(nbi  ,nbj)
-            zwfl = bdytmask(nbi+1,nbj)
-            IF( zefl + zwfl == 2 ) THEN
-               icount = icount + 1
-            ELSE
-               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
-            ENDIF
+         DO igrd = 1,jpbgrd 
+            SELECT CASE( igrd )
+               CASE( 1 )
+                  cgrid = 'T' 
+                  pmask => umask(:,:,1)
+                  i_offset = 0
+               CASE( 2 ) 
+                  cgrid = 'U'
+                  pmask => bdytmask
+                  i_offset = 1
+               CASE( 3 ) 
+                  cgrid = 'V'
+                  pmask => fmask(:,:,1)
+                  i_offset = 0
+            END SELECT 
+            icount = 0
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
+               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+               zefl = pmask(nbi+i_offset-1,nbj)
+               zwfl = pmask(nbi+i_offset,nbj)
+               IF( zefl + zwfl == 2 * i_offset ) THEN
+                  icount = icount + 1
+                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
+               ELSE
+                  idx_bdy(ib_bdy)%flagu(ib,igrd)=-zefl+zwfl
+               ENDIF
+            END DO
+            IF( icount /= 0 ) THEN
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid,' grid points,',   &
+                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
+               IF(lwp) WRITE(numout,*) ' ========== '
+               IF(lwp) WRITE(numout,*)
+               nstop = nstop + 1
+            ENDIF 
          END DO
 
-         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
-         !flagv =  0 : u is tangential
-         !flagv =  1 : u is normal to the boundary and is direction is inward
-
-         igrd = 3      ! v-component
-         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
-            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
-            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
-            znfl = bdytmask(nbi,nbj  )
-            zsfl = bdytmask(nbi,nbj+1)
-            IF( znfl + zsfl == 2 ) THEN
-               icount = icount + 1
-            ELSE
-               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
-            END IF
+         ! Calculate relationship of V direction to the local orientation of the boundary
+         ! flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
+         ! flagv =  0 : u is tangential
+         ! flagv =  1 : u is normal to the boundary and is direction is inward
+
+         DO igrd = 1,jpbgrd 
+            SELECT CASE( igrd )
+               CASE( 1 )
+                  cgrid = 'T'
+                  pmask => vmask(:,:,1)
+                  j_offset = 0
+               CASE( 2 )
+                  cgrid = 'U'
+                  pmask => fmask(:,:,1)
+                  j_offset = 0
+               CASE( 3 )
+                  cgrid = 'V'
+                  pmask => bdytmask
+                  j_offset = 1
+            END SELECT 
+            icount = 0
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
+               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)
+               IF( znfl + zsfl == 2 * j_offset ) THEN
+                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
+                  icount = icount + 1
+               ELSE
+                  idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl
+               END IF
+            END DO
+            IF( icount /= 0 ) THEN
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid,' grid points,',   &
+                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
+               IF(lwp) WRITE(numout,*) ' ========== '
+               IF(lwp) WRITE(numout,*)
+               nstop = nstop + 1
+            ENDIF 
          END DO
 
-         IF( icount /= 0 ) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
-               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
-            IF(lwp) WRITE(numout,*) ' ========== '
-            IF(lwp) WRITE(numout,*)
-            nstop = nstop + 1
-         ENDIF 
-    
-      ENDDO
+      END DO
 
       ! Compute total lateral surface for volume correction:
@@ -1161,5 +1295,5 @@
                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
-               flagu => idx_bdy(ib_bdy)%flagu(ib)
+               flagu => idx_bdy(ib_bdy)%flagu(ib,igrd)
                bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
                   &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
@@ -1174,5 +1308,5 @@
                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
-               flagv => idx_bdy(ib_bdy)%flagv(ib)
+               flagv => idx_bdy(ib_bdy)%flagv(ib,igrd)
                bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
                   &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
@@ -1584,7 +1718,7 @@
       itest = 0
 
-      IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1
-      IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1
-      IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1
+      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
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 3991)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 3991)
@@ -0,0 +1,229 @@
+MODULE bdylib
+   !!======================================================================
+   !!                       ***  MODULE  bdylib  ***
+   !! Unstructured Open Boundary Cond. :  Library module of generic boundary algorithms.
+   !!======================================================================
+   !! History :  3.6  !  2013     (D. Storkey) new module
+   !!----------------------------------------------------------------------
+#if defined key_bdy 
+   !!----------------------------------------------------------------------
+   !!   'key_bdy' :                    Unstructured Open Boundary Condition
+   !!----------------------------------------------------------------------
+   !!   bdy_orlanski_2d
+   !!   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  !
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   bdy_orlanski_2d     ! routine called where?
+   PUBLIC   bdy_orlanski_3d     ! routine called where?
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, ll_npo )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_orlanski_2d  ***
+      !!             
+      !!              - Apply Orlanski radiation condition adaptively to 2D fields:
+      !!                  - radiation plus weak nudging at outflow points
+      !!                  - no radiation and strong nudging at inflow points
+      !! 
+      !!
+      !! 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
+
+      INTEGER  ::   jb                                     ! dummy loop indices
+      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
+      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
+      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
+      INTEGER  ::   flagu, flagv                           ! short cuts
+      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations
+      REAL(wp) ::   zout, zwgt, zdy_centred, zsign_ups
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_2d')
+
+      ! ----------------------------------!
+      ! Orlanski boundary conditions     :!
+      ! ----------------------------------! 
+     
+      !
+      DO jb = 1, idx%nblenrim(igrd)
+         ii  = idx%nbi(jb,igrd)
+         ij  = idx%nbj(jb,igrd) 
+         flagu = int( idx%flagu(jb,igrd) )
+         flagv = int( idx%flagv(jb,igrd) )
+         !
+         ! calculate positions of b-1 and b-2 points for this rim point
+         ! also (b-1,j-1) and (b-1,j+1) points
+         iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
+         ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv
+         !
+         iijm1 = ii - abs(flagv) ; iijp1 = ii + abs(flagv) 
+         ijjm1 = ij - abs(flagu) ; ijjp1 = ij + abs(flagu)
+         !
+         iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
+         ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
+         !
+         ! calculate normal (zcx) and tangential (zcy) components of radiation velocities:
+         zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1)
+         zdx = phia(iibm1,ijbm1) - phia(iibm2,ijbm2)
+         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1)
+         ! upstream differencing for tangential derivatives
+         zsign_ups = sign( 1., zdt * zdy_centred )
+         zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
+         zdy = zsign_ups        * ( phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1) ) &
+       &        + (1. - zsign_ups) * ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1   ) )
+         znor2 = zdx * zdx + zdy * zdy
+         zcx = zdt * zdx / znor2
+         zcy = zdt * zdy / znor2
+         !
+         ! update boundary value:
+         zout = sign( 1., zcx )
+         zout = 0.5*( zout + abs(zout) )
+         zwgt = (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd)
+         ! only apply radiation on outflow points 
+         if( ll_npo ) then     !! NPO version !!
+            phia(ii,ij) = (1.-zout) * phia(ii,ij)                                                   &
+           &            + zout      * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1) ) / ( 1. + zcx ) 
+         else                  !! full oblique radiation !!
+            zsign_ups = sign( 1., zcy )
+            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
+            phia(ii,ij) = (1.-zout) * phia(ii,ij)                                                   &
+           &            + zout      * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1)                         &
+           &                    - zsign_ups      * zcy * ( phib(ii   ,ij   ) - phib(iijm1,ijjm1 ) ) &
+           &                    - (1.-zsign_ups) * zcy * ( phib(iijp1,ijjp1) - phib(ii   ,ij    ) ) ) / ( 1. + zcx ) 
+         end if
+         phia(ii,ij) = phia(ii,ij) + zwgt * ( phi_ext(jb) - phia(ii,ij) ) 
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_2d')
+
+   END SUBROUTINE bdy_orlanski_2d
+
+
+   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, ll_npo )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_orlanski_3d  ***
+      !!             
+      !!              - Apply Orlanski radiation condition adaptively to 3D fields:
+      !!                  - radiation plus weak nudging at outflow points
+      !!                  - no radiation and strong nudging at inflow points
+      !! 
+      !!
+      !! 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
+
+      INTEGER  ::   jb, jk                                 ! dummy loop indices
+      INTEGER  ::   ii, ij, iibm1, iibm2, ijbm1, ijbm2     ! 2D addresses
+      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
+      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
+      INTEGER  ::   flagu, flagv                           ! short cuts
+      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations
+      REAL(wp) ::   zout, zwgt, zdy_centred, zsign_ups
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_orlanski_3d')
+
+      ! ----------------------------------!
+      ! Orlanski boundary conditions     :!
+      ! ----------------------------------! 
+     
+      DO jk = 1, jpk
+         !            
+         DO jb = 1, idx%nblenrim(igrd)
+            ii  = idx%nbi(jb,igrd)
+            ij  = idx%nbj(jb,igrd) 
+            flagu = int( idx%flagu(jb,igrd) )
+            flagv = int( idx%flagv(jb,igrd) )
+            !
+            ! calculate positions of b-1 and b-2 points for this rim point
+            ! also (b-1,j-1) and (b-1,j+1) points
+            iibm1 = ii + flagu ; iibm2 = ii + 2*flagu 
+            ijbm1 = ij + flagv ; ijbm2 = ij + 2*flagv
+            !
+            iijm1 = ii - abs(flagv) ; iijp1 = ii + abs(flagv) 
+            ijjm1 = ij - abs(flagu) ; ijjp1 = ij + abs(flagu)
+            !
+            iibm1jm1 = ii + flagu - abs(flagv) ; iibm1jp1 = ii + flagu + abs(flagv) 
+            ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
+            !
+            ! calculate normal (zcx) and tangential (zcy) components of radiation velocities:
+            zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk)
+            zdx = phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk)
+            zdy_centred = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1jm1,ijbm1jm1,jk)
+            ! upstream differencing for tangential derivatives
+            zsign_ups = sign( 1., zdt * zdy_centred )
+            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
+            zdy = zsign_ups        * ( phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) &
+           &    + (1. - zsign_ups) * ( phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk) )
+            znor2 = zdx * zdx + zdy * zdy
+            zcx = zdt * zdx / znor2
+            zcy = zdt * zdy / znor2
+            !
+            ! update boundary value:
+            zout = sign( 1., zcx )
+            zout = 0.5*( zout + abs(zout) )
+            zwgt = (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd)
+            ! only apply radiation on outflow points 
+            if( ll_npo ) then     !! NPO version !!
+               phia(ii,ij,jk) = (1.-zout) * phia(ii,ij,jk)                                                   &
+              &               + zout      * ( phib(ii,ij,jk) + zcx*phia(iibm1,ijbm1,jk) ) / ( 1. + zcx ) 
+            else                  !! full oblique radiation !!
+               zsign_ups = sign( 1., zcy )
+               zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) )
+               phia(ii,ij,jk) = (1.-zout) * phia(ii,ij,jk)                                                   &
+              &               + zout      * ( phib(ii,ij,jk) + zcx*phia(iibm1,ijbm1,jk)                      &
+              &                    - zsign_ups      * zcy * ( phib(ii   ,ij   ,jk) - phib(iijm1,ijjm1,jk ) ) &
+              &                    - (1.-zsign_ups) * zcy * ( phib(iijp1,ijjp1,jk) - phib(ii   ,ij   ,jk ) ) ) / ( 1. + zcx ) 
+            end if
+            phia(ii,ij,jk) = phia(ii,ij,jk) + zwgt * ( phi_ext(jb,jk) - phia(ii,ij,jk) ) 
+         END DO
+         !
+      END DO
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_orlanski_3d')
+
+   END SUBROUTINE bdy_orlanski_3d
+
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                   NO Unstruct Open Boundary Conditions
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext  )      ! Empty routine
+      WRITE(*,*) 'bdy_orlanski_2d: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_orlanski_2d
+   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext  )      ! Empty routine
+      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_orlanski_3d
+#endif
+
+   !!======================================================================
+END MODULE bdylib
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3991)
@@ -128,5 +128,5 @@
             ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
             ! relaxation area      
-            IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
+            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN
                ilen0(:)=nblen(:)
             ELSE
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 3991)
@@ -20,4 +20,5 @@
    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)
@@ -51,14 +52,18 @@
       DO ib_bdy=1, nb_bdy
 
-         SELECT CASE( nn_tra(ib_bdy) )
-         CASE(jp_none)
+         SELECT CASE( cn_tra(ib_bdy) )
+         CASE('none')
             CYCLE
-         CASE(jp_frs)
+         CASE('frs')
             CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
-         CASE(2)
+         CASE('specified')
             CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
-         CASE(3)
+         CASE('neumann')
             CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
-         CASE(4)
+         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
@@ -196,4 +201,38 @@
       !
    END SUBROUTINE bdy_tra_nmn
+ 
+
+   SUBROUTINE bdy_tra_orlanski( idx, dta, ll_npo )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_orlanski  ***
+      !!             
+      !!              - Apply Orlanski radiation to temperature and salinity. 
+      !!              - Wrapper routine for bdy_orlanski_3d
+      !! 
+      !!
+      !! 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
+
+      INTEGER  ::   igrd                                    ! grid index
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_orlanski')
+      !
+      igrd = 1      ! Orlanski bc on temperature; 
+      !            
+      CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_tem), tsa(:,:,:,jp_tem), dta%tem, ll_npo )
+
+      igrd = 1      ! Orlanski bc on salinity;
+      !  
+      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')
+      !
+
+   END SUBROUTINE bdy_tra_orlanski
+
 
    SUBROUTINE bdy_tra_rnf( idx, dta, kt )
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 3991)
@@ -104,5 +104,5 @@
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               zubtpecor = zubtpecor + idx%flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
+               zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
             END DO
          END DO
@@ -112,5 +112,5 @@
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               zubtpecor = zubtpecor + idx%flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
+               zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
             END DO
          END DO
@@ -136,6 +136,6 @@
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb) * zubtpecor * umask(ii,ij,jk)
-               ztranst = ztranst + idx%flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
+               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)
             END DO
          END DO
@@ -145,6 +145,6 @@
                ii = idx%nbi(jb,jgrd)
                ij = idx%nbj(jb,jgrd)
-               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb) * zubtpecor * vmask(ii,ij,jk)
-               ztranst = ztranst + idx%flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk)
+               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)
             END DO
          END DO
Index: /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3990)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3991)
@@ -594,6 +594,8 @@
          phur => hur_e
          phvr => hvr_e
-         pu2d => ua_e
-         pv2d => va_e
+         pua2d => ua_e
+         pva2d => va_e
+         pub2d => zub_e
+         pvb2d => zvb_e
 
          IF( lk_bdy )   CALL bdy_dyn2d( kt ) 
