pkg/mom_vecinv/mom_vi_u_coriolis.F

cec2469d72 Alis*0001 #include "MOM_VECINV_OPTIONS.h"
aea29c8517 Alis*0002
355718db32 Jean*0003 CBOP
                0004 C     !ROUTINE: MOM_VI_U_CORIOLIS
                0005 C     !INTERFACE:
5d7e0a8948 Jean*0006       SUBROUTINE MOM_VI_U_CORIOLIS(
31fb0e0e6d Jean*0007      I                     bi, bj, k,
3370e71df9 Mart*0008      I                     selectVortScheme, useJamartMomAdv,
355718db32 Jean*0009      I                     vFld, omega3, hFacZ, r_hFacZ,
                0010      O                     uCoriolisTerm,
                0011      I                     myThid )
                0012 C     !DESCRIPTION: \bv
cab1a69b8a Jean*0013 C     *==========================================================*
                0014 C     | S/R MOM_VI_U_CORIOLIS
355718db32 Jean*0015 C     |==========================================================*
                0016 C     | o Calculate flux (in Y-dir.) of vorticity at U point
                0017 C     |   using 2nd order interpolation
cab1a69b8a Jean*0018 C     *==========================================================*
355718db32 Jean*0019 C     \ev
                0020
                0021 C     !USES:
                0022       IMPLICIT NONE
aea29c8517 Alis*0023
                0024 C     == Global variables ==
                0025 #include "SIZE.h"
                0026 #include "EEPARAMS.h"
                0027 #include "GRID.h"
                0028
355718db32 Jean*0029 C     !INPUT/OUTPUT PARAMETERS:
aea29c8517 Alis*0030 C     == Routine arguments ==
355718db32 Jean*0031       INTEGER bi, bj, k
                0032       _RL     vFld     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0033       _RL     omega3   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0034       _RS     hFacZ    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0035       _RS     r_hFacZ  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
aea29c8517 Alis*0036       _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
3370e71df9 Mart*0037       INTEGER selectVortScheme
                0038       LOGICAL useJamartMomAdv
aea29c8517 Alis*0039       INTEGER myThid
355718db32 Jean*0040 CEOP
aea29c8517 Alis*0041
                0042 C     == Local variables ==
31fb0e0e6d Jean*0043 C     msgBuf :: Informational/error message buffer
355718db32 Jean*0044       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0045       INTEGER i, j
                0046       _RL     vBarXY, vBarXm, vBarXp
                0047       _RL     vort3u
4e5f84a272 Jean*0048       _RL     vort3mj, vort3ij, vort3mp, vort3ip
                0049       _RL     oneThird, tmpFac
355718db32 Jean*0050       _RS     epsil
aea29c8517 Alis*0051
5d7e0a8948 Jean*0052       epsil = 1. _d -9
4e5f84a272 Jean*0053       tmpFac = 1. _d 0
                0054 c     oneThird = 1. _d 0 / ( 1. _d 0 + 2.*tmpFac )
                0055       oneThird = 1. _d 0 / 3. _d 0
aea29c8517 Alis*0056
355718db32 Jean*0057       IF ( selectVortScheme.EQ.0 ) THEN
                0058 C--   using enstrophy conserving scheme (Shallow-Water Eq.) by Sadourny, JAS 75
                0059
31fb0e0e6d Jean*0060        DO j=1-OLy,sNy+OLy-1
                0061         DO i=2-OLx,sNx+OLx
                0062          vBarXY = 0.25 _d 0*(
                0063      &      ( vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
                0064      &      + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
                0065      &     +( vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
                0066      &      + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
                0067      &                      )
                0068          vort3u = halfRL*( omega3(i,j)*r_hFacZ(i,j)
                0069      &                   + omega3(i,j+1)*r_hFacZ(i,j+1) )
                0070          uCoriolisTerm(i,j) = +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
                0071      &                               * _maskW(i,j,k,bi,bj)
355718db32 Jean*0072         ENDDO
                0073        ENDDO
                0074
                0075       ELSEIF ( selectVortScheme.EQ.1 ) THEN
                0076 C--   same as above, with different formulation (relatively to hFac)
                0077
31fb0e0e6d Jean*0078        DO j=1-OLy,sNy+OLy-1
                0079         DO i=2-OLx,sNx+OLx
                0080          vBarXY = halfRL*(
                0081      &      ( vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i, j )
                0082      &      + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacZ(i, j ) )
                0083      &     +( vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacZ(i,j+1)
                0084      &      + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacZ(i,j+1) )
                0085      &                     )/MAX( epsil, hFacZ(i,j)+hFacZ(i,j+1) )
                0086          vort3u = halfRL*( omega3(i,j) + omega3(i,j+1) )
                0087          uCoriolisTerm(i,j) = +vort3u*vBarXY*recip_dxC(i,j,bi,bj)
                0088      &                               * _maskW(i,j,k,bi,bj)
355718db32 Jean*0089         ENDDO
                0090        ENDDO
                0091
                0092       ELSEIF ( selectVortScheme.EQ.2 ) THEN
                0093 C--   using energy conserving scheme (used by Sadourny in JAS 75 paper)
                0094
31fb0e0e6d Jean*0095        DO j=1-OLy,sNy+OLy-1
                0096         DO i=2-OLx,sNx+OLx
                0097          vBarXm = halfRL*(
355718db32 Jean*0098      &       vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
31fb0e0e6d Jean*0099      &     + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
                0100          vBarXp = halfRL*(
355718db32 Jean*0101      &       vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
31fb0e0e6d Jean*0102      &     + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
                0103          vort3u = ( vBarXm*r_hFacZ(i, j )*omega3(i, j )
                0104      &            + vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1)
                0105      &            )*halfRL
                0106          uCoriolisTerm(i,j) = +vort3u*recip_dxC(i,j,bi,bj)
                0107      &                               * _maskW(i,j,k,bi,bj)
355718db32 Jean*0108         ENDDO
                0109        ENDDO
                0110
4e5f84a272 Jean*0111       ELSEIF ( selectVortScheme.EQ.3 ) THEN
                0112 C--   using energy & enstrophy conserving scheme
                0113 C     (from Sadourny, described by Burridge & Haseler, ECMWF Rep.4, 1977)
                0114
                0115 C     domain where uCoriolisTerm is valid :
                0116 C     [ 3-Olx : sNx+Olx-1 ] x [ 2-Oly : sNy+Oly-1 ]
                0117 C     (=> might need overlap of 3 if using CD-scheme)
31fb0e0e6d Jean*0118        DO j=1-OLy,sNy+OLy-1
                0119         DO i=2-OLx,sNx+OLx-1
                0120          vort3mj = ( r_hFacZ(i, j )*omega3(i, j )
                0121      &             +(r_hFacZ(i,j+1)*omega3(i,j+1)
                0122      &              +r_hFacZ(i-1,j)*omega3(i-1,j)
                0123      &             ))*oneThird
                0124 c    &             )*tmpFac)*oneThird
4e5f84a272 Jean*0125      &      *vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj)
31fb0e0e6d Jean*0126          vort3ij = ( r_hFacZ(i, j )*omega3(i, j )
                0127      &             +(r_hFacZ(i,j+1)*omega3(i,j+1)
                0128      &              +r_hFacZ(i+1,j)*omega3(i+1,j)
                0129      &             ))*oneThird
                0130 c    &             )*tmpFac)*oneThird
4e5f84a272 Jean*0131      &      *vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
31fb0e0e6d Jean*0132          vort3mp = ( r_hFacZ(i,j+1)*omega3(i,j+1)
                0133      &             +(r_hFacZ(i, j )*omega3(i, j )
                0134      &              +r_hFacZ(i-1,j+1)*omega3(i-1,j+1)
                0135      &             ))*oneThird
                0136 c    &             )*tmpFac)*oneThird
4e5f84a272 Jean*0137      &      *vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj)
31fb0e0e6d Jean*0138          vort3ip = ( r_hFacZ(i,j+1)*omega3(i,j+1)
                0139      &             +(r_hFacZ(i, j )*omega3(i, j )
                0140      &              +r_hFacZ(i+1,j+1)*omega3(i+1,j+1)
                0141      &             ))*oneThird
                0142 c    &             )*tmpFac)*oneThird
4e5f84a272 Jean*0143      &      *vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
                0144 C---
31fb0e0e6d Jean*0145          uCoriolisTerm(i,j) = +( (vort3mj+vort3ij)+(vort3mp+vort3ip) )
                0146      &                      *0.25 _d 0 *recip_dxC(i,j,bi,bj)
                0147      &                                 * _maskW(i,j,k,bi,bj)
                0148         ENDDO
                0149        ENDDO
                0150
                0151       ELSEIF ( selectVortScheme.EQ.4 ) THEN
                0152 C--   using energy conserving scheme without r_hFacZ factor
                0153
                0154        DO j=1-OLy,sNy+OLy-1
                0155         DO i=2-OLx,sNx+OLx
                0156          vBarXm = halfRL*(
                0157      &       vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj)
                0158      &     + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) )
                0159          vBarXp = halfRL*(
                0160      &       vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj)
                0161      &     + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) )
                0162          vort3u = ( vBarXm*omega3(i, j )
                0163      &            + vBarXp*omega3(i,j+1)
                0164      &            )*halfRL
                0165          uCoriolisTerm(i,j) = +vort3u*recip_dxC(i,j,bi,bj)
                0166      &                               *recip_hFacW(i,j,k,bi,bj)
4e5f84a272 Jean*0167         ENDDO
                0168        ENDDO
                0169
355718db32 Jean*0170       ELSE
                0171         WRITE(msgBuf,'(A,I5,A)')
                0172      &   'MOM_VI_U_CORIOLIS: selectVortScheme=', selectVortScheme,
                0173      &   ' not implemented'
                0174         CALL PRINT_ERROR( msgBuf, myThid )
                0175         STOP 'ABNORMAL END: S/R MOM_VI_U_CORIOLIS'
                0176
                0177       ENDIF
                0178
                0179       IF ( useJamartMomAdv ) THEN
31fb0e0e6d Jean*0180        DO j=1-OLy,sNy+OLy-1
                0181         DO i=2-OLx,sNx+OLx-1
355718db32 Jean*0182          uCoriolisTerm(i,j) = uCoriolisTerm(i,j)
                0183      &           * 4. _d 0 * _hFacW(i,j,k,bi,bj)
                0184      &           / MAX( epsil,
                0185      &                  (_hFacS(i, j ,k,bi,bj)+_hFacS(i-1, j ,k,bi,bj))
                0186      &                 +(_hFacS(i,j+1,k,bi,bj)+_hFacS(i-1,j+1,k,bi,bj))
                0187      &                )
                0188         ENDDO
aea29c8517 Alis*0189        ENDDO
355718db32 Jean*0190       ENDIF
aea29c8517 Alis*0191
                0192       RETURN
                0193       END