Index: /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90	(revision 14282)
+++ /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90	(revision 14283)
@@ -715,5 +715,5 @@
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
-          IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(jj,ji) == 1 ) lpyc(ji,jj) = .FALSE. 
+          IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE. 
         END DO
       END DO        
@@ -915,9 +915,9 @@
              IF ( j_ddh(ji,jj) == 0 ) THEN
 ! Place holding code. Parametrization needs checking for these conditions.
-               zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird
+               zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird
                zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
                zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
              ELSE
-               zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird
+               zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird
                zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
                zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
@@ -928,5 +928,6 @@
              DO jk = imld(ji,jj), ibld(ji,jj)
                 zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
-                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
+                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse * &
+                     & ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
              END DO
              zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) )
@@ -935,5 +936,6 @@
              DO jk = imld(ji,jj), ibld(ji,jj)
                zznd_pyc = -( gdepw_n(ji,jj,jk) -zhbl(ji,jj) ) / zdh(ji,jj)
-               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)**2 ) * zvw_bse *  &
+                    & ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
              END DO
            ENDIF  ! lpyc
@@ -1067,5 +1069,5 @@
             IF ( .not. lconv(ji,jj) ) THEN
                DO jk = 2, ibld(ji,jj)
-                  znd = ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about
+                  znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about
                   IF ( znd >= 0.0 ) THEN
                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
@@ -1086,5 +1088,4 @@
              IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
                 DO jk= 2, ibld(ji,jj)
-                   znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
                    ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
@@ -1407,6 +1408,6 @@
                  zbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln )
                ELSE
-                 zbeta_d_sc = 1.0
-                 zbeta_v_sc = 1.0
+                 zbeta_d_sc(ji,jj) = 1.0
+                 zbeta_v_sc(ji,jj) = 1.0
                ENDIF
              ELSE
@@ -1448,5 +1449,5 @@
                    zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj)
                    zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj)  )  / MAX(zvispyc_n_sc(ji,jj), 1.e-8)
-                   zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic )
+                   zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic )
                    zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
                    DO jk = imld(ji,jj) , ibld(ji,jj)
@@ -1604,5 +1605,5 @@
 
      DO jj = 2, jpjm1
-       DO ji = 2, jpjm1
+       DO ji = 2, jpim1
          IF ( lshear(ji,jj) ) THEN
            IF ( lconv(ji,jj) ) THEN
@@ -1779,4 +1780,7 @@
                 zpe_mle_layer = 0._wp
                 zpe_mle_ref = 0._wp
+                zthermal = rab_n(ji,jj,1,jp_tem)
+                zbeta    = rab_n(ji,jj,1,jp_sal)
+
                 DO jk = ibld(ji,jj), mld_prof(ji,jj)
                   zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) )
@@ -1816,5 +1820,5 @@
                      lpyc(ji,jj) = .FALSE.
                    ELSE
-                      lpyc = .TRUE.
+                      lpyc(ji,jj) = .TRUE.
                    ENDIF
                  ELSE
@@ -2503,4 +2507,8 @@
       END DO
 
+      ! ensure salinity > 0 in unset values so EOS doesn't give FP error with fpe0 on
+      ztsm_midu(:,jpj,jp_sal) = 10.
+      ztsm_midv(jpi,:,jp_sal) = 10.
+
       CALL eos_rab(ztsm_midu, zmld_midu, zabu)
       CALL eos_rab(ztsm_midv, zmld_midv, zabv)
