diff --git a/model/src/wav_import_export.F90 b/model/src/wav_import_export.F90
index 64c577ecec..1880ad2024 100644
--- a/model/src/wav_import_export.F90
+++ b/model/src/wav_import_export.F90
@@ -846,7 +846,7 @@ subroutine export_fields (gcomp, rc)
         call init_get_isea(isea, jsea)
         ix  = mapsf(isea,1)
         iy  = mapsf(isea,2)
-        if (mapsta(iy,ix) == 1) then
+        if (mapsta(iy,ix) /= 0) then
           sw_tauox(jsea) = tauox(jsea)
           sw_tauoy(jsea) = tauoy(jsea)
         endif
@@ -883,7 +883,7 @@ subroutine export_fields (gcomp, rc)
         call init_get_isea(isea, jsea)
         ix  = mapsf(isea,1)
         iy  = mapsf(isea,2)
-        if (mapsta(iy,ix) == 1) then
+        if (mapsta(iy,ix) /= 0) then
           sw_wnmean(jsea) = wnmean(jsea)
         endif
       enddo
@@ -902,7 +902,7 @@ subroutine export_fields (gcomp, rc)
         call init_get_isea(isea, jsea)
         ix  = mapsf(isea,1)
         iy  = mapsf(isea,2)
-        if (mapsta(iy,ix) == 1) then
+        if (mapsta(iy,ix) /= 0) then
           sw_taubblx(jsea) = taubbl(jsea,1)
           sw_taubbly(jsea) = taubbl(jsea,2)
         endif
@@ -1243,7 +1243,7 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn, fval)
       call init_get_isea(isea, jsea)
       ix  = mapsf(isea,1)                   ! global ix
       iy  = mapsf(isea,2)                   ! global iy
-      if (mapsta(iy,ix) == 1) then          ! active sea point
+      if (mapsta(iy,ix) /= 0) then          ! active sea point
         sxx1 = 0.0
         syy1 = 0.0
         sxy1 = 0.0
@@ -1271,7 +1271,7 @@ subroutine CalcRadstr2D ( a, sxxn, sxyn, syyn, fval)
         syy1 = syy1 + fte * abyy/cg(nk,isea)
         sxy1 = sxy1 + fte * abxy/cg(nk,isea)
       end if
-      if (mapsta(iy,ix) == 1) then          ! active sea point
+      if (mapsta(iy,ix) /= 0) then          ! active sea point
         sxxn(jsea) = sxx1*dwat*grav
         syyn(jsea) = syy1*dwat*grav
         sxyn(jsea) = sxy1*dwat*grav
@@ -1367,7 +1367,7 @@ subroutine CalcHS (a, hs, fval)
       call init_get_isea(isea, jsea)
       ix  = mapsf(isea,1)                   ! global ix
       iy  = mapsf(isea,2)                   ! global iy
-      if (mapsta(iy,ix) == 1) then          ! active sea point
+      if (mapsta(iy,ix) /= 0) then          ! active sea point
         et = 0.0
         do ik = 1,nk
           factor = dden(ik) / cg(ik,isea)
@@ -1425,7 +1425,7 @@ subroutine CalcBHD (a, bhd, fval)
       call init_get_isea(isea, jsea)
       ix  = mapsf(isea,1)                   ! global ix
       iy  = mapsf(isea,2)                   ! global iy
-      if (mapsta(iy,ix) == 1) then          ! active sea point
+      if (mapsta(iy,ix) /= 0) then          ! active sea point
         ebd = 0.0
         bhd1 = 0.0
         do ik = 1,nk
@@ -1481,7 +1481,7 @@ subroutine CalcStokes(a, us, vs, fval)
       call init_get_isea(isea, jsea)
       ix  = mapsf(isea,1)                   ! global ix
       iy  = mapsf(isea,2)                   ! global iy
-      if (mapsta(iy,ix) == 1) then          ! active sea point
+      if (mapsta(iy,ix) /= 0) then          ! active sea point
         us1 = 0.0
         vs1 = 0.0
         do ik = 1,nk
@@ -1492,8 +1492,8 @@ subroutine CalcStokes(a, us, vs, fval)
             abx = abx + a(ith,ik,jsea)*ecos(ith)
             aby = aby + a(ith,ik,jsea)*esin(ith)
           end do
-          kd = max ( 0.001 , wn(ik,isea) * dw(isea) )
-          if (kd .lt. 6.0) then
+          kd = max ( 1.0e-7 , wn(ik,isea) * dw(isea) ) !0.001
+          if (kd .lt. 18.0) then !6.0
             fkd =  factor / sinh(kd)**2
             ussco = fkd*sig(ik)*wn(ik,isea)*cosh(2.0*kd)
           else
@@ -1545,7 +1545,7 @@ subroutine CalcUVBed(a, ubrx, ubry, fval)
       call init_get_isea(isea, jsea)
       ix  = mapsf(isea,1)                   ! global ix
       iy  = mapsf(isea,2)                   ! global iy
-      if (mapsta(iy,ix) == 1) then          ! active sea point
+      if (mapsta(iy,ix) /= 0) then          ! active sea point
         uba1 = 0.0
         ubd1 = 0.0
         ubr1 = 0.0
@@ -1559,23 +1559,19 @@ subroutine CalcUVBed(a, ubrx, ubry, fval)
             abx = abx + a(ith,ik,jsea)*ecos(ith)
             aby = aby + a(ith,ik,jsea)*esin(ith)
           end do
-          kd = max ( 0.001 , wn(ik,isea) * dw(isea) )
-          if (kd .lt. 6.0) then
+          !kd = max ( 0.001 , wn(ik,isea) * dw(isea) )
+          kd = max ( 1.0e-7, min(18.,wn(ik,isea)) * dw(isea) )
+          !if (kd .lt. 6.0) then
             fkd =  factor / sinh(kd)**2
             ubr1 = ubr1 + ab*sig(ik)**2 * fkd
             uba1 = uba1 + abx*sig(ik)**2 * fkd
             ubd1 = ubd1 + aby*sig(ik)**2 * fkd
-          end if
+          !end if
         end do !ik
-        ubr1 = sqrt(2.0*max(0.0,ubr1))
-        if (ubr1 .ge. 1.0e-7) then
-          ubd1 = atan2(ubd1,uba1)
-        else
-          ubd1 = 0.0
-        end if
-        uba1 = ubr1
-        ubrx(jsea) = uba1*cos(ubd1)
-        ubry(jsea) = uba1*sin(ubd1)
+        ubr1 = sqrt(2.0*max(0.0,ubr1)) !make ubr1 >=0 already
+        ubd1 = atan2(ubd1,uba1)
+        ubrx(jsea) = ubr1*cos(ubd1)
+        ubry(jsea) = ubr1*sin(ubd1)
       else
         ubrx(jsea) = fval
         ubry(jsea) = fval
@@ -1615,7 +1611,7 @@ subroutine CalcTHM (a, thm, fval)
       call init_get_isea(isea, jsea)
       ix  = mapsf(isea,1)                   ! global ix
       iy  = mapsf(isea,2)                   ! global iy
-      if (mapsta(iy,ix) == 1) then          ! active sea point
+      if (mapsta(iy,ix) /= 0) then          ! active sea point
         etx = 0.0
         ety = 0.0
         do ik = 1,nk
@@ -1676,7 +1672,7 @@ subroutine CalcT0M1 (a, t0m1, fval)
       call init_get_isea(isea, jsea)
       ix  = mapsf(isea,1)                   ! global ix
       iy  = mapsf(isea,2)                   ! global iy
-      if (mapsta(iy,ix) == 1) then          ! active sea point
+      if (mapsta(iy,ix) /= 0) then          ! active sea point
         etr = 0.0
         et = 0.0
         do ik = 1,nk