C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_common/mom_calc_3d_strain.F,v 1.1 2013/11/05 13:31:50 jmc Exp $
C $Name:  $

#include "MOM_COMMON_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: MOM_CALC_3D_STRAIN

C     !INTERFACE:
      SUBROUTINE MOM_CALC_3D_STRAIN(
     O        str11, str22, str33, str12, str13, str23,
     I        bi, bj, myThid )

C     !DESCRIPTION:
C     Calculates the strain tensor of the 3-D flow field

C     !USES:
      IMPLICIT NONE

C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"

C     !INPUT PARAMETERS:
C     bi, bj      :: tile indices
C     myThid      :: my Thread Id number
      INTEGER bi, bj
      INTEGER myThid

C     !OUTPUT PARAMETERS:
C     str11       :: strain component Vxx @ grid-cell center
C     str22       :: strain component Vyy @ grid-cell center
C     str33       :: strain component Vzz @ grid-cell center
C     str12       :: strain component Vxy @ grid-cell corner
C     str13       :: strain component Vxz @ above uVel
C     str23       :: strain component Vyz @ above vVel
      _RL str11(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL str22(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL str12(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
      _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)

C     !LOCAL VARIABLES:
C      i, j, k    :: loop indices
      INTEGER i, j, k
      INTEGER kp1
      _RL maskp1
      LOGICAL freeSlip3d
CEOP
      freeSlip3d = .NOT.( no_slip_sides .AND. no_slip_bottom )

      DO k=1,Nr
       kp1 = MIN(k+1,Nr)
       maskp1 = oneRL
       IF ( k.EQ.Nr ) maskp1 = zeroRL

C-    Fills up array edges:
       i = sNx+OLx
       DO j=1-OLy,sNy+OLy
         str11(i,j,k) = 0. _d 0
       ENDDO
       j = sNy+OLy
       DO i=1-OLx,sNx+OLx
         str22(i,j,k) = 0. _d 0
       ENDDO
       i = 1-OLx
       DO j=1-OLy,sNy+OLy
         str12(i,j,k) = 0. _d 0
         str13(i,j,k) = 0. _d 0
       ENDDO
       j = 1-OLy
       DO i=1-OLx,sNx+OLx
         str12(i,j,k) = 0. _d 0
         str23(i,j,k) = 0. _d 0
       ENDDO

C     str11 = u_x
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx-1
          str11(i,j,k) =
     &        recip_dxF(i,j,bi,bj)*(
     &              uVel(i+1, j , k ,bi,bj)-uVel( i , j , k ,bi,bj) )
        ENDDO
       ENDDO

C     str22 = v_y
       DO j=1-OLy,sNy+OLy-1
        DO i=1-OLx,sNx+OLx
          str22(i,j,k) =
     &        recip_dyF(i,j,bi,bj)*(
     &              vVel( i ,j+1, k ,bi,bj)-vVel( i , j , k ,bi,bj) )
        ENDDO
       ENDDO

C     str33 = w_z
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
          str33(i,j,k) =
     &        recip_drF(k)*rkSign*(
     &       maskp1*wVel( i , j ,kp1,bi,bj)-wVel( i , j , k ,bi,bj) )
        ENDDO
       ENDDO

C     str12 = ( u_y + v_x )/2
       DO j=2-OLy,sNy+OLy
        DO i=2-OLx,sNx+OLx
          str12(i,j,k) = halfRL*(
     &      recip_dyU(i,j,bi,bj)*(
     &              uVel( i , j , k ,bi,bj)-uVel( i ,j-1, k ,bi,bj) )
     &     +recip_dxV(i,j,bi,bj)*(
     &              vVel( i , j , k ,bi,bj)-vVel(i-1, j , k ,bi,bj) )
     &                          )
        ENDDO
       ENDDO

C     str13 & str23 special case: k=1
       IF ( k.EQ.1 .AND. freeSlip3d ) THEN
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          str13(i,j,k) = 0. _d 0
          str23(i,j,k) = 0. _d 0
         ENDDO
        ENDDO
       ELSEIF ( k.EQ.1 ) THEN
C--    should put surface wind-stress if z-coords; but right in p-coords:
        DO j=1-OLy,sNy+OLy
         DO i=2-OLx,sNx+OLx
          str13(i,j,k) = halfRL*(
     &      recip_drC(k)*rkSign*(
     &              uVel( i , j , k ,bi,bj)*twoRL )
     &     +recip_dxC(i,j,bi,bj)*(
     &              wVel( i , j , k ,bi,bj)-wVel(i-1, j , k ,bi,bj) )
     &                          )
         ENDDO
        ENDDO
        DO j=2-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          str23(i,j,k) = halfRL*(
     &      recip_drC(k)*rkSign*(
     &              vVel( i , j , k ,bi,bj)*twoRL )
     &     +recip_dyC(i,j,bi,bj)*(
     &              wVel( i , j , k ,bi,bj)-wVel( i ,j-1, k ,bi,bj) )
     &                          )
         ENDDO
        ENDDO
       ELSE
C     str13 = ( u_z + w_x )/2
        DO j=1-OLy,sNy+OLy
         DO i=2-OLx,sNx+OLx
          str13(i,j,k) = halfRL*(
     &      recip_drC(k)*rkSign*(
     &              uVel( i , j , k ,bi,bj)-uVel( i , j ,k-1 ,bi,bj) )
     &     +recip_dxC(i,j,bi,bj)*(
     &              wVel( i , j , k ,bi,bj)-wVel(i-1, j , k ,bi,bj) )
     &                          )
         ENDDO
        ENDDO
C     str23 = ( v_z + w_y )/2
        DO j=2-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          str23(i,j,k) = halfRL*(
     &      recip_drC(k)*rkSign*(
     &              vVel( i , j , k ,bi,bj)-vVel( i , j ,k-1,bi,bj) )
     &     +recip_dyC(i,j,bi,bj)*(
     &              wVel( i , j , k ,bi,bj)-wVel( i ,j-1, k ,bi,bj) )
     &                          )
         ENDDO
        ENDDO
       ENDIF

       IF ( freeSlip3d ) THEN
        DO j=2-OLy,sNy+OLy
         DO i=2-OLx,sNx+OLx
           str12(i,j,k) = str12(i,j,k)
     &                  *maskW(i,j-1,k,bi,bj)*maskW(i,j,k,bi,bj)
         ENDDO
        ENDDO
        IF ( k.GE.2 ) THEN
         DO j=1-OLy,sNy+OLy
          DO i=2-OLx,sNx+OLx
           str13(i,j,k) = str13(i,j,k)
     &                  *maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)
          ENDDO
         ENDDO
         DO j=2-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           str23(i,j,k) = str23(i,j,k)
     &                  *maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDIF
       ENDIF

C--   end k loop
      ENDDO

C--   fill-up strain tensor component at the very bottom (k=Nr+1)
      k = Nr+1

      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
         str13(i,j,k) = 0. _d 0
         str23(i,j,k) = 0. _d 0
       ENDDO
      ENDDO

      IF ( .NOT.freeSlip3d ) THEN

C     str13 = ( u_z + w_x )/2
       DO j=1-OLy,sNy+OLy
        DO i=2-OLx,sNx+OLx
          str13(i,j,k) =
     &         recip_drF(Nr)*rkSign*(
c    &      recip_drC(k)*rkSign*(
     &               0. _d 0 - uVel( i , j ,k-1 ,bi,bj) )
        ENDDO
       ENDDO

C     str23 = ( v_z + w_y )/2
       DO j=2-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
        str23(i,j,k) =
     &         recip_drF(Nr)*rkSign*(
c    &      recip_drC(k)*rkSign*(
     &               0. _d 0 - vVel( i , j ,k-1,bi,bj) )
        ENDDO
       ENDDO

      ENDIF

C     Special stuff for Cubed Sphere
c     IF (useCubedSphereExchange) THEN
c      STOP 'S/R MOM_CALC_3D_STRAIN: should not be used on the cube!'
c     ENDIF

      RETURN
      END