--- MITgcm/pkg/seaice/seaice_calc_strainrates.F	2009/06/24 08:23:38	1.14
+++ MITgcm/pkg/seaice/seaice_calc_strainrates.F	2009/10/23 08:10:16	1.15
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_calc_strainrates.F,v 1.14 2009/06/24 08:23:38 mlosch Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_calc_strainrates.F,v 1.15 2009/10/23 08:10:16 mlosch Exp $
 C $Name:  $
 
 #include "SEAICE_OPTIONS.h"
@@ -51,9 +51,17 @@
 C     === Local variables ===
 C     i,j,bi,bj - Loop counters
       INTEGER i, j, bi, bj
-C  hFacU, hFacV - determine the no-slip boundary condition
+C     hFacU, hFacV - determine the no-slip boundary condition
       INTEGER k
       _RS hFacU, hFacV, noSlipFac
+C     auxillary variables that help writing code that
+C     vectorizes even after TAFization
+      _RL dudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL dvdy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL dudy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL dvdx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL uave (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL vave (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 
       k = 1
       noSlipFac = 0. _d 0
@@ -62,54 +70,69 @@
 #ifndef SEAICE_OLD_AND_BAD_DISCRETIZATION
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
+C     abbreviations on C-points, need to do them in separate loops
+C     for vectorization
         DO j=1-Oly,sNy+Oly-1
          DO i=1-Olx,sNx+Olx-1
-C     evaluate strain rates
-          e11Loc(I,J,bi,bj) = _recip_dxF(I,J,bi,bj) *
+          dudx(I,J) = _recip_dxF(I,J,bi,bj) *
      &         (uFld(I+1,J,bi,bj)-uFld(I,J,bi,bj))
-     &         +HALF*
-     &         (vFld(I,J,bi,bj)+vFld(I,J+1,bi,bj))
-     &         * k2AtC(I,J,bi,bj)
-          e22Loc(I,J,bi,bj) = _recip_dyF(I,J,bi,bj) *
+          uave(I,J) = 0.5 _d 0 * (uFld(I,J,bi,bj)+uFld(I+1,J,bi,bj))
+         ENDDO
+        ENDDO
+        DO j=1-Oly,sNy+Oly-1
+         DO i=1-Olx,sNx+Olx-1
+          dvdy(I,J) = _recip_dyF(I,J,bi,bj) *
      &         (vFld(I,J+1,bi,bj)-vFld(I,J,bi,bj))
-     &         +HALF*
-     &         (uFld(I,J,bi,bj)+uFld(I+1,J,bi,bj))
-     &         * k1AtC(I,J,bi,bj)
+          vave(I,J) = 0.5 _d 0 * (vFld(I,J,bi,bj)+vFld(I,J+1,bi,bj))
+         ENDDO
+        ENDDO
+C     evaluate strain rates at C-points
+        DO j=1-Oly,sNy+Oly-1
+         DO i=1-Olx,sNx+Olx-1
+          e11Loc(I,J,bi,bj) = dudx(I,J) + vave(I,J) * k2AtC(I,J,bi,bj)
+          e22Loc(I,J,bi,bj) = dvdy(I,J) + uave(I,J) * k1AtC(I,J,bi,bj)
          ENDDO
         ENDDO
+C     abbreviations at Z-points, need to do them in separate loops
+C     for vectorization
         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)
-          e12Loc(I,J,bi,bj) = HALF*(
-     &           ( uFld(I,J,bi,bj) - uFld(I  ,J-1,bi,bj) )
+          dudy(I,J) = ( uFld(I,J,bi,bj) - uFld(I  ,J-1,bi,bj) )
      &         * _recip_dyU(I,J,bi,bj)
-     &         + ( vFld(I,J,bi,bj) - vFld(I-1,J  ,bi,bj) )
+          uave(I,J) = 0.5 _d 0 * (uFld(I,J,bi,bj)+uFld(I  ,J-1,bi,bj))
+         ENDDO
+        ENDDO
+        DO j=1-Oly+1,sNy+Oly
+         DO i=1-Olx+1,sNx+Olx
+          dvdx(I,J) = ( vFld(I,J,bi,bj) - vFld(I-1,J  ,bi,bj) )
      &         * _recip_dxV(I,J,bi,bj) 
-     &         - k1AtZ(I,J,bi,bj)
-     &         * 0.5 _d 0 * (vFld(I,J,bi,bj)+vFld(I-1,J  ,bi,bj))
-     &         - k2AtZ(I,J,bi,bj)
-     &         * 0.5 _d 0 * (uFld(I,J,bi,bj)+uFld(I  ,J-1,bi,bj))
+          vave(I,J) = 0.5 _d 0 * (vFld(I,J,bi,bj)+vFld(I-1,J  ,bi,bj))
+         ENDDO
+        ENDDO
+C     evaluate strain rates at Z-points
+        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)
+          e12Loc(I,J,bi,bj) = 0.5 _d 0 * (
+     &         dudy(I,J) + dvdx(I,J)
+     &         - k1AtZ(I,J,bi,bj) * vave(I,J)
+     &         - k2AtZ(I,J,bi,bj) * uave(I,J)
      &         )
      &         *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,bi,bj) + uFld(I  ,J-1,bi,bj) )
-     &         * _recip_dyU(I,J,bi,bj) * hFacU
-     &         + ( vFld(I,J,bi,bj) + vFld(I-1,J  ,bi,bj) )
-     &         * _recip_dxV(I,J,bi,bj) * hFacV
+     &           2.0 _d 0 * uave(I,J) * _recip_dyU(I,J,bi,bj) * hFacU
+     &         + 2.0 _d 0 * vave(I,J) * _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,bi,bj)+vFld(I-1,J  ,bi,bj))
-c$$$     &         - hFacU * k2AtZ(I,J,bi,bj)
-c$$$     &         * 0.5 _d 0 * (uFld(I,J,bi,bj)+uFld(I  ,J-1,bi,bj))
+c$$$     &         - hFacV * k1AtZ(I,J,bi,bj) * vave(I,J)
+c$$$     &         - hFacU * k2AtZ(I,J,bi,bj) * uave(I,J)
          ENDDO
         ENDDO
 
-c$$$        ENDIF
        ENDDO
       ENDDO
 #else