C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_calc_3d_strain.F,v 1.2 2013/11/06 00:37:11 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)
#ifdef ALLOW_SMAG_3D
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
#endif /* ALLOW_SMAG_3D */
RETURN
END