Index: /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
===================================================================
--- /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 6990)
+++ /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 6991)
@@ -59,6 +59,10 @@
          CASE('specified')
             CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+         CASE('zerograd')
+            CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
          CASE('zero')
             CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+         CASE('neumann')
+            CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy )
          CASE('orlanski')
             CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
@@ -117,4 +121,52 @@
 
    END SUBROUTINE bdy_dyn3d_spe
+
+   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***
+      !!
+      !! ** Purpose : - Enforce a zero gradient of normal velocity
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                     ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
+      !!
+      INTEGER  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   fu, fv
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad')
+      !
+      igrd = 2                      ! Copying tangential velocity into bdy points
+      DO jb = 1, idx%nblenrim(igrd)
+         DO jk = 1, jpkm1
+            ii   = idx%nbi(jb,igrd)
+            ij   = idx%nbj(jb,igrd)
+            fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 )
+            ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) &
+                        &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu )
+         END DO
+      END DO
+      !
+      igrd = 3                      ! Copying tangential velocity into bdy points
+      DO jb = 1, idx%nblenrim(igrd)
+         DO jk = 1, jpkm1
+            ii   = idx%nbi(jb,igrd)
+            ij   = idx%nbj(jb,igrd)
+            fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 )
+            va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) &
+                        &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv )
+         END DO
+      END DO
+      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 )
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad')
+   END SUBROUTINE bdy_dyn3d_zgrad
 
    SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
@@ -303,4 +355,35 @@
    END SUBROUTINE bdy_dyn3d_dmp
 
+   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***
+      !!             
+      !!              - Apply Neumann condition to baroclinic velocities.
+      !!              - Wrapper routine for bdy_nmn
+      !!
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
+      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
+      INTEGER  ::   jb, igrd                               ! dummy loop indices
+      !!----------------------------------------------------------------------
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn')
+      !
+      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.
+      !
+      igrd = 2      ! Neumann bc on u-velocity;
+      !           
+      CALL bdy_nmn( idx, igrd, ua )
+      igrd = 3      ! Neumann bc on v-velocity
+      ! 
+      CALL bdy_nmn( idx, igrd, va )
+      !
+      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_nmn')
+      !
+   END SUBROUTINE bdy_dyn3d_nmn
+   
 #else
    !!----------------------------------------------------------------------
Index: /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 6990)
+++ /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 6991)
@@ -213,4 +213,12 @@
              dta_bdy(ib_bdy)%ll_u3d = .true.
              dta_bdy(ib_bdy)%ll_v3d = .true.
+          CASE('neumann')             
+             IF(lwp) WRITE(numout,*) '      Neumann conditions'
+             dta_bdy(ib_bdy)%ll_u3d = .false.
+             dta_bdy(ib_bdy)%ll_v3d = .false.
+          CASE('zerograd')
+             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities'
+             dta_bdy(ib_bdy)%ll_u3d = .false.
+             dta_bdy(ib_bdy)%ll_v3d = .false.
           CASE('zero')
              IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
Index: /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
===================================================================
--- /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 6990)
+++ /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 6991)
@@ -26,4 +26,5 @@
    PUBLIC   bdy_orlanski_2d     ! routine called where?
    PUBLIC   bdy_orlanski_3d     ! routine called where?
+   PUBLIC   bdy_nmn     ! routine called where?
 
    !!----------------------------------------------------------------------
@@ -354,4 +355,76 @@
    END SUBROUTINE bdy_orlanski_3d
 
+   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) * pmask(ii,ij,ik)
+            ENDIF
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn')
+      !
+   END SUBROUTINE bdy_nmn
 
 #else
@@ -366,4 +439,7 @@
       WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_orlanski_3d
+   SUBROUTINE bdy_nmn( idx, igrd, phia )      ! Empty routine
+      WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_nmn
 #endif
 
Index: /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 6990)
+++ /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 6991)
@@ -160,6 +160,7 @@
       !! 
       REAL(wp) ::   zwgt           ! boundary weight
-      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
-      INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses
+      REAL(wp) ::   zcoef, zcoef1,zcoef2
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses
       !!----------------------------------------------------------------------
       !
@@ -174,5 +175,5 @@
             zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
             zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
-            IF ( zcoef1+zcoef2 == 0) THEN
+            IF ( NINT(zcoef1+zcoef2) == 0) THEN
                ! corner
                zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
@@ -181,15 +182,15 @@
                  &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + &
                  &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik)
-               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik)
                tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + &
                  &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + &
                  &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + &
                  &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik)
-               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik)
             ELSE
-               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
-               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
-               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik)
-               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik)
+               ip = NINT(bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ))
+               jp = NINT(bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1))
+               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) * tmask(ii,ij,ik)
             ENDIF
          END DO
Index: /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6990)
+++ /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6991)
@@ -189,5 +189,5 @@
                zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
                       &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
-               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
+               z2d(ji,jj) = 0.5_wp * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
                !
             ENDDO
Index: /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 6990)
+++ /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 6991)
@@ -24,5 +24,8 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
-
+#if defined key_bdy
+   USE bdy_oce         ! ocean open boundary conditions
+#endif
+   
    IMPLICIT NONE
    PRIVATE
@@ -78,4 +81,9 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
+#if defined key_bdy
+      INTEGER  ::   jb                 ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers
+      INTEGER  ::   fu, fv
+#endif
       !!----------------------------------------------------------------------
       !
@@ -97,4 +105,32 @@
       
       zhke(:,:,jpk) = 0._wp
+
+#if defined key_bdy
+      ! Maria Luneva & Fred Wobus: July-2016
+      ! compensate for lack of turbulent kinetic energy on liquid bdy points
+      DO ib_bdy = 1, nb_bdy
+         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
+            igrd = 2           ! Copying normal velocity into points outside bdy
+            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+               DO jk = 1, jpkm1
+                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+                  fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
+                  un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
+               END DO
+            END DO
+            !
+            igrd = 3           ! Copying normal velocity into points outside bdy
+            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+               DO jk = 1, jpkm1
+                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+                  fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
+                  vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
+               END DO
+            END DO
+         ENDIF
+      ENDDO 
+#endif 
       
       SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
@@ -134,4 +170,11 @@
       END SELECT
       !
+
+#if defined key_bdy
+      ! restore velocity masks at points outside boundary
+      un(:,:,:) = un(:,:,:) * umask(:,:,:)
+      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
+#endif     
+
       DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
          DO jj = 2, jpjm1
Index: anches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/TRD/trdtrc.F90
===================================================================
--- /branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/TRD/trdtrc.F90	(revision 6990)
+++ 	(revision )
@@ -1,23 +1,0 @@
-MODULE trdtrc
-   !!======================================================================
-   !!                       ***  MODULE trdtrc  ***
-   !!  Dummy module
-   !!======================================================================
-   !!----------------------------------------------------------------------
-   !!   Dummy module                                             NO TOP use
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE trd_trc( ptrtrd, kjn, ktrd, kt )
-      INTEGER ::   kt, kjn, ktrd   
-      REAL    ::   ptrtrd(:,:,:)  
-      WRITE(*,*) 'trd_trc : You should not have seen this print! error?', ptrtrd(1,1,1)
-      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn, ktrd, kt
-   END SUBROUTINE trd_trc
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!======================================================================
-END MODULE trdtrc
