Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:42:03 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
af960ebfb4 Jean*0001 #include "MOM_COMMON_OPTIONS.h"
                0002 
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: MOM_CALC_3D_STRAIN
                0006 
                0007 C     !INTERFACE:
                0008       SUBROUTINE MOM_CALC_3D_STRAIN(
                0009      O        str11, str22, str33, str12, str13, str23,
                0010      I        bi, bj, myThid )
                0011 
                0012 C     !DESCRIPTION:
                0013 C     Calculates the strain tensor of the 3-D flow field
                0014 
                0015 C     !USES:
                0016       IMPLICIT NONE
                0017 
                0018 C     == Global variables ==
                0019 #include "SIZE.h"
                0020 #include "EEPARAMS.h"
                0021 #include "PARAMS.h"
                0022 #include "GRID.h"
                0023 #include "DYNVARS.h"
                0024 
                0025 C     !INPUT PARAMETERS:
                0026 C     bi, bj      :: tile indices
                0027 C     myThid      :: my Thread Id number
                0028       INTEGER bi, bj
                0029       INTEGER myThid
                0030 
                0031 C     !OUTPUT PARAMETERS:
                0032 C     str11       :: strain component Vxx @ grid-cell center
                0033 C     str22       :: strain component Vyy @ grid-cell center
                0034 C     str33       :: strain component Vzz @ grid-cell center
                0035 C     str12       :: strain component Vxy @ grid-cell corner
                0036 C     str13       :: strain component Vxz @ above uVel
                0037 C     str23       :: strain component Vyz @ above vVel
                0038       _RL str11(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0039       _RL str22(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0040       _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0041       _RL str12(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0042       _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0043       _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
                0044 
245cf5c24c Jean*0045 #ifdef ALLOW_SMAG_3D
af960ebfb4 Jean*0046 C     !LOCAL VARIABLES:
                0047 C      i, j, k    :: loop indices
                0048       INTEGER i, j, k
                0049       INTEGER kp1
                0050       _RL maskp1
                0051       LOGICAL freeSlip3d
                0052 CEOP
                0053       freeSlip3d = .NOT.( no_slip_sides .AND. no_slip_bottom )
                0054 
                0055       DO k=1,Nr
                0056        kp1 = MIN(k+1,Nr)
                0057        maskp1 = oneRL
                0058        IF ( k.EQ.Nr ) maskp1 = zeroRL
                0059 
                0060 C-    Fills up array edges:
                0061        i = sNx+OLx
                0062        DO j=1-OLy,sNy+OLy
                0063          str11(i,j,k) = 0. _d 0
                0064        ENDDO
                0065        j = sNy+OLy
                0066        DO i=1-OLx,sNx+OLx
                0067          str22(i,j,k) = 0. _d 0
                0068        ENDDO
                0069        i = 1-OLx
                0070        DO j=1-OLy,sNy+OLy
                0071          str12(i,j,k) = 0. _d 0
                0072          str13(i,j,k) = 0. _d 0
                0073        ENDDO
                0074        j = 1-OLy
                0075        DO i=1-OLx,sNx+OLx
                0076          str12(i,j,k) = 0. _d 0
                0077          str23(i,j,k) = 0. _d 0
                0078        ENDDO
                0079 
                0080 C     str11 = u_x
                0081        DO j=1-OLy,sNy+OLy
                0082         DO i=1-OLx,sNx+OLx-1
245cf5c24c Jean*0083           str11(i,j,k) = recip_dxF(i,j,bi,bj)
                0084      &           *( uVel(i+1, j , k ,bi,bj)-uVel( i , j , k ,bi,bj) )
af960ebfb4 Jean*0085         ENDDO
                0086        ENDDO
                0087 
                0088 C     str22 = v_y
                0089        DO j=1-OLy,sNy+OLy-1
                0090         DO i=1-OLx,sNx+OLx
245cf5c24c Jean*0091           str22(i,j,k) = recip_dyF(i,j,bi,bj)
                0092      &           *( vVel( i ,j+1, k ,bi,bj)-vVel( i , j , k ,bi,bj) )
af960ebfb4 Jean*0093         ENDDO
                0094        ENDDO
                0095 
                0096 C     str33 = w_z
                0097        DO j=1-OLy,sNy+OLy
                0098         DO i=1-OLx,sNx+OLx
245cf5c24c Jean*0099           str33(i,j,k) = recip_drF(k)*rkSign
                0100      &    *( maskp1*wVel( i , j ,kp1,bi,bj)-wVel( i , j , k ,bi,bj) )
af960ebfb4 Jean*0101         ENDDO
                0102        ENDDO
                0103 
                0104 C     str12 = ( u_y + v_x )/2
                0105        DO j=2-OLy,sNy+OLy
                0106         DO i=2-OLx,sNx+OLx
                0107           str12(i,j,k) = halfRL*(
245cf5c24c Jean*0108      &      recip_dyU(i,j,bi,bj)
                0109      &           *( uVel( i , j , k ,bi,bj)-uVel( i ,j-1, k ,bi,bj) )
                0110      &     +recip_dxV(i,j,bi,bj)
                0111      &           *( vVel( i , j , k ,bi,bj)-vVel(i-1, j , k ,bi,bj) )
af960ebfb4 Jean*0112      &                          )
                0113         ENDDO
                0114        ENDDO
                0115 
                0116 C     str13 & str23 special case: k=1
                0117        IF ( k.EQ.1 .AND. freeSlip3d ) THEN
                0118         DO j=1-OLy,sNy+OLy
                0119          DO i=1-OLx,sNx+OLx
                0120           str13(i,j,k) = 0. _d 0
                0121           str23(i,j,k) = 0. _d 0
                0122          ENDDO
                0123         ENDDO
                0124        ELSEIF ( k.EQ.1 ) THEN
                0125 C--    should put surface wind-stress if z-coords; but right in p-coords:
                0126         DO j=1-OLy,sNy+OLy
                0127          DO i=2-OLx,sNx+OLx
                0128           str13(i,j,k) = halfRL*(
245cf5c24c Jean*0129      &      recip_drC(k)*rkSign
                0130      &           *( uVel( i , j , k ,bi,bj)*twoRL )
                0131      &     +recip_dxC(i,j,bi,bj)
                0132      &           *( wVel( i , j , k ,bi,bj)-wVel(i-1, j , k ,bi,bj) )
af960ebfb4 Jean*0133      &                          )
                0134          ENDDO
                0135         ENDDO
                0136         DO j=2-OLy,sNy+OLy
                0137          DO i=1-OLx,sNx+OLx
                0138           str23(i,j,k) = halfRL*(
245cf5c24c Jean*0139      &      recip_drC(k)*rkSign
                0140      &           *( vVel( i , j , k ,bi,bj)*twoRL )
                0141      &     +recip_dyC(i,j,bi,bj)
                0142      &           *( wVel( i , j , k ,bi,bj)-wVel( i ,j-1, k ,bi,bj) )
af960ebfb4 Jean*0143      &                          )
                0144          ENDDO
                0145         ENDDO
                0146        ELSE
                0147 C     str13 = ( u_z + w_x )/2
                0148         DO j=1-OLy,sNy+OLy
                0149          DO i=2-OLx,sNx+OLx
                0150           str13(i,j,k) = halfRL*(
245cf5c24c Jean*0151      &      recip_drC(k)*rkSign
                0152      &           *( uVel( i , j , k ,bi,bj)-uVel( i , j ,k-1 ,bi,bj) )
                0153      &     +recip_dxC(i,j,bi,bj)
                0154      &           *( wVel( i , j , k ,bi,bj)-wVel(i-1, j , k ,bi,bj) )
af960ebfb4 Jean*0155      &                          )
                0156          ENDDO
                0157         ENDDO
                0158 C     str23 = ( v_z + w_y )/2
                0159         DO j=2-OLy,sNy+OLy
                0160          DO i=1-OLx,sNx+OLx
                0161           str23(i,j,k) = halfRL*(
245cf5c24c Jean*0162      &      recip_drC(k)*rkSign
                0163      &           *( vVel( i , j , k ,bi,bj)-vVel( i , j ,k-1,bi,bj) )
                0164      &     +recip_dyC(i,j,bi,bj)
                0165      &           *( wVel( i , j , k ,bi,bj)-wVel( i ,j-1, k ,bi,bj) )
af960ebfb4 Jean*0166      &                          )
                0167          ENDDO
                0168         ENDDO
                0169        ENDIF
                0170 
                0171        IF ( freeSlip3d ) THEN
                0172         DO j=2-OLy,sNy+OLy
                0173          DO i=2-OLx,sNx+OLx
                0174            str12(i,j,k) = str12(i,j,k)
                0175      &                  *maskW(i,j-1,k,bi,bj)*maskW(i,j,k,bi,bj)
                0176          ENDDO
                0177         ENDDO
                0178         IF ( k.GE.2 ) THEN
                0179          DO j=1-OLy,sNy+OLy
                0180           DO i=2-OLx,sNx+OLx
                0181            str13(i,j,k) = str13(i,j,k)
                0182      &                  *maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)
                0183           ENDDO
                0184          ENDDO
                0185          DO j=2-OLy,sNy+OLy
                0186           DO i=1-OLx,sNx+OLx
                0187            str23(i,j,k) = str23(i,j,k)
                0188      &                  *maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)
                0189           ENDDO
                0190          ENDDO
                0191         ENDIF
                0192        ENDIF
                0193 
                0194 C--   end k loop
                0195       ENDDO
                0196 
                0197 C--   fill-up strain tensor component at the very bottom (k=Nr+1)
                0198       k = Nr+1
                0199 
                0200       DO j=1-OLy,sNy+OLy
                0201        DO i=1-OLx,sNx+OLx
                0202          str13(i,j,k) = 0. _d 0
                0203          str23(i,j,k) = 0. _d 0
                0204        ENDDO
                0205       ENDDO
                0206 
                0207       IF ( .NOT.freeSlip3d ) THEN
                0208 
                0209 C     str13 = ( u_z + w_x )/2
                0210        DO j=1-OLy,sNy+OLy
                0211         DO i=2-OLx,sNx+OLx
                0212           str13(i,j,k) =
245cf5c24c Jean*0213      &         recip_drF(Nr)*rkSign
                0214 c    &      recip_drC(k)*rkSign
                0215      &            *( 0. _d 0 - uVel( i , j ,k-1 ,bi,bj) )
af960ebfb4 Jean*0216         ENDDO
                0217        ENDDO
                0218 
                0219 C     str23 = ( v_z + w_y )/2
                0220        DO j=2-OLy,sNy+OLy
                0221         DO i=1-OLx,sNx+OLx
                0222         str23(i,j,k) =
245cf5c24c Jean*0223      &         recip_drF(Nr)*rkSign
                0224 c    &      recip_drC(k)*rkSign
                0225      &            *( 0. _d 0 - vVel( i , j ,k-1,bi,bj) )
af960ebfb4 Jean*0226         ENDDO
                0227        ENDDO
                0228 
                0229       ENDIF
                0230 
                0231 C     Special stuff for Cubed Sphere
                0232 c     IF (useCubedSphereExchange) THEN
                0233 c      STOP 'S/R MOM_CALC_3D_STRAIN: should not be used on the cube!'
                0234 c     ENDIF
                0235 
245cf5c24c Jean*0236 #endif /* ALLOW_SMAG_3D */
af960ebfb4 Jean*0237       RETURN
                0238       END