Index: NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/DIA/diawri.F90	(revision 12324)
+++ NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/DIA/diawri.F90	(revision 12386)
@@ -43,4 +43,5 @@
    USE zdfdrg         ! ocean vertical physics: top/bottom friction
    USE zdfmxl         ! mixed layer
+   USE zdfosm         ! mixed layer
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -926,4 +927,18 @@
          CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
       ENDIF
+
+      IF( ln_zdfosm ) THEN
+         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)   )    ! now boundary-layer depth
+         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)    )    ! now mixed-layer depth
+         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask       )    ! w-level diffusion
+         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask       )    ! now w-level viscosity
+         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask       )    ! non-local t forcing
+         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask       )    ! non-local s forcing
+         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*wmask       )    ! non-local u forcing
+         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamu*wmask       )    ! non-local v forcing
+         IF( ln_osm_mle ) THEN
+            CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1)   )    ! now transition-layer depth
+         END IF
+      ENDIF
  
 #if defined key_si3
Index: NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/TRA/tramle.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/TRA/tramle.F90	(revision 12324)
+++ NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/TRA/tramle.F90	(revision 12386)
@@ -254,6 +254,6 @@
     DO jj = 1, jpjm1
        DO ji = 1, fs_jpim1   ! vector opt.
-          zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
-          zhv(ji,jj) = 1._wp / zhv(ji,jj)
+          zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu
+          zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 
        END DO
     END DO
Index: NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/ZDF/zdfosm.F90	(revision 12324)
+++ NEMO/branches/UKMO/NEMO_4.0.1_FKOSM/src/OCE/ZDF/zdfosm.F90	(revision 12386)
@@ -1358,5 +1358,5 @@
      REAL(wp) :: ztgrad, zsgrad, zbgrad
      REAL(wp) :: zgamma_b_nd, zgamma_c, znd
-     REAL(wp) :: zzeta_s=0.3
+     REAL(wp) :: zzeta_s=0.3, ztmp
 
      DO jj = 2, jpjm1
@@ -1367,18 +1367,19 @@
                          & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh &
                          &  .OR.  zdb_bl(ji,jj) > 0._wp)) THEN  ! zdhdt could be <0 due to FK, hence check zdhdt>0
-                       ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj)
-                       zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj)
-                       zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj)
-                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
+                       ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
+                       ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj)
+                       zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj)
+                       zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj)
+                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
                        zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /&
                             & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check
                        DO jk = 2, ibld(ji,jj)+ibld_ext
-                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
+                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp
                           IF ( znd <= zzeta_s ) THEN
-                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * &
+                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * &
                                   & EXP( -6.0 * ( znd -zzeta_s )**2 )
-                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * &
+                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * &
                                   & EXP( -6.0 * ( znd -zzeta_s )**2 )
-                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * &
+                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * &
                                   & EXP( -6.0 * ( znd -zzeta_s )**2 )
                           ELSE
@@ -1395,9 +1396,10 @@
                     IF ( zdb_bl(ji,jj) > 0._wp ) THEN
                        IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
-                          ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj)
-                          zsgrad = zds_bl(ji,jj) / zhbl(ji,jj)
-                          zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
+                          ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
+                          ztgrad = zdt_bl(ji,jj) * ztmp
+                          zsgrad = zds_bl(ji,jj) * ztmp
+                          zbgrad = zdb_bl(ji,jj) * ztmp
                           DO jk = 2, ibld(ji,jj)
-                             znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
+                             znd = gdepw_n(ji,jj,jk) * ztmp
                              zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
                              zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
@@ -1405,9 +1407,10 @@
                           END DO
                        ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
-                          ztgrad = zdt_bl(ji,jj) / zdh(ji,jj)
-                          zsgrad = zds_bl(ji,jj) / zdh(ji,jj)
-                          zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
+                          ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
+                          ztgrad = zdt_bl(ji,jj) * ztmp
+                          zsgrad = zds_bl(ji,jj) * ztmp
+                          zbgrad = zdb_bl(ji,jj) * ztmp
                           DO jk = 2, ibld(ji,jj)
-                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
+                             znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp
                              zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
                              zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
