Index: /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90	(revision 12690)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DIA/diaar5.F90	(revision 12691)
@@ -136,6 +136,5 @@
          !                                         ! steric sea surface height
          CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
-         !CW: commenting out the following line, as seems weird to zero the array before output & only way to get rhop output
-         !zrhop(:,:,jpk) = 0._wp
+         zrhop(:,:,jpk) = 0._wp
          CALL iom_put( 'rhop', zrhop )
          !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90	(revision 12690)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90	(revision 12691)
@@ -17,4 +17,5 @@
    USE ldfslp         ! iso-neutral slopes 
    USE zdf_oce        ! ocean vertical physics
+   USE phycst
    !
    USE in_out_manager ! I/O manager
@@ -211,6 +212,7 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       !
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulapdz, zvlapdz   ! -1*bhm * d/dz(del^2 u) at u- and v-point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap     ! laplacian at u- and v-point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulapdz, zvlapdz ! -1*bhm * d/dz(del^2 u) at u- and v-point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmu              ! = bhm / avm
       !!----------------------------------------------------------------------
       !
@@ -222,4 +224,17 @@
       !
   
+      ! Calculate ratio of bhm and avm for stabilising correction terms
+      ! and add to avm to be included in the vertical diffusion calculation later.
+      DO jk = 2, jpkm1 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, jpim1   ! vector opt.
+               zmu(ji,jj,jk) = ( bhm(ji,jj,jk) / (avm(ji,jj,jk) + rsmall) ) * wmask(ji,jj,jk)
+               avm(ji,jj,jk) = avm(ji,jj,jk) + zmu(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+      CALL lbc_lnk_multi( 'dyn_ldf_bgm', zmu, 'W', 1. , avm, 'W', 1. ) 
+      
+      ! Calculate (del2 u)
       CALL dyn_ldf_lap_no_ahm( kt, pub, pvb, zulap, zvlap, 1 )   
 
@@ -230,10 +245,11 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, jpim1   ! vector opt.
-                zulapdz(ji,jj,jk) = -0.5_wp*(bhm(ji+1,jj  ,jk)+bhm(ji,jj,jk))*(zulap(ji,jj,jk-1) - zulap(ji,jj,jk)) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk)
-                zvlapdz(ji,jj,jk) = -0.5_wp*(bhm(ji  ,jj+1,jk)+bhm(ji,jj,jk))*(zvlap(ji,jj,jk-1) - zvlap(ji,jj,jk)) * wvmask(ji,jj,jk) / e3vw_n(ji,jj,jk)
+                zulapdz(ji,jj,jk) = ( -0.5_wp*(bhm(ji+1,jj  ,jk)+bhm(ji,jj,jk))*(zulap(ji,jj,jk-1) - zulap(ji,jj,jk)) &
+     &                                -0.5_wp*(zmu(ji+1,jj  ,jk)+zmu(ji,jj,jk))*(pub  (ji,jj,jk-1) - pub  (ji,jj,jk)) ) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk) 
+                zvlapdz(ji,jj,jk) = ( -0.5_wp*(bhm(ji  ,jj+1,jk)+bhm(ji,jj,jk))*(zvlap(ji,jj,jk-1) - zvlap(ji,jj,jk)) &
+     &                                -0.5_wp*(zmu(ji  ,jj+1,jk)+zmu(ji,jj,jk))*(pvb  (ji,jj,jk-1) - pvb  (ji,jj,jk)) ) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk) 
             ENDDO
          ENDDO
       ENDDO
-
 
 !CW: set boundary conditions: d/dz(del^2 u) = 0 at top and bottom, so that eddy-induced velocity, w*=0
@@ -252,7 +268,4 @@
         ENDDO
 
-
-
-
 !!    calculate d/dz(-bhm * d/dz(del^2 u))
                                                        ! ===============
Index: /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90	(revision 12690)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/LDF/ldfdyn.F90	(revision 12691)
@@ -521,4 +521,5 @@
             ELSEIF( ln_dynldf_bgm ) THEN                 !bi-Laplacian GM operator (mask only)
                bhm(:,:,1:jpkm1) =       bhm(:,:,1:jpkm1)   * wmask(:,:,1:jpkm1)
+               CALL lbc_lnk( 'ldf_dyn_init', bhm , 'W', 1. ) 
                !
             ENDIF
