--- MITgcm/pkg/seaice/seaice_calc_strainrates.F	2009/03/18 10:26:49	1.10
+++ MITgcm/pkg/seaice/seaice_calc_strainrates.F	2009/03/18 12:58:17	1.11
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_calc_strainrates.F,v 1.10 2009/03/18 10:26:49 mlosch Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_calc_strainrates.F,v 1.11 2009/03/18 12:58:17 mlosch Exp $
 C $Name:  $
 
 #include "SEAICE_OPTIONS.h"
@@ -22,6 +22,7 @@
 #include "PARAMS.h"
 #include "GRID.h"
 #include "SEAICE_PARAMS.h"
+#include "SEAICE.h"
 
 #ifdef ALLOW_AUTODIFF_TAMC
 # include "tamc.h"
@@ -54,10 +55,69 @@
       INTEGER i, j, bi, bj
 C  hFacU, hFacV - determine the no-slip boundary condition
       INTEGER k
-      _RS hFacU, hFacV
+      _RS hFacU, hFacV, noSlipFac
 
       k = 1
+      noSlipFac = 0. _d 0
+      IF ( SEAICE_no_slip ) noSlipFac = 1. _d 0
 C
+#ifndef SEAICE_OLD_AND_BAD_DISCRETIZATION
+      DO bj=myByLo(myThid),myByHi(myThid)
+       DO bi=myBxLo(myThid),myBxHi(myThid)
+        DO j=1-Oly,sNy+Oly-1
+         DO i=1-Olx,sNx+Olx-1
+C     evaluate strain rates
+          e11(I,J,bi,bj) = _recip_dxF(I,J,bi,bj) *
+     &         (uFld(I+1,J,1,bi,bj)-uFld(I,J,1,bi,bj))
+     &         +HALF*
+     &         (vFld(I,J,1,bi,bj)+vFld(I,J+1,1,bi,bj))
+     &         * k2AtC(I,J,bi,bj)
+          e22(I,J,bi,bj) = _recip_dyF(I,J,bi,bj) *
+     &         (vFld(I,J+1,1,bi,bj)-vFld(I,J,1,bi,bj))
+     &         +HALF*
+     &         (uFld(I,J,1,bi,bj)+uFld(I+1,J,1,bi,bj))
+     &         * k1AtC(I,J,bi,bj)
+C     one metric term is missing
+         ENDDO
+        ENDDO
+        DO j=1-Oly+1,sNy+Oly
+         DO i=1-Olx+1,sNx+Olx
+          hFacU = _maskW(i,j,k,bi,bj) - _maskW(i,j-1,k,bi,bj)
+          hFacV = _maskS(i,j,k,bi,bj) - _maskS(i-1,j,k,bi,bj)
+          e12(I,J,bi,bj) = HALF*(
+     &           ( uFld(I,J,1,bi,bj) - uFld(I  ,J-1,1,bi,bj) )
+     &         * _recip_dyU(I,J,bi,bj)
+     &         + ( vFld(I,J,1,bi,bj) - vFld(I-1,J  ,1,bi,bj) )
+     &         * _recip_dxV(I,J,bi,bj) 
+     &         - k1AtZ(I,J,bi,bj)
+     &         * 0.5 _d 0 * (vFld(I,J,1,bi,bj)+vFld(I-1,J  ,1,bi,bj))
+     &         - k2AtZ(I,J,bi,bj)
+     &         * 0.5 _d 0 * (uFld(I,J,1,bi,bj)+uFld(I  ,J-1,1,bi,bj))
+     &         )
+     &         *maskC(I  ,J  ,k,bi,bj)*maskC(I-1,J  ,k,bi,bj)
+     &         *maskC(I  ,J-1,k,bi,bj)*maskC(I-1,J-1,k,bi,bj)
+     &         + 2.0 _d 0 * noSlipFac * (
+     &           ( uFld(I,J,1,bi,bj) + uFld(I  ,J-1,1,bi,bj) )
+     &         * _recip_dyU(I,J,bi,bj) * hFacU
+     &         + ( vFld(I,J,1,bi,bj) + vFld(I-1,J  ,1,bi,bj) )
+     &         * _recip_dxV(I,J,bi,bj) * hFacV
+     &         )
+C     no slip at the boundary implies u(j)+u(j-1)=0 and v(i)+v(i-1)=0
+C     accross the boundary; this is already accomplished by masking so
+C     that the following lines are not necessary
+c$$$     &         - hFacV * k1AtZ(I,J,bi,bj)
+c$$$     &         * 0.5 _d 0 * (vFld(I,J,1,bi,bj)+vFld(I-1,J  ,1,bi,bj))
+c$$$     &         - hFacU * k2AtZ(I,J,bi,bj)
+c$$$     &         * 0.5 _d 0 * (uFld(I,J,1,bi,bj)+uFld(I  ,J-1,1,bi,bj))
+         ENDDO
+        ENDDO
+
+c$$$        ENDIF
+       ENDDO
+      ENDDO
+#else 
+C     this the old and incomplete discretization, here I also erroneously
+C     used finite-volumes to discretize the strain rates
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-Oly,sNy+Oly-1
@@ -115,6 +175,7 @@
         ENDIF
        ENDDO
       ENDDO
+#endif /* SEAICE_OLD_AND_BAD_DISCRETIZATION */
 #endif /* SEAICE_ALLOW_DYNAMICS */
 #endif /* SEAICE_CGRID */
       RETURN