Index: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 6140)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 7646)
@@ -5,8 +5,6 @@
    !!======================================================================
    !! History :  3.6  !  2013     (D. Storkey) original code
+   !!            4.0  !  2014     (T. Lovato) Generalize OBC structure
    !!----------------------------------------------------------------------
-#if defined key_bdy 
-   !!----------------------------------------------------------------------
-   !!   'key_bdy' :                    Unstructured Open Boundary Condition
    !!----------------------------------------------------------------------
    !!   bdy_orlanski_2d
@@ -25,13 +23,105 @@
    PRIVATE
 
-   PUBLIC   bdy_orlanski_2d     ! routine called where?
-   PUBLIC   bdy_orlanski_3d     ! routine called where?
+   PUBLIC   bdy_frs, bdy_spe, bdy_nmn, bdy_orl
+   PUBLIC   bdy_orlanski_2d
+   PUBLIC   bdy_orlanski_3d
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
    !! $Id$ 
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
+
+   SUBROUTINE bdy_frs( idx, pta, dta )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_frs  ***
+      !!
+      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries.
+      !!
+      !! Reference : Engedahl H., 1995, Tellus, 365-382.
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
+      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      !!
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_frs')
+      ! 
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblen(igrd)
+         DO ik = 1, jpkm1
+            ii = idx%nbi(ib,igrd) 
+            ij = idx%nbj(ib,igrd)
+            zwgt = idx%nbw(ib,igrd)
+            pta(ii,ij,ik) = ( pta(ii,ij,ik) + zwgt * (dta(ib,ik) - pta(ii,ij,ik) ) ) * tmask(ii,ij,ik)
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_frs')
+      !
+   END SUBROUTINE bdy_frs
+
+   SUBROUTINE bdy_spe( idx, pta, dta )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_spe  ***
+      !!
+      !! ** Purpose : Apply a specified value for tracers at open boundaries.
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
+      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      !!
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_spe')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            pta(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik)
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_spe')
+      !
+   END SUBROUTINE bdy_spe
+
+   SUBROUTINE bdy_orl( idx, ptb, pta, dta, ll_npo )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_orl  ***
+      !!
+      !! ** Purpose : Apply Orlanski radiation for tracers at open boundaries.
+      !!              This is a wrapper routine for bdy_orlanski_3d below
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
+      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptb  ! before tracer field
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      LOGICAL,                             INTENT(in) ::   ll_npo  ! switch for NPO version
+      !!
+      INTEGER  ::   igrd                                    ! grid index
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_orl')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      !
+      CALL bdy_orlanski_3d( idx, igrd, ptb(:,:,:), pta(:,:,:), dta, ll_npo )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_orl')
+      !
+   END SUBROUTINE bdy_orl
 
    SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, ll_npo )
@@ -355,17 +445,76 @@
    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
+   SUBROUTINE bdy_nmn( idx, igrd, phia )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_nmn  ***
+      !!                    
+      !! ** Purpose : Duplicate the value at open boundaries, zero gradient.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,                    INTENT(in)     ::   igrd     ! grid index
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated)
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      !! 
+      REAL(wp) ::   zcoef, zcoef1, zcoef2
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field
+      REAL(wp), POINTER, DIMENSION(:,:)        :: bdypmask      ! land/sea mask for field
+      INTEGER  ::   ib, ik   ! dummy loop indices
+      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses
+      !!----------------------------------------------------------------------
+      !!
+      IF( nn_timing == 1 ) CALL timing_start('bdy_nmn')
+      !
+      SELECT CASE(igrd)
+         CASE(1)
+            pmask => tmask(:,:,:)
+            bdypmask => bdytmask(:,:)
+         CASE(2)
+            pmask => umask(:,:,:)
+            bdypmask => bdyumask(:,:)
+         CASE(3)
+            pmask => vmask(:,:,:)
+            bdypmask => bdyvmask(:,:)
+         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' )
+      END SELECT
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ! search the sense of the gradient
+            zcoef1 = bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik) +  bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik)
+            zcoef2 = bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik) +  bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik)
+            IF ( nint(zcoef1+zcoef2) == 0) THEN
+               ! corner **** we probably only want to set the tangentail component for the dynamics here
+               zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) +  pmask(ii,ij-1,ik) +  pmask(ii,ij+1,ik)
+               IF (zcoef > .5_wp) THEN ! Only set none isolated points.
+                 phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik) + &
+                   &              phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik) + &
+                   &              phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik) + &
+                   &              phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)
+                 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik)
+               ELSE
+                 phia(ii,ij,ik) = phia(ii,ij  ,ik) * pmask(ii,ij  ,ik)
+               ENDIF
+            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN
+               ! oblique corner **** we probably only want to set the normal component for the dynamics here
+               zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij  ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij  ) + &
+                   &   pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) +  pmask(ii,ij+1,ik)*bdypmask(ii,ij+1  )
+               phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik)*bdypmask(ii-1,ij  ) + &
+                   &            phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik)*bdypmask(ii+1,ij  )  + &
+                   &            phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik)*bdypmask(ii,ij -1 ) + &
+                   &            phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)*bdypmask(ii,ij+1  )
+ 
+               phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik)
+            ELSE
+               ip = nint(bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik))
+               jp = nint(bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik))
+               phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik)
+            ENDIF
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn')
+      !
+   END SUBROUTINE bdy_nmn
 
    !!======================================================================
