Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdy_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdy_oce.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdy_oce.F90	(revision 11083)
@@ -85,5 +85,5 @@
    !
    INTEGER                    ::   nb_bdy                   !: number of open boundary sets
-   INTEGER                    ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run)
+   INTEGER, DIMENSION(jp_bdy) ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run)
    INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme
    INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P 
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydta.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydta.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydta.F90	(revision 11083)
@@ -243,5 +243,5 @@
                         IF( ln_full_vel_array(jbdy) ) THEN
                            CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
-                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy,   &
+                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy(jbdy),   &
                                      & fvl=ln_full_vel_array(jbdy)  )
                         ELSE
@@ -313,5 +313,5 @@
                      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, jpk_bdy=nb_jpk_bdy,   &
+                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy(jbdy),   &
                                   & fvl=ln_full_vel_array(jbdy) )
                   ENDIF
@@ -446,5 +446,5 @@
       NAMELIST/nambdy_dta/ bn_a_i, bn_h_i, bn_h_s
 #endif
-      NAMELIST/nambdy_dta/ ln_full_vel, nb_jpk_bdy
+      NAMELIST/nambdy_dta/ ln_full_vel
       !!---------------------------------------------------------------------------
       !
@@ -508,9 +508,9 @@
       ! Read namelists
       ! --------------
-      REWIND(numnam_ref)
       REWIND(numnam_cfg)
       jfld = 0 
       DO jbdy = 1, nb_bdy         
          IF( nn_dta(jbdy) == 1 ) THEN
+            REWIND(numnam_ref)
             READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
 901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydyn2d.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydyn2d.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdydyn2d.F90	(revision 11083)
@@ -187,5 +187,5 @@
          ! Use characteristics method instead
          zflag = ABS(flagu)
-         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij)
+         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(ii+NINT(flagu),ij)
          pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
       END DO
@@ -205,5 +205,5 @@
          ! Use characteristics method instead
          zflag = ABS(flagv)
-         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1)
+         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ij+NINT(flagv))
          pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
       END DO
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyice.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyice.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyice.F90	(revision 11083)
@@ -57,6 +57,7 @@
       INTEGER ::   jbdy   ! BDY set index
       !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('bdy_ice_thd')
+      ! controls
+      IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing
+      IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
       CALL ice_var_glo2eqv
@@ -78,6 +79,8 @@
       CALL ice_var_agg(1)
       !
-      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
-      IF( ln_timing )   CALL timing_stop('bdy_ice_thd')
+      ! controls
+      IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints
+      IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing
       !
    END SUBROUTINE bdy_ice
@@ -148,8 +151,8 @@
             jpbound = 0   ;   ib = ji   ;   jb = jj
             !
-            IF( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji+1 ; jb = jj
-            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji-1 ; jb = jj
-            IF( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. )   jpbound = 1 ; ib = ji   ; jb = jj+1
-            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. )   jpbound = 1 ; ib = ji   ; jb = jj-1
+            IF( u_ice(ji  ,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. )   jpbound = 1 ; ib = ji+1
+            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji  ,jj  ,1) == 0. )   jpbound = 1 ; ib = ji-1
+            IF( v_ice(ji  ,jj  ) < 0. .AND. vmask(ji  ,jj-1,1) == 0. )   jpbound = 1 ; jb = jj+1
+            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj  ,1) == 0. )   jpbound = 1 ; jb = jj-1
             !
             IF( nn_ice_dta(jbdy) == 0 )   jpbound = 0 ; ib = ji ; jb = jj   ! case ice boundaries = initial conditions
@@ -306,5 +309,5 @@
                      ! one of the two zmsk is always 0 (because of zflag)
                      zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
-                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
+                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) )   ! 0 if no ice
                      !  
                      ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice)
@@ -329,5 +332,5 @@
                      ! one of the two zmsk is always 0 (because of zflag)
                      zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
-                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
+                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) )   ! 0 if no ice
                      !  
                      ! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice)
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyini.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyini.F90	(revision 11083)
@@ -140,5 +140,5 @@
       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
       CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid
-      INTEGER :: com_east, com_west, com_south, com_north          ! Flags for boundaries sending
+      INTEGER :: com_east, com_west, com_south, com_north, jpk_max  ! Flags for boundaries sending
       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
@@ -397,5 +397,5 @@
           IF(lwp) WRITE(numout,*)
         ENDIF
-        IF( nb_jpk_bdy > 0 ) THEN
+        IF( nb_jpk_bdy(ib_bdy) > 0 ) THEN
            IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'
         ELSE
@@ -516,25 +516,16 @@
          ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
             &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
-
-         IF( nb_jpk_bdy>0 ) THEN
-            ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) )
-            ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) )
-            ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) )
-         ELSE
-            ALLOCATE( dta_global(jpbdtau, 1, jpk) )
-            ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO
-            ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO
-         ENDIF
+         
+         jpk_max = MAXVAL(nb_jpk_bdy)
+         jpk_max = MAX(jpk_max, jpk)
+
+         ALLOCATE( dta_global(jpbdtau, 1, jpk_max) )
+         ALLOCATE( dta_global_z(jpbdtau, 1, jpk_max) ) ! needed ?? TODO
+         ALLOCATE( dta_global_dz(jpbdtau, 1, jpk_max) )! needed ?? TODO
 
          IF ( icount>0 ) THEN
-            IF( nb_jpk_bdy>0 ) THEN
-               ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) )
-               ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) )
-               ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) )
-            ELSE
-               ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
-               ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO
-               ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO  
-            ENDIF
+            ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk_max) )
+            ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk_max) ) ! needed ?? TODO
+            ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk_max) )! needed ?? TODO  
          ENDIF
          ! 
@@ -960,5 +951,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
-                       if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
+                       if( ii == (nlcit(nowe+1)-1) ) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -974,5 +965,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
-                       if((com_east_b .ne. 1) .and. (ii == 2)) then
+                       if( ii == 2 ) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -989,5 +980,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
-                       if((com_west_b .ne. 1) .and. (ii == (nlcit(nowe+1)-1))) then
+                       if( ii == (nlcit(nowe+1)-1) ) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
@@ -1004,5 +995,5 @@
                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
-                       if((com_east_b .ne. 1) .and. (ii == 2)) then
+                       if( ii == 2 ) then
                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
                           if((ij == 2) .and. (nbondj == 0 .or. nbondj == 1)) then
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DIA/diacfl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DIA/diacfl.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DIA/diacfl.F90	(revision 11083)
@@ -29,9 +29,4 @@
    REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number 
 
-!!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc !
-!!gm          I don't understand why.
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
-!!gm end
-
    PUBLIC   dia_cfl       ! routine called by step.F90
    PUBLIC   dia_cfl_init  ! routine called by nemogcm
@@ -55,8 +50,8 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER                ::   ji, jj, jk                            ! dummy loop indices
-      REAL(wp)               ::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars
-      INTEGER , DIMENSION(3) ::   iloc_u , iloc_v , iloc_w , iloc       ! workspace
-!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
+      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices
+      REAL(wp)                         ::   z2dt, zCu_max, zCv_max, zCw_max  ! local scalars
+      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
       !!----------------------------------------------------------------------
       !
@@ -71,5 +66,5 @@
       DO jk = 1, jpk       ! calculate Courant numbers
          DO jj = 1, jpj
-            DO ji = 1, fs_jpim1   ! vector opt.
+            DO ji = 1, jpi
                zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction
                zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction
@@ -111,5 +106,5 @@
       !                    ! write out to file
       IF( lwp ) THEN
-         WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
+         WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
          WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
          WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
@@ -172,8 +167,4 @@
       rCw_max = 0._wp
       !
-!!gm required to work
-      ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) )
-!!gm end
-      !      
    END SUBROUTINE dia_cfl_init
 
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/dynkeg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/dynkeg.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/dynkeg.F90	(revision 11083)
@@ -74,10 +74,10 @@
       INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
       !
-      INTEGER  ::   ji, jj, jk, jb    ! dummy loop indices
-      INTEGER  ::   ii, ifu, ib_bdy   ! local integers
-      INTEGER  ::   ij, ifv, igrd     !   -       -
-      REAL(wp) ::   zu, zv            ! local scalars
+      INTEGER  ::   ji, jj, jk, jb           ! dummy loop indices
+      INTEGER  ::   ifu, ifv, igrd, ib_bdy   ! local integers
+      REAL(wp) ::   zu, zv                   ! local scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
+      REAL(wp)  :: zweightu, zweightv
       !!----------------------------------------------------------------------
       !
@@ -97,32 +97,4 @@
       
       zhke(:,:,jpk) = 0._wp
-      
-      IF (ln_bdy) THEN
-         ! 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)
-                     ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
-                     un(ii-ifu,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)
-                     ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
-                     vn(ii,ij-ifv,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  ==!
@@ -140,4 +112,28 @@
             END DO
          END DO
+         !
+         IF (ln_bdy) THEN
+            ! 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 = 1           ! compensating null velocity on the bdy
+                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1
+                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1
+                     DO jk = 1, jpkm1
+                        zhke(ji,jj,jk) = 0._wp
+                        zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk)
+                        zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)
+                        zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
+                        zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
+                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu) 
+                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv) 
+                     END DO
+                  END DO
+               END IF
+               CALL lbc_bdy_lnk( 'dynkeg', zhke, 'T', 1., ib_bdy )   ! send 2 and recv jpi, jpj used in the computation of the speed tendencies
+            END DO
+         END IF
          !
       CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
@@ -158,14 +154,37 @@
             END DO
          END DO
+         IF (ln_bdy) THEN
+            ! 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 = 1           ! compensation null velocity on land at the bdy
+                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1
+                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1
+                     DO jk = 1, jpkm1
+                        zhke(ji,jj,jk) = 0._wp
+                        zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) &
+                             &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) ) 
+                        zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) &
+                             &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) )
+                        zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
+                           &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
+                           &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
+                           &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
+                        zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
+                           &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
+                           &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
+                           &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
+                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu )
+                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv )
+                     END DO
+                  END DO
+               END IF
+            END DO
+         END IF
          CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
          !
-      END SELECT
-
-      IF (ln_bdy) THEN
-         ! restore velocity masks at points outside boundary
-         un(:,:,:) = un(:,:,:) * umask(:,:,:)
-         vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
-      ENDIF      
-
+      END SELECT 
       !
       DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/sshwzv.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/DYN/sshwzv.F90	(revision 11083)
@@ -297,7 +297,6 @@
          IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
          IF(lwp) WRITE(numout,*) '~~~~~ '
-         !
-         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all)
-      ENDIF
+      ENDIF
+      !
       !
       DO jk = 1, jpkm1            ! calculate Courant numbers
@@ -305,27 +304,27 @@
             DO ji = 2, fs_jpim1   ! vector opt.
                z1_e3w = 1._wp / e3w_n(ji,jj,jk)
-               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &
-                  &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
-                  &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
-                  &                        * r1_e1e2t(ji,jj)                                                  &
-                  &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
-                  &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
-                  &                        * r1_e1e2t(ji,jj)                                                  &
-                  &                      ) * z1_e3w
+               Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &  ! 2*rdt and not r2dt (for restartability)
+                  &                             + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
+                  &                                 MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
+                  &                               * r1_e1e2t(ji,jj)                                                  &
+                  &                             + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
+                  &                                 MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
+                  &                               * r1_e1e2t(ji,jj)                                                  &
+                  &                             ) * z1_e3w
             END DO
          END DO
       END DO
+      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
       !
       CALL iom_put("Courant",Cu_adv)
       !
-      wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero
       IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
          DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
-            DO jj = 2, jpjm1                            ! w where necessary
-               DO ji = 2, fs_jpim1   ! vector opt.
+            DO jj = 1, jpj                              ! w where necessary
+               DO ji = 1, jpi
                   !
                   zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
                   !
-                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit
+                  IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
                      zcff = 0._wp
                   ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
@@ -346,5 +345,6 @@
       ELSE
          ! Fully explicit everywhere
-         Cu_adv = 0.0_wp                                ! Reuse array to output coefficient
+         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient
+         wi    (:,:,:) = 0._wp
       ENDIF
       CALL iom_put("wimp",wi) 
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/LBC/lib_mpp.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/LBC/lib_mpp.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/LBC/lib_mpp.F90	(revision 11083)
@@ -1480,6 +1480,7 @@
       LOGICAL         , OPTIONAL, INTENT(in   ) ::   ld_lbc, ld_glb, ld_dlg
       !!
+      CHARACTER(len=128)                        ::   ccountname  ! name of a subroutine to count communications
       LOGICAL ::   ll_lbc, ll_glb, ll_dlg
-      INTEGER ::    ji,  jj,  jk,  jh, jf   ! dummy loop indices
+      INTEGER ::    ji,  jj,  jk,  jh, jf, jcount   ! dummy loop indices
       !!----------------------------------------------------------------------
       !
@@ -1538,13 +1539,21 @@
          WRITE(numcom,*) ' '
          WRITE(numcom,*) ' lbc_lnk called'
-         jj = 1
-         DO ji = 2, n_sequence_lbc
-            IF( crname_lbc(ji-1) /= crname_lbc(ji) ) THEN
-               WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_lbc(ji-1))
-               jj = 0
+         DO ji = 1, n_sequence_lbc - 1
+            IF ( crname_lbc(ji) /= 'already counted' ) THEN
+               ccountname = crname_lbc(ji)
+               crname_lbc(ji) = 'already counted'
+               jcount = 1
+               DO jj = ji + 1, n_sequence_lbc
+                  IF ( ccountname ==  crname_lbc(jj) ) THEN
+                     jcount = jcount + 1
+                     crname_lbc(jj) = 'already counted'
+                  END IF
+               END DO
+               WRITE(numcom,'(A, I4, A, A)') ' - ', jcount,' times by subroutine ', TRIM(ccountname)
             END IF
-            jj = jj + 1 
          END DO
-         WRITE(numcom,'(A, I4, A, A)') ' - ', jj,' times by subroutine ', TRIM(crname_lbc(n_sequence_lbc))
+         IF ( crname_lbc(n_sequence_lbc) /= 'already counted' ) THEN
+            WRITE(numcom,'(A, I4, A, A)') ' - ', 1,' times by subroutine ', TRIM(crname_lbc(ncom_rec_max))
+         END IF
          WRITE(numcom,*) ' '
          IF ( n_sequence_glb > 0 ) THEN
Index: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/ZDF/zdfphy.F90	(revision 10901)
+++ NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/ZDF/zdfphy.F90	(revision 11083)
@@ -132,6 +132,7 @@
       IF( ln_zad_Aimp ) THEN
          IF( zdf_phy_alloc() /= 0 )   &
-        &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
-         wi(:,:,:) = 0._wp
+            &       CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' )
+         Cu_adv(:,:,:) = 0._wp
+         wi    (:,:,:) = 0._wp
       ENDIF
       !                          !==  Background eddy viscosity and diffusivity  ==!
