Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 4291)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 4292)
@@ -6,4 +6,5 @@
    !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite
    !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!----------------------------------------------------------------------
 #if defined key_bdy 
@@ -11,6 +12,9 @@
    !!   'key_bdy' :                    Unstructured Open Boundary Condition
    !!----------------------------------------------------------------------
-   !!   bdy_dyn2d      : Apply open boundary conditions to barotropic variables.
-   !!   bdy_dyn2d_fla    : Apply Flather condition
+   !!   bdy_dyn2d          : Apply open boundary conditions to barotropic variables.
+   !!   bdy_dyn2d_frs      : Apply Flow Relaxation Scheme 
+   !!   bdy_dyn2d_fla      : Apply Flather condition
+   !!   bdy_dyn2d_orlanski : Orlanski Radiation
+   !!   bdy_ssh            : Duplicate sea level across open boundaries
    !!----------------------------------------------------------------------
    USE timing          ! Timing
@@ -18,4 +22,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 +31,6 @@
    PRIVATE
 
-   PUBLIC   bdy_dyn2d     ! routine called in dynspg_ts and bdy_dyn
+   PUBLIC   bdy_dyn2d   ! routine called in dynspg_ts and bdy_dyn
+   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv
 
    !!----------------------------------------------------------------------
@@ -48,11 +54,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 +99,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 +107,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,9 +143,13 @@
       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
+      REAL(wp) ::   zflag, z1_2                      !    "        "
       !!----------------------------------------------------------------------
 
       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla')
+
+      z1_2 = 0.5_wp
 
       ! ---------------------------------!
@@ -160,10 +174,15 @@
          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) )
-         zforc = dta%u2d(jb)
-         pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
+         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
+
+         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
+         ! Use characteristics method instead
+         zflag = ABS(flagu)
+         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
+         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
       END DO
       !
@@ -173,17 +192,110 @@
          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) )
-         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 )   !
+         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
+
+         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
+         ! Use characteristics method instead
+         zflag = ABS(flagv)
+         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
+         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * 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
+
+   SUBROUTINE bdy_ssh( zssh )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_ssh  ***
+      !!
+      !! ** Purpose : Duplicate sea level across open boundaries
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level
+      !!
+      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers
+      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       "
+
+      igrd = 1                       ! Everything is at T-points here
+
+      DO ib_bdy = 1, nb_bdy
+         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+            ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+            ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+            ! Set gradient direction:
+            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
+            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
+            IF ( zcoef1+zcoef2 == 0 ) THEN
+               ! corner
+!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1)
+!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + &
+!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + &
+!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + &
+!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1)
+               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1)
+               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + &
+                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + &
+                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + &
+                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1)
+               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1)
+            ELSE
+               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
+               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1)
+            ENDIF
+         END DO
+
+         ! Boundary points should be updated
+         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy )
+      END DO
+
+   END SUBROUTINE bdy_ssh
+
 #else
    !!----------------------------------------------------------------------
@@ -192,8 +304,11 @@
 CONTAINS
    SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
-      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt
+      INTEGER, intent(in) :: kt
+      WRITE(*,*) 'bdy_dyn2d: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_dyn2d
+
 #endif
 
    !!======================================================================
 END MODULE bdydyn2d
+
