Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:43:52 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
6d54cf9ca1 Ed H*0001 #include "SHAP_FILT_OPTIONS.h"
bd0c662004 Jean*0002 #undef CALC_CS_CORNER_EXTENDED
d88f4c0b16 Jean*0003 
d4399095cb Jean*0004       SUBROUTINE SHAP_FILT_RELVORT3(
d88f4c0b16 Jean*0005      I        bi,bj,k,
                0006      I        uFld, vFld, hFacZ,
                0007      O        vort3,
                0008      I        myThid)
                0009       IMPLICIT NONE
                0010 C     *==========================================================*
                0011 C     | S/R SHAP_FILT_RELVORT3
                0012 C     *==========================================================*
d4399095cb Jean*0013 C     | copied from S/R MOM_CALC_RELVORT3, except:
                0014 C     |  a) do not apply any masking
                0015 C     |  b) use more uniform rAz at CS-grid corners
d88f4c0b16 Jean*0016 C     *==========================================================*
                0017 
                0018 C     == Global variables ==
                0019 #include "SIZE.h"
                0020 #include "EEPARAMS.h"
                0021 #include "PARAMS.h"
                0022 #include "GRID.h"
d4399095cb Jean*0023 #ifdef ALLOW_EXCH2
f9f661930b Jean*0024 #include "W2_EXCH2_SIZE.h"
d4399095cb Jean*0025 #include "W2_EXCH2_TOPOLOGY.h"
                0026 #endif /* ALLOW_EXCH2 */
d88f4c0b16 Jean*0027 C     == Routine arguments ==
                0028 C     myThid - Instance number for this innvocation of CALC_MOM_RHS
                0029       INTEGER bi,bj,k
                0030       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0031       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0032       _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0033       _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0034       INTEGER myThid
                0035 
                0036 C     == Local variables ==
                0037       INTEGER i,j
d4399095cb Jean*0038       LOGICAL northWestCorner, northEastCorner,
                0039      &        southWestCorner, southEastCorner
bd0c662004 Jean*0040       INTEGER myFace
d4399095cb Jean*0041 #ifdef ALLOW_EXCH2
                0042       INTEGER myTile
                0043 #endif /* ALLOW_EXCH2 */
                0044       _RL AZcorner
                0045 
41b8eb6e45 Jean*0046 #ifdef ALLOW_AUTODIFF
                0047       DO J=1-OLy,sNy+OLy
                0048        DO I=1-OLx,sNx+OLx
d4399095cb Jean*0049         vort3(I,J) = 0. _d 0
                0050        ENDDO
                0051       ENDDO
                0052 #endif
d88f4c0b16 Jean*0053 
41b8eb6e45 Jean*0054       DO J=2-OLy,sNy+OLy
                0055        DO I=2-OLx,sNx+OLx
d88f4c0b16 Jean*0056 
                0057 C       Horizontal curl of flow field - ignoring lopping factors
                0058         vort3(I,J)=
                0059      &      recip_rAz(I,J,bi,bj)*(
bd0c662004 Jean*0060      &      ( vFld(I,J)*dyC(I,J,bi,bj)
                0061      &       -vFld(I-1,J)*dyC(I-1,J,bi,bj) )
                0062      &     -( uFld(I,J)*dxC(I,J,bi,bj)
                0063      &       -uFld(I,J-1)*dxC(I,J-1,bi,bj) )
d88f4c0b16 Jean*0064      &                           )
                0065 
                0066 C       Horizontal curl of flow field - including lopping factors
d4399095cb Jean*0067 c       IF (hFacZ(i,j).NE.0.) THEN
d88f4c0b16 Jean*0068 c        vort3(I,J)=
                0069 c    &      recip_rAz(I,J,bi,bj)*(
                0070 c    &      vFld(I,J)*dyc(I,J,bi,bj)*_hFacW(i,j,k,bi,bj)
                0071 c    &     -vFld(I-1,J)*dyc(I-1,J,bi,bj)*_hFacW(i-1,j,k,bi,bj)
                0072 c    &     -uFld(I,J)*dxc(I,J,bi,bj)*_hFacS(i,j,k,bi,bj)
                0073 c    &     +uFld(I,J-1)*dxc(I,J-1,bi,bj)*_hFacS(i,j-1,k,bi,bj)
                0074 c    &                           )
                0075 c    &                            /hFacZ(i,j)
d4399095cb Jean*0076 c       ELSE
                0077 c        vort3(I,J)=0.
                0078 c       ENDIF
d88f4c0b16 Jean*0079 
                0080        ENDDO
                0081       ENDDO
d4399095cb Jean*0082 
                0083 C     RETURN
                0084 
d88f4c0b16 Jean*0085 C     Special stuff for Cubed Sphere
                0086       IF (useCubedSphereExchange) THEN
d4399095cb Jean*0087        AZcorner = 0.75 _d 0
                0088 #ifdef ALLOW_EXCH2
c424ee7cc7 Jean*0089        myTile = W2_myTileList(bi,bj)
bd0c662004 Jean*0090        myFace = exch2_myFace(myTile)
                0091          southWestCorner = exch2_isWedge(myTile).EQ.1
                0092      &               .AND. exch2_isSedge(myTile).EQ.1
                0093          southEastCorner = exch2_isEedge(myTile).EQ.1
                0094      &               .AND. exch2_isSedge(myTile).EQ.1
                0095          northEastCorner = exch2_isEedge(myTile).EQ.1
                0096      &               .AND. exch2_isNedge(myTile).EQ.1
                0097          northWestCorner = exch2_isWedge(myTile).EQ.1
                0098      &               .AND. exch2_isNedge(myTile).EQ.1
d4399095cb Jean*0099 #else
bd0c662004 Jean*0100        myFace = bi
d4399095cb Jean*0101        southWestCorner = .TRUE.
                0102        southEastCorner = .TRUE.
                0103        northWestCorner = .TRUE.
                0104        northEastCorner = .TRUE.
                0105 #endif /* ALLOW_EXCH2 */
                0106 
                0107        IF ( southWestCorner ) THEN
                0108 C               U(0,1)     D(0,1)      U(1,1)     TILE
                0109 C                |                      |
                0110 C   V(-1,1) --- Z(0,1) --- V(0,1) ---  Z(1,1) --- V(1,1) ---
                0111 C                |                      |
                0112 C               U(0,0)     D(0,0)      U(1,0)     D(1,0)
                0113 C                |                      |
                0114 C                      --- V(0,0) ---  Z(1,0) --- V(1,0) ---
                0115 C                                       |
                0116 C                                      U(1,-1)
d88f4c0b16 Jean*0117          I=1
                0118          J=1
bd0c662004 Jean*0119 C-    to get the same truncation, independent from the face Nb,
                0120 C     do (1+2)+3, and always in the same order (exch3 convention order):
d88f4c0b16 Jean*0121          vort3(I,J)=
d4399095cb Jean*0122      &     +recip_rA(I,J,bi,bj)/AZcorner*(
bd0c662004 Jean*0123      &      ( vFld(I,J)*dyC(I,J,bi,bj)
                0124      &       -uFld(I,J)*dxC(I,J,bi,bj) )
                0125      &      + uFld(I,J-1)*dxC(I,J-1,bi,bj)
d88f4c0b16 Jean*0126      &     )
bd0c662004 Jean*0127 C-    the quick way, but do not get the same truncation on the 3 faces:
                0128 c        vort3(I,J)=
                0129 c    &     +recip_rA(I,J,bi,bj)/AZcorner*(
                0130 c    &      vFld(I,J)*dyC(I,J,bi,bj)
                0131 c    &     -uFld(I,J)*dxC(I,J,bi,bj)
                0132 c    &     +uFld(I,J-1)*dxC(I,J-1,bi,bj)
                0133 c    &     )
                0134 c        IF (hFacZ(i,j).EQ.0. _d 0) vort3(i,j) = 0. _d 0
d4399095cb Jean*0135 #ifdef CALC_CS_CORNER_EXTENDED
bd0c662004 Jean*0136          vort3(I-1,J)=
d4399095cb Jean*0137      &      recip_rAz(I-1,J,bi,bj)*(
bd0c662004 Jean*0138      &      vFld(I-1,J)*dyC(I-1,J,bi,bj)
                0139      &     -vFld(I-2,J)*dyC(I-2,J,bi,bj)
                0140      &     -uFld(I-1,J)*dxC(I-1,J,bi,bj)
                0141      &     +vFld(I+0,J-1)*dyC(I+0,J-1,bi,bj)
d4399095cb Jean*0142      &                             )
                0143 c    &     *maskS(i-1,j,k,bi,bj)*maskS(i-2,j,k,bi,bj)
                0144 c    &     *maskW(i-1,j,k,bi,bj)*maskS(i,j-1,k,bi,bj)
bd0c662004 Jean*0145          vort3(I,J-1)=vort3(I-1,J)
d4399095cb Jean*0146 #endif
                0147        ENDIF
                0148 
                0149        IF ( southEastCorner ) THEN
                0150 C   TILE       U(N+1,1)     D(N+1,1)      U(N+2,1)
                0151 C               |                          |
                0152 C   V(N,1) --- Z(N+1,1) --- V(N+1,1) ---  Z(N+2,1) --- V(N+3,1) ---
                0153 C               |                          |
                0154 C   D(N,0)     U(N+1,0)     D(N+1,0)      U(N+2,0)
                0155 C               |                          |
                0156 C   V(N,0) --- Z(N+1,0) --- V(N+1,0) ---
                0157 C               |                          |
                0158 C              U(N+1,-1)
d88f4c0b16 Jean*0159          I=sNx+1
                0160          J=1
bd0c662004 Jean*0161 C-    to get the same truncation, independent from the face Nb,
                0162 C      (exch3 convention order):
                0163          IF ( myFace.EQ.2 ) THEN
                0164           vort3(I,J)=
                0165      &     +recip_rA(I-1,J,bi,bj)/AZcorner*(
                0166      &      (-uFld(I,J)*dxC(I,J,bi,bj)
                0167      &       -vFld(I-1,J)*dyC(I-1,J,bi,bj) )
                0168      &      + uFld(I,J-1)*dxC(I,J-1,bi,bj)
                0169      &     )
                0170          ELSEIF ( myFace.EQ.4 ) THEN
                0171           vort3(I,J)=
d4399095cb Jean*0172      &     +recip_rA(I-1,J,bi,bj)/AZcorner*(
bd0c662004 Jean*0173      &      (-vFld(I-1,J)*dyC(I-1,J,bi,bj)
                0174      &       +uFld(I,J-1)*dxC(I,J-1,bi,bj) )
                0175      &      - uFld(I,J)*dxC(I,J,bi,bj)
d88f4c0b16 Jean*0176      &     )
bd0c662004 Jean*0177          ELSE
                0178           vort3(I,J)=
                0179      &     +recip_rA(I-1,J,bi,bj)/AZcorner*(
                0180      &      (+uFld(I,J-1)*dxC(I,J-1,bi,bj)
                0181      &       -uFld(I,J)*dxC(I,J,bi,bj)     )
                0182      &      - vFld(I-1,J)*dyC(I-1,J,bi,bj)
                0183      &     )
                0184          ENDIF
                0185 C-    the quick way, but do not get the same truncation on the 3 faces:
                0186 c        vort3(I,J)=
                0187 c    &     +recip_rA(I-1,J,bi,bj)/AZcorner*(
                0188 c    &     -vFld(I-1,J)*dyC(I-1,J,bi,bj)
                0189 c    &     -uFld(I,J)*dxC(I,J,bi,bj)
                0190 c    &     +uFld(I,J-1)*dxC(I,J-1,bi,bj)
                0191 c    &     )
                0192 c        IF (hFacZ(i,j).EQ.0. _d 0) vort3(i,j) = 0. _d 0
d4399095cb Jean*0193 #ifdef CALC_CS_CORNER_EXTENDED
bd0c662004 Jean*0194          vort3(I+1,J)=
d4399095cb Jean*0195      &      recip_rAz(I+1,J,bi,bj)*(
bd0c662004 Jean*0196      &      vFld(I+1,J)*dyC(I+1,J,bi,bj)
                0197      &     -vFld(I-0,J)*dyC(I-0,J,bi,bj)
                0198      &     -uFld(I+1,J)*dxC(I+1,J,bi,bj)
                0199      &     -vFld(I-1,J-1)*dyC(I-1,J-1,bi,bj)
d4399095cb Jean*0200      &                           )
                0201 c    &     *maskS(i+1,j,k,bi,bj)*maskS(i-0,j,k,bi,bj)
                0202 c    &     *maskW(i+1,j,k,bi,bj)*maskS(i-1,j-1,k,bi,bj)
bd0c662004 Jean*0203          vort3(I,J-1)=vort3(I+1,J)
d4399095cb Jean*0204 #endif
                0205        ENDIF
                0206 
                0207        IF ( northWestCorner ) THEN
                0208 C                                            U(1,N+2)
                0209 C                                             |
                0210 C                          --- V(0,N+1) ---  Z(1,N+2) --- V(1,N+2) ---
                0211 C                  |                          |
                0212 C                 U(0,N+1)     D(0,N+1)      U(1,N+1)     D(1,N+1)
                0213 C                  |                          |
                0214 C   V(-1,N+1) --- Z(0,N+1) --- V(0,N+1) ---  Z(1,N+1) --- V(1,N+1) ---
                0215 C                  |                          |
                0216 C                 U(0,N)       D(0,N)        U(1,N)       TILE
d88f4c0b16 Jean*0217          I=1
                0218          J=sNy+1
bd0c662004 Jean*0219 C-    to get the same truncation, independent from the face Nb,
                0220 C      (exch3 convention order):
                0221          IF ( myFace.EQ.1 ) THEN
                0222           vort3(I,J)=
                0223      &     +recip_rA(I,J-1,bi,bj)/AZcorner*(
                0224      &      (+uFld(I,J-1)*dxC(I,J-1,bi,bj)
                0225      &       +vFld(I,J)*dyC(I,J,bi,bj)     )
                0226      &       -uFld(I,J)*dxC(I,J,bi,bj)
                0227      &     )
                0228          ELSEIF ( myFace.EQ.3 ) THEN
                0229           vort3(I,J)=
                0230      &     +recip_rA(I,J-1,bi,bj)/AZcorner*(
                0231      &      (-uFld(I,J)*dxC(I,J,bi,bj)
                0232      &       +uFld(I,J-1)*dxC(I,J-1,bi,bj) )
                0233      &      + vFld(I,J)*dyC(I,J,bi,bj)
                0234      &     )
                0235          ELSE
                0236           vort3(I,J)=
d4399095cb Jean*0237      &     +recip_rA(I,J-1,bi,bj)/AZcorner*(
bd0c662004 Jean*0238      &      (+vFld(I,J)*dyC(I,J,bi,bj)
                0239      &       -uFld(I,J)*dxC(I,J,bi,bj)     )
                0240      &      + uFld(I,J-1)*dxC(I,J-1,bi,bj)
d88f4c0b16 Jean*0241      &     )
bd0c662004 Jean*0242          ENDIF
                0243 C-    the quick way, but do not get the same truncation on the 3 faces:
                0244 c        vort3(I,J)=
                0245 c    &     +recip_rA(I,J-1,bi,bj)/AZcorner*(
                0246 c    &      vFld(I,J)*dyC(I,J,bi,bj)
                0247 c    &     -uFld(I,J)*dxC(I,J,bi,bj)
                0248 c    &     +uFld(I,J-1)*dxC(I,J-1,bi,bj)
                0249 c    &     )
                0250 c        IF (hFacZ(i,j).EQ.0. _d 0) vort3(i,j) = 0. _d 0
d4399095cb Jean*0251 #ifdef CALC_CS_CORNER_EXTENDED
bd0c662004 Jean*0252          vort3(I-1,J)=
d4399095cb Jean*0253      &      recip_rAz(I-1,J,bi,bj)*(
bd0c662004 Jean*0254      &      vFld(I-1,J)*dyC(I-1,J,bi,bj)
                0255      &     -vFld(I-2,J)*dyC(I-2,J,bi,bj)
                0256      &     +vFld(I-0,J+1)*dyC(I-0,J+1,bi,bj)
                0257      &     +uFld(I-1,J-1)*dxC(I-1,J-1,bi,bj)
d4399095cb Jean*0258      &                           )
                0259 c    &     *maskS(i-1,j,k,bi,bj)*maskS(i-2,j,k,bi,bj)
                0260 c    &     *maskS(i,j+1,k,bi,bj)*maskW(i-1,j-1,k,bi,bj)
bd0c662004 Jean*0261          vort3(I,J+1)=vort3(I-1,J)
d4399095cb Jean*0262 #endif
                0263        ENDIF
                0264 
                0265        IF ( northEastCorner ) THEN
                0266 C                U(N+1,N+2)
                0267 C                 |                              |
                0268 C   V(N,N+2) --- Z(N+1,N+2) --- V(N+1,N+2) ---
                0269 C                 |                              |
                0270 C   D(N,N+1)     U(N+1,N+1)     D(N+1,N+1)      U(N+2,N+1)
                0271 C                 |                              |
                0272 C   V(N,N+1) --- Z(N+1,N+1) --- V(N+1,N+1) ---  Z(N+2,N+1) --- V(N+3,N+1) ---
                0273 C                 |                              |
                0274 C   TILE         U(N+1,N)       D(N+1,N)        U(N+2,N)
d88f4c0b16 Jean*0275          I=sNx+1
                0276          J=sNy+1
bd0c662004 Jean*0277 C-    to get the same truncation, independent from the face Nb:
                0278 C      (exch3 convention order):
                0279          IF ( MOD(myFace,2).EQ.1 ) THEN
                0280           vort3(I,J)=
d4399095cb Jean*0281      &     +recip_rA(I-1,J-1,bi,bj)/AZcorner*(
bd0c662004 Jean*0282      &      (-uFld(I,J)*dxC(I,J,bi,bj)
                0283      &       -vFld(I-1,J)*dyC(I-1,J,bi,bj) )
                0284      &      + uFld(I,J-1)*dxC(I,J-1,bi,bj)
d88f4c0b16 Jean*0285      &     )
bd0c662004 Jean*0286          ELSE
                0287           vort3(I,J)=
                0288      &     +recip_rA(I-1,J-1,bi,bj)/AZcorner*(
                0289      &      (+uFld(I,J-1)*dxC(I,J-1,bi,bj)
                0290      &       -uFld(I,J)*dxC(I,J,bi,bj)     )
                0291      &      - vFld(I-1,J)*dyC(I-1,J,bi,bj)
                0292      &     )
                0293          ENDIF
                0294 C-    the quick way, but do not get the same truncation on the 3 faces:
                0295 c        vort3(I,J)=
                0296 c    &     +recip_rA(I-1,J-1,bi,bj)/AZcorner*(
                0297 c    &     -vFld(I-1,J)*dyC(I-1,J,bi,bj)
                0298 c    &     -uFld(I,J)*dxC(I,J,bi,bj)
                0299 c    &     +uFld(I,J-1)*dxC(I,J-1,bi,bj)
                0300 c    &     )
                0301 c        IF (hFacZ(i,j).EQ.0. _d 0) vort3(i,j) = 0. _d 0
d4399095cb Jean*0302 #ifdef CALC_CS_CORNER_EXTENDED
bd0c662004 Jean*0303          vort3(I+1,J)=
d4399095cb Jean*0304      &      recip_rAz(I+1,J,bi,bj)*(
bd0c662004 Jean*0305      &      vFld(I+1,J)*dyC(I+1,J,bi,bj)
                0306      &     -vFld(I-0,J)*dyC(I-0,J,bi,bj)
                0307      &     -vFld(I-1,J+1)*dyC(I-1,J+1,bi,bj)
                0308      &     +uFld(I+1,J-1)*dxC(I+1,J-1,bi,bj)
d4399095cb Jean*0309      &                           )
                0310 c    &     *maskS(i+1,j,k,bi,bj)*maskS(i-0,j,k,bi,bj)
                0311 c    &     *maskS(i-1,j+1,k,bi,bj)*maskW(i+1,j-1,k,bi,bj)
bd0c662004 Jean*0312          vort3(I,J+1)=vort3(I+1,J)
d4399095cb Jean*0313 #endif
d88f4c0b16 Jean*0314        ENDIF
d4399095cb Jean*0315       ENDIF
d88f4c0b16 Jean*0316 
                0317       RETURN
                0318       END