Index: NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90	(revision 14699)
+++ NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90	(revision 14704)
@@ -858,4 +858,6 @@
     ! Rate of change of hbl
     CALL zdf_osm_calculate_dhdt( zdhdt )
+    ! Test if surface boundary layer coupled to bottom.
+    lcoup(:,:) = .FALSE.                                    ! ag 19/03
     DO jj = 2, jpjm1
        DO ji = 2, jpim1
@@ -863,13 +865,15 @@
           ! adjustment to represent limiting by ocean bottom
           IF ( mbkt(ji,jj) >2 ) THEN ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access
-             IF( zhbl_t(ji,jj) >= gdepw_n(ji, jj, mbkt(ji,jj) - 2 ) ) THEN
-                zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) - 2) - depth_tol)! ht_n(:,:))
+             IF( zhbl_t(ji,jj) > gdepw_n(ji, jj, mbkt(ji,jj) - 2 ) ) THEN
+                zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) - 2))! ht_n(:,:))
                 lpyc(ji,jj) = .FALSE.
+                lcoup(ji,jj) = .TRUE.                            ! ag 19/03
              ENDIF
           ENDIF
 #ifdef key_osm_debug
           IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-             WRITE(narea+100,'(2(a,g11.3),/,2(a,g11.3))')'after zdf_osm_calculate_dhdt: zhbl_t=',zhbl_t(ji,jj), 'hbl=', hbl(ji,jj),&
-                  & 'delta hbl from dzdhdt', zdhdt(ji,jj)*rn_rdt,' delta hbl from w ', wn(ji,jj,ibld(ji,jj))*rn_rdt
+             WRITE(narea+100,'(2(a,g11.3),/,2(a,g11.3),2(a,l7))')'after zdf_osm_calculate_dhdt: zhbl_t=',zhbl_t(ji,jj), 'hbl=', hbl(ji,jj),&
+                  & 'delta hbl from dzdhdt', zdhdt(ji,jj)*rn_rdt,' delta hbl from w ', wn(ji,jj,ibld(ji,jj))*rn_rdt, &
+                  & ' lcoup= ',lcoup(ji,jj), ' lpyc= ', lpyc(ji,jj)
              FLUSH(narea+100)
           END IF
@@ -890,18 +894,4 @@
        END DO
     END DO
-
-! Test if surface boundary layer coupled to bottom.
-
-    lcoup(:,:) = .FALSE.                                    ! ag 19/03
-    DO jj = 2, jpjm1                                        ! ag 19/03
-      DO ji = 2, jpim1                                      ! ag 19/03
-        IF ( mbkt(ji,jj) >2 ) THEN ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access
-           IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,mbkt(ji,jj) - 2) ) THEN ! ag 19/03
-              zhbl_t(ji,jj) = gdepw_n(ji,jj,mbkt(ji,jj)-2)     ! ag 19/03
-              lcoup(ji,jj) = .TRUE.                            ! ag 19/03
-           ENDIF                                               ! ag 19/03
-        ENDIF                                               ! ag 19/03
-      END DO                                                ! ag 19/03
-    END DO                                                  ! ag 19/03
 
     !
@@ -932,7 +922,7 @@
 #ifdef key_osm_debug
                 IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                   WRITE(narea+100,'(a)')'After setting pycnocline thickness BL running aground: lpyc= Fl: ibld(ji,jj) >= mbkt(ji,jj) -2'
+                   WRITE(narea+100,'(a)')'After setting pycnocline thickness BL running aground: lpyc= F5: ibld(ji,jj) >= mbkt(ji,jj) -2'
                    WRITE(narea+100,'(2(a,i7),2(a,g11.3))')' ibld=',ibld(ji,jj),' imld=',imld(ji,jj), ' zdh=',zdh(ji,jj), ' zhml=',zhml(ji,jj)
-                   WRITE(narea+100,'(2(a,g11.3))')'dh=',dh,' hml=',hml(ji,jj)
+                   WRITE(narea+100,'(2(a,g11.3))')'dh=',dh(ji,jj),' hml=',hml(ji,jj)
                    FLUSH(narea+100)
                 END IF
@@ -1559,6 +1549,6 @@
           DO ji = 2, jpim1
              DO jk = ibld(ji,jj) + 1, jpkm1
-                zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
-                zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
+                zdiffut(ji,jj,jk) = MAX(zdiffut(ji,jj,jk), zrimix(ji,jj,jk)*rn_difri)
+                zviscos(ji,jj,jk) = MAX(zviscos(ji,jj,jk), zrimix(ji,jj,jk)*rn_difri)
              END DO
           END DO
@@ -1572,5 +1562,5 @@
           DO ji = 2, jpim1
              DO jk = ibld(ji,jj) + 1, jpkm1
-                IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
+                IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = MAX(rn_difconv,zdiffut(ji,jj,jk))
              END DO
           END DO
@@ -1870,11 +1860,11 @@
                   pvispyc_s_sc(ji,jj) = 0._wp  ! ag 19/03
                   IF(lcoup(ji,jj) ) THEN   ! ag 19/03
-                    ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = Cd|ub|
+                    ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub|
                     !     already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90
                     !  Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub)
                     ! wet-cell averaging ..
-                    zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
-                    zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
-                    pb_coup(ji,jj) = 0.4 * SQRT(rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
+                    zmsku = 0.5 * ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
+                    zmskv = 0.5 * ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
+                    pb_coup(ji,jj) = 0.4 * SQRT(-rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
                   &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  ) )
 
@@ -1883,10 +1873,11 @@
 #ifdef key_osm_debug
                     IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                       WRITE(narea+100,'(3(a,g11.3))')' lcoup = T; 1st pz_b= ', pz_b, 'pb_coup', pb_coup(ji,jj), ' pc_coup_vis', pc_coup_vis(ji,jj)
-                       FLUSH(narea+100)
+                        WRITE(narea+100,'(4(a,g11.3))')' lcoup = T; 1st pz_b= ', pz_b, ' pb_coup ', pb_coup(ji,jj), ' pc_coup_vis ', pc_coup_vis(ji,jj), ' rCdU_bot ',rCdU_bot(ji,jj)
+                        WRITE(narea+100,'(2(a,g11.3))')' zmsku ', zmsku, ' zmskv ', zmskv
+                      FLUSH(narea+100)
                     END IF
 #endif
 #ifdef key_osm_debug
-                    WRITE(narea+400,'(4(a,i7))') ' lcoup = T at ji=',ji,' jj= ',jj,' jig= ', mig(ji), 'jjg= ', mjg(jj)
+                    WRITE(narea+400,'(4(a,i7))') ' lcoup = T at ji=',ji,' jj= ',jj,' jig= ', mig(ji), ' jjg= ', mjg(jj)
                     WRITE(narea+400,'(3(a,g11.3))') '1st pz_b= ', pz_b, 'pb_coup', pb_coup(ji,jj), ' pc_coup_vis', pc_coup_vis(ji,jj)
                     FLUSH(narea+400)
@@ -1903,5 +1894,5 @@
 #endif
 #ifdef key_osm_debug
-                    WRITE(narea+400,'(3(a,g11.3))') '2nd pz_b= ', pz_b,' pc_coup_vis', pc_coup_vis(ji,jj)
+                    WRITE(narea+400,'(3(a,g11.3))') '2nd pz_b= ', pz_b,' pc_coup_dif', pc_coup_dif(ji,jj)
                     FLUSH(narea+400)
 #endif
