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 4006)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4007)
@@ -109,5 +109,4 @@
    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
@@ -144,5 +143,5 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj), bdyfmask(jpi,jpj),     &  
+      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),     &  
          &      STAT=bdy_oce_alloc )
          !
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 4006)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 4007)
@@ -220,9 +220,9 @@
       igrd = 2      ! Orlanski bc on u-velocity; 
       !            
-      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, umask(:,:,1), ll_npo )
+      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, vmask(:,:,1), ll_npo )
+      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo )
       !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski')
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 4006)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 4007)
@@ -239,9 +239,9 @@
       igrd = 2      ! Orlanski bc on u-velocity; 
       !            
-      CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, umask, ll_npo )
+      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, vmask, ll_npo )
+      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
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 4006)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4007)
@@ -21,4 +21,5 @@
    !!   bdy_init       : Initialization of unstructured open boundaries
    !!----------------------------------------------------------------------
+   USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
    USE oce             ! ocean dynamics and tracers variables
@@ -92,4 +93,5 @@
       INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving
       INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates
+      REAL(wp), POINTER, DIMENSION(:,:)       ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat)
 
       !!
@@ -1178,14 +1180,16 @@
       ENDDO
 
-      ! bdyfmask required for flagu, flagv calculations below even though F-points 
-      ! not defined for BDY grid. 
+      ! For the flagu/flagv calculation below we require a version of fmask without
+      ! the land boundary condition (shlat) included:
+      CALL wrk_alloc(jpi,jpj,zfmask) 
       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)
+            zfmask(ii,ij) = tmask(ii,ij  ,1) * tmask(ii+1,ij  ,1)   &
+           &              * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1)
          END DO      
       END DO
 
       ! Lateral boundary conditions
+      CALL lbc_lnk( zfmask       , 'F', 1. )
       CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
       CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
@@ -1205,14 +1209,11 @@
             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)
+                  pmask => zfmask(:,:)
                   i_offset = 0
             END SELECT 
@@ -1223,14 +1224,15 @@
                zefl = pmask(nbi+i_offset-1,nbj)
                zwfl = pmask(nbi+i_offset,nbj)
-               IF( zefl + zwfl == 2 * i_offset ) THEN
+               ! This error check only works if you are using the bdyXmask arrays
+               IF( i_offset == 1 .and. zefl + zwfl == 2 ) 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
+                  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,',   &
+               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   &
                   ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
                IF(lwp) WRITE(numout,*) ' ========== '
@@ -1241,20 +1243,17 @@
 
          ! 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
+         ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
+         ! flagv =  0 : v is tangential
+         ! flagv =  1 : v 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)
+                  pmask => zfmask(:,:)
                   j_offset = 0
                CASE( 3 )
-                  cgrid = 'V'
                   pmask => bdytmask
                   j_offset = 1
@@ -1266,5 +1265,6 @@
                znfl = pmask(nbi,nbj+j_offset-1  )
                zsfl = pmask(nbi,nbj+j_offset)
-               IF( znfl + zsfl == 2 * j_offset ) THEN
+               ! This error check only works if you are using the bdyXmask arrays
+               IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN
                   IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj)
                   icount = icount + 1
@@ -1275,5 +1275,5 @@
             IF( icount /= 0 ) THEN
                IF(lwp) WRITE(numout,*)
-               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid,' grid points,',   &
+               IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   &
                   ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
                IF(lwp) WRITE(numout,*) ' ========== '
@@ -1324,4 +1324,6 @@
          DEALLOCATE(nbidta, nbjdta, nbrdta)
       ENDIF
+
+      CALL wrk_dealloc(jpi,jpj,zfmask) 
 
       IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
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 4006)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 4007)
@@ -34,5 +34,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, mask, ll_npo )
+   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, ll_npo )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_orlanski_2d  ***
@@ -50,5 +50,4 @@
       REAL(wp), DIMENSION(:,:),   INTENT(inout)  ::   phia     ! model after 2D field (to be updated)
       REAL(wp), DIMENSION(:),     INTENT(in)     ::   phi_ext  ! external forcing data
-      REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   mask     ! land/sea mask
       LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version
 
@@ -57,7 +56,12 @@
       INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
       INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
+      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices
       INTEGER  ::   flagu, flagv                           ! short cuts
       REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations
-      REAL(wp) ::   zout, zwgt, zdy_centred, zsign_ups
+      REAL(wp) ::   zout, zwgt, zdy_centred
+      REAL(wp) ::   zdy_left, zdy_right, zsign_ups
+      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask      ! land/sea mask for field
+      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_xdiv ! land/sea mask for x-derivatives
+      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_ydiv ! land/sea mask for y-derivatives
       !!----------------------------------------------------------------------
 
@@ -68,4 +72,25 @@
       ! ----------------------------------! 
      
+      SELECT CASE(igrd)
+         CASE(1)
+            pmask => tmask(:,:,1)
+            pmask_xdiv => umask(:,:,1)
+            ii_offset = 0
+            pmask_ydiv => vmask(:,:,1)
+            ij_offset = 0
+         CASE(2)
+            pmask => umask(:,:,1)
+            pmask_xdiv => tmask(:,:,1)
+            ii_offset = 1
+            pmask_ydiv => fmask(:,:,1)
+            ij_offset = 0
+         CASE(3)
+            pmask => vmask(:,:,1)
+            pmask_xdiv => fmask(:,:,1)
+            ii_offset = 0
+            pmask_ydiv => tmask(:,:,1)
+            ij_offset = 1
+         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' )
+      END SELECT
       !
       DO jb = 1, idx%nblenrim(igrd)
@@ -86,13 +111,24 @@
          ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
          !
-         ! calculate normal (zcx) and tangential (zcy) components of radiation velocities:
+         ! Calculate normal (zcx) and tangential (zcy) components of radiation velocities.
+         ! Mask derivatives to ensure correct land boundary conditions for each variable.
+         ! Centred derivative is calculated as average of "left" and "right" derivatives for 
+         ! this reason. 
          zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1)
-         zdx = phia(iibm1,ijbm1) - phia(iibm2,ijbm2)
-         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1)
+         zdx = ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) )                          &
+        &    * ( abs(iibm1-iibm2) * pmask_xdiv(iibm2+ii_offset,ijbm2          )   &
+        &      + abs(ijbm1-ijbm2) * pmask_ydiv(iibm2          ,ijbm2+ij_offset) ) 
+         zdy_left  = phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1)                  &
+        &    * ( (iibm1-iibm1jm1) * pmask_xdiv(iibm1jm1+ii_offset,ijbm1jm1          )   & 
+        &      + (ijbm1-ijbm1jm1) * pmask_ydiv(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
+         zdy_right = phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1)                     &
+        &    * ( (iibm1jp1-iibm1) * pmask_xdiv(iibm1+ii_offset,ijbm1)                   &
+        &      + (ijbm1jp1-ijbm1) * pmask_ydiv(iibm1          ,ijbm1+ij_offset) ) 
+         zdy_centred = 0.5 * ( zdy_left + zdy_right )
+!!$         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   ) )
+         zdy = zsign_ups * zdy_left + (1. - zsign_ups) * zdy_right
          znor2 = zdx * zdx + zdy * zdy
          znor2 = max(znor2,rsmall)
@@ -106,16 +142,16 @@
          ! only apply radiation on outflow points 
          if( ll_npo ) then     !! NPO version !!
-            phia(ii,ij) = (1.-zout) * phia(ii,ij)                                                   &
+            phia(ii,ij) = (1.-zout) * phib(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)                                                   &
+            phia(ii,ij) = (1.-zout) * phib(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) ) 
-         phia(ii,ij) = phia(ii,ij) * mask(ii,ij)
+!!$         phia(ii,ij) = phia(ii,ij) + zwgt * ( phi_ext(jb) - phib(ii,ij) ) 
+         phia(ii,ij) = phia(ii,ij) * pmask(ii,ij)
       END DO
       !
@@ -125,5 +161,5 @@
 
 
-   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, mask, ll_npo )
+   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, ll_npo )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_orlanski_3d  ***
@@ -141,5 +177,4 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated)
       REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   phi_ext  ! external forcing data
-      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   mask     ! land/sea mask
       LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version
 
@@ -148,7 +183,12 @@
       INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses
       INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses
+      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices
       INTEGER  ::   flagu, flagv                           ! short cuts
       REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations
-      REAL(wp) ::   zout, zwgt, zdy_centred, zsign_ups
+      REAL(wp) ::   zout, zwgt, zdy_centred
+      REAL(wp) ::   zdy_left, zdy_right,  zsign_ups
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_xdiv ! land/sea mask for x-derivatives
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_ydiv ! land/sea mask for y-derivatives
       !!----------------------------------------------------------------------
 
@@ -159,4 +199,26 @@
       ! ----------------------------------! 
      
+      SELECT CASE(igrd)
+         CASE(1)
+            pmask => tmask(:,:,:)
+            pmask_xdiv => umask(:,:,:)
+            ii_offset = 0
+            pmask_ydiv => vmask(:,:,:)
+            ij_offset = 0
+         CASE(2)
+            pmask => umask(:,:,:)
+            pmask_xdiv => tmask(:,:,:)
+            ii_offset = 1
+            pmask_ydiv => fmask(:,:,:)
+            ij_offset = 0
+         CASE(3)
+            pmask => vmask(:,:,:)
+            pmask_xdiv => fmask(:,:,:)
+            ii_offset = 0
+            pmask_ydiv => tmask(:,:,:)
+            ij_offset = 1
+         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' )
+      END SELECT
+
       DO jk = 1, jpk
          !            
@@ -178,13 +240,24 @@
             ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu) 
             !
-            ! calculate normal (zcx) and tangential (zcy) components of radiation velocities:
+            ! Calculate normal (zcx) and tangential (zcy) components of radiation velocities.
+            ! Mask derivatives to ensure correct land boundary conditions for each variable.
+            ! Centred derivative is calculated as average of "left" and "right" derivatives for 
+            ! this reason. 
             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)
+            zdx = phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk)                     &
+        &    * ( (iibm1-iibm2) * pmask_xdiv(iibm2+ii_offset,ijbm2          ,jk)   &
+        &      + (ijbm1-ijbm2) * pmask_ydiv(iibm2          ,ijbm2+ij_offset,jk) ) 
+            zdy_left  = phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk)            &
+        &    * ( (iibm1-iibm1jm1) * pmask_xdiv(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   & 
+        &      + (ijbm1-ijbm1jm1) * pmask_ydiv(iibm1jm1          ,ijbm1jm1+ij_offset,jk) ) 
+            zdy_right = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk)      &
+        &    * ( (iibm1jp1-iibm1) * pmask_xdiv(iibm1+ii_offset,ijbm1          ,jk)   &
+        &      + (ijbm1jp1-ijbm1) * pmask_ydiv(iibm1          ,ijbm1+ij_offset,jk) ) 
+            zdy_centred = 0.5 * ( zdy_left + zdy_right )
+!!$            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) )
+            zdy = zsign_ups * zdy_left + (1. - zsign_ups) * zdy_right
             znor2 = zdx * zdx + zdy * zdy
             znor2 = max(znor2,rsmall)
@@ -198,16 +271,16 @@
             ! only apply radiation on outflow points 
             if( ll_npo ) then     !! NPO version !!
-               phia(ii,ij,jk) = (1.-zout) * phia(ii,ij,jk)                                                   &
+               phia(ii,ij,jk) = (1.-zout) * phib(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)                                                   &
+               phia(ii,ij,jk) = (1.-zout) * phib(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) ) 
-            phia(ii,ij,jk) = phia(ii,ij,jk) * mask(ii,ij,jk)
+!!$            phia(ii,ij,jk) = phia(ii,ij,jk) + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) 
+            phia(ii,ij,jk) = phia(ii,ij,jk) * pmask(ii,ij,jk)
          END DO
          !
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 4006)
+++ /branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 4007)
@@ -224,9 +224,9 @@
       igrd = 1      ! Orlanski bc on temperature; 
       !            
-      CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_tem), tsa(:,:,:,jp_tem), dta%tem, tmask, ll_npo )
+      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, tmask, ll_npo )
+      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')
