Back to home page

MITgcm

 
 

    


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

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
aea29c8517 Alis*0001 #include "SHAP_FILT_OPTIONS.h"
                0002 
15c70d7cd1 Jean*0003 CBOP
                0004 C     !ROUTINE: SHAP_FILT_UV_S4
                0005 C     !INTERFACE:
                0006       SUBROUTINE SHAP_FILT_UV_S4(
                0007      U           uFld, vFld, tmpFldU, tmpFldV,
                0008      I           kSize, myTime, myThid )
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | S/R SHAP_FILT_UV_S4
                0012 C     | o Applies Shapiro filter to velocity field (u & v).
                0013 C     | o use filtering function "S4" = [1 - d_xx^n][1- d_yy^n]
                0014 C     |     with no grid spacing (computational Filter) ;
                0015 C     |     include No-Slip option
                0016 C     *==========================================================*
                0017 C     \ev
                0018 
                0019 C     !USES:
aea29c8517 Alis*0020       IMPLICIT NONE
                0021 
                0022 C     == Global variables ===
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "SHAP_FILT.h"
                0028 
15c70d7cd1 Jean*0029 C     !INPUT/OUTPUT PARAMETERS:
aea29c8517 Alis*0030 C     == Routine arguments
15c70d7cd1 Jean*0031 C     uFld :: velocity field (U component) on which filter applies
                0032 C     vFld :: velocity field (V component) on which filter applies
                0033 C     tmpFldU :: working temporary array
                0034 C     tmpFldV :: working temporary array
                0035 C     kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
                0036 C     myTime :: Current time in simulation
                0037 C     myThid :: Thread number for this instance of SHAP_FILT_UV_S4
                0038       INTEGER kSize
                0039       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
                0040       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
                0041       _RL tmpFldU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
                0042       _RL tmpFldV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
aea29c8517 Alis*0043       _RL     myTime
                0044       INTEGER myThid
                0045 
                0046 #ifdef ALLOW_SHAP_FILT
                0047 
15c70d7cd1 Jean*0048 C     !LOCAL VARIABLES: 
aea29c8517 Alis*0049 C     == Local variables ==
15c70d7cd1 Jean*0050       INTEGER bi,bj,k,i,j,N
aea29c8517 Alis*0051       _RL tmpGrdU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0052       _RL tmpGrdV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0053       _RL maskZj,maskZp
15c70d7cd1 Jean*0054       _RL noSlipFact
                0055 CEOP
                0056 
                0057       noSlipFact = Shap_noSlip*2. _d 0
aea29c8517 Alis*0058 
                0059       IF (nShapUV.gt.0) THEN
                0060 
                0061         DO bj=myByLo(myThid),myByHi(myThid)
                0062          DO bi=myBxLo(myThid),myBxHi(myThid)
15c70d7cd1 Jean*0063           DO K=1,kSize
                0064            DO J=1-OLy,sNy+OLy
                0065             DO I=1-OLx,sNx+OLx
aea29c8517 Alis*0066              tmpFldU(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
                0067      &                *_maskW(i,j,k,bi,bj)
                0068              tmpFldV(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
                0069      &                *_maskS(i,j,k,bi,bj)
                0070             ENDDO
                0071            ENDDO
                0072           ENDDO
                0073          ENDDO
                0074         ENDDO
                0075 
                0076 
                0077 C      d_xx^n tmpFld 
                0078 
                0079        DO N=1,nShapUV
                0080 
15c70d7cd1 Jean*0081         IF (kSize.EQ.Nr) THEN
                0082           CALL EXCH_UV_XYZ_RL(tmpFldU,tmpFldV,.TRUE.,myThid)
                0083         ELSE
                0084           CALL EXCH_UV_XY_RL(tmpFldU,tmpFldV,.TRUE.,myThid)
                0085         ENDIF
aea29c8517 Alis*0086 
                0087         DO bj=myByLo(myThid),myByHi(myThid)
                0088          DO bi=myBxLo(myThid),myBxHi(myThid)
15c70d7cd1 Jean*0089           DO K=1,kSize
aea29c8517 Alis*0090 
                0091 C          Uxx
                0092            DO J=1,sNy
                0093             DO I=1,sNx+1
                0094              tmpGrdU(i,j) = -0.25*(
                0095      &          tmpFldU(i-1,j,k,bi,bj) + tmpFldU(i+1,j,k,bi,bj)
                0096      &             - 2.*tmpFldU(i,j,k,bi,bj)
                0097      &            )*_maskW(i,j,k,bi,bj)
                0098             ENDDO
                0099            ENDDO
                0100 
                0101            DO J=1,sNy
                0102             DO I=1,sNx+1
                0103              tmpFldU(i,j,k,bi,bj) = tmpGrdU(i,j)
                0104             ENDDO
                0105            ENDDO
                0106 
                0107 C          Vyy
                0108            DO J=1,sNy+1
                0109             DO I=1,sNx
                0110              tmpGrdV(i,j) = -0.25*(
                0111      &          tmpFldV(i,j-1,k,bi,bj) + tmpFldV(i,j+1,k,bi,bj)
                0112      &             - 2.*tmpFldV(i,j,k,bi,bj)
                0113      &            )*_maskS(i,j,k,bi,bj)
                0114             ENDDO
                0115            ENDDO
                0116 
                0117            DO J=1,sNy+1
                0118             DO I=1,sNx
                0119              tmpFldV(i,j,k,bi,bj) = tmpGrdV(i,j)
                0120             ENDDO
                0121            ENDDO
                0122 
                0123           ENDDO
                0124          ENDDO
                0125         ENDDO
                0126 
                0127        ENDDO
                0128 
15c70d7cd1 Jean*0129 C      F <-  [1 - d_xx^n *deltaT/tau].F
aea29c8517 Alis*0130        DO bj=myByLo(myThid),myByHi(myThid)
                0131         DO bi=myBxLo(myThid),myBxHi(myThid)
15c70d7cd1 Jean*0132          DO K=1,kSize
aea29c8517 Alis*0133           DO J=1,sNy
                0134            DO I=1,sNx+1
15c70d7cd1 Jean*0135             uFld(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
                0136      &             -tmpFldU(i,j,k,bi,bj)*deltaTmom/Shap_uvtau
                0137             tmpFldU(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
aea29c8517 Alis*0138            ENDDO
                0139           ENDDO
                0140           DO J=1,sNy+1
                0141            DO I=1,sNx
15c70d7cd1 Jean*0142             vFld(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
                0143      &             -tmpFldV(i,j,k,bi,bj)*deltaTmom/Shap_uvtau
                0144             tmpFldV(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
aea29c8517 Alis*0145            ENDDO
                0146           ENDDO
                0147          ENDDO
                0148         ENDDO
                0149        ENDDO
                0150 
                0151 
                0152 C      d_yy^n tmpFld 
                0153 
                0154        DO N=1,nShapUV
                0155 
15c70d7cd1 Jean*0156         IF (kSize.EQ.Nr) THEN
                0157           CALL EXCH_UV_XYZ_RL(tmpFldU,tmpFldV,.TRUE.,myThid)
                0158         ELSE
                0159           CALL EXCH_UV_XY_RL(tmpFldU,tmpFldV,.TRUE.,myThid)
                0160         ENDIF
aea29c8517 Alis*0161 
                0162         DO bj=myByLo(myThid),myByHi(myThid)
                0163          DO bi=myBxLo(myThid),myBxHi(myThid)
15c70d7cd1 Jean*0164           DO K=1,kSize
aea29c8517 Alis*0165 
                0166 C          Uyy
                0167            DO J=1,sNy
                0168             DO I=1,sNx+1
                0169              maskZj=_maskS(i-1, j ,k,bi,bj)
                0170      &             *_maskS( i , j ,k,bi,bj)
                0171              maskZp=_maskS(i-1,j+1,k,bi,bj)
                0172      &             *_maskS( i ,j+1,k,bi,bj)
                0173              tmpGrdU(i,j) = -0.25*(
                0174      &        (tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i, j ,k,bi,bj))*maskZp
                0175      &       -(tmpFldU(i, j ,k,bi,bj)-tmpFldU(i,j-1,k,bi,bj))*maskZj
15c70d7cd1 Jean*0176      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
aea29c8517 Alis*0177      &             )*_maskW(i,j,k,bi,bj)
                0178             ENDDO
                0179            ENDDO
                0180 
                0181            IF (useCubedSphereExchange) THEN
                0182             J=1
                0183             DO I=1,sNx+1,sNx
                0184              maskZj=maskS(i-1, j ,k,bi,bj)*maskS( i , j ,k,bi,bj)
                0185              maskZp=maskS(i-1,j+1,k,bi,bj)*maskS( i ,j+1,k,bi,bj)
                0186              tmpGrdU(i,j) = -0.25*(
                0187      &        (tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i, j ,k,bi,bj))*maskZp
                0188      &       -(tmpFldU(i, j ,k,bi,bj)-0*tmpFldU(i,j-1,k,bi,bj))*maskZj
15c70d7cd1 Jean*0189      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
aea29c8517 Alis*0190      &             )*_maskW(i,j,k,bi,bj)
                0191             ENDDO
                0192             J=sNy
                0193             DO I=1,sNx+1,sNx
                0194              maskZj=maskS(i-1, j ,k,bi,bj)*maskS( i , j ,k,bi,bj)
                0195              maskZp=maskS(i-1,j+1,k,bi,bj)*maskS( i ,j+1,k,bi,bj)
                0196              tmpGrdU(i,j) = -0.25*(
                0197      &        (0*tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i, j ,k,bi,bj))*maskZp
                0198      &       -(tmpFldU(i, j ,k,bi,bj)-tmpFldU(i,j-1,k,bi,bj))*maskZj
15c70d7cd1 Jean*0199      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
aea29c8517 Alis*0200      &             )*_maskW(i,j,k,bi,bj)
                0201             ENDDO
                0202            ENDIF
                0203 
                0204            DO J=1,sNy
                0205             DO I=1,sNx+1
                0206              tmpFldU(i,j,k,bi,bj) = tmpGrdU(i,j)
                0207             ENDDO
                0208            ENDDO
                0209 
                0210 C          Vxx
                0211            DO J=1,sNy+1
                0212             DO I=1,sNx
                0213              maskZj=_maskW( i ,j-1,k,bi,bj)
                0214      &             *_maskW( i , j ,k,bi,bj)
                0215              maskZp=_maskW(i+1,j-1,k,bi,bj)
                0216      &             *_maskW(i+1, j ,k,bi,bj)
                0217              tmpGrdV(i,j) = -0.25*(
                0218      &        (tmpFldV(i+1,j,k,bi,bj)-tmpFldV( i ,j,k,bi,bj))*maskZp
                0219      &       -(tmpFldV( i ,j,k,bi,bj)-tmpFldV(i-1,j,k,bi,bj))*maskZj
15c70d7cd1 Jean*0220      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
aea29c8517 Alis*0221      &             )*_maskS(i,j,k,bi,bj)
                0222             ENDDO
                0223            ENDDO
                0224 
                0225            IF (useCubedSphereExchange) THEN
                0226             DO J=1,sNy+1,sNy
                0227             I=1
                0228              maskZj=maskW( i ,j-1,k,bi,bj)*maskW( i , j ,k,bi,bj)
                0229              maskZp=maskW(i+1,j-1,k,bi,bj)*maskW(i+1, j ,k,bi,bj)
                0230              tmpGrdV(i,j) = -0.25*(
                0231      &        (tmpFldV(i+1,j,k,bi,bj)-tmpFldV( i ,j,k,bi,bj))*maskZp
                0232      &       -(tmpFldV( i ,j,k,bi,bj)-0*tmpFldV(i-1,j,k,bi,bj))*maskZj
15c70d7cd1 Jean*0233      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
aea29c8517 Alis*0234      &       -2.*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
                0235      &             )*_maskS(i,j,k,bi,bj)
                0236             ENDDO
                0237             DO J=1,sNy+1,sNy
                0238             I=sNx
                0239              maskZj=maskW( i ,j-1,k,bi,bj)*maskW( i , j ,k,bi,bj)
                0240              maskZp=maskW(i+1,j-1,k,bi,bj)*maskW(i+1, j ,k,bi,bj)
                0241              tmpGrdV(i,j) = -0.25*(
                0242      &        (0*tmpFldV(i+1,j,k,bi,bj)-tmpFldV( i ,j,k,bi,bj))*maskZp
                0243      &       -(tmpFldV( i ,j,k,bi,bj)-tmpFldV(i-1,j,k,bi,bj))*maskZj
15c70d7cd1 Jean*0244      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
aea29c8517 Alis*0245      &             )*_maskS(i,j,k,bi,bj)
                0246             ENDDO
                0247            ENDIF
                0248 
                0249            DO J=1,sNy+1
                0250             DO I=1,sNx
                0251              tmpFldV(i,j,k,bi,bj) = tmpGrdV(i,j)
                0252             ENDDO
                0253            ENDDO
                0254 
                0255           ENDDO
                0256          ENDDO
                0257         ENDDO
                0258 
                0259        ENDDO
                0260 
15c70d7cd1 Jean*0261 C      F <-  [1 - d_yy^n *deltaT/tau].F
aea29c8517 Alis*0262        DO bj=myByLo(myThid),myByHi(myThid)
                0263         DO bi=myBxLo(myThid),myBxHi(myThid)
15c70d7cd1 Jean*0264          DO K=1,kSize
aea29c8517 Alis*0265           DO J=1,sNy
                0266            DO I=1,sNx+1
15c70d7cd1 Jean*0267             uFld(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
                0268      &             -tmpFldU(i,j,k,bi,bj)*deltaTmom/Shap_uvtau
aea29c8517 Alis*0269            ENDDO
                0270           ENDDO
                0271           DO J=1,sNy+1
                0272            DO I=1,sNx
15c70d7cd1 Jean*0273             vFld(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
                0274      &             -tmpFldV(i,j,k,bi,bj)*deltaTmom/Shap_uvtau
aea29c8517 Alis*0275            ENDDO
                0276           ENDDO
                0277          ENDDO
                0278         ENDDO
                0279        ENDDO
                0280 
15c70d7cd1 Jean*0281         IF (kSize.EQ.Nr) THEN
                0282           CALL EXCH_UV_XYZ_RL(uFld,vFld,.TRUE.,myThid)
                0283         ELSEIF (kSize.EQ.1) THEN
                0284           CALL EXCH_UV_XY_RL(uFld,vFld,.TRUE.,myThid)
                0285         ELSE
                0286           STOP 'S/R SHAP_FILT_UV_S4: kSize is wrong'
                0287         ENDIF
                0288 
aea29c8517 Alis*0289       ENDIF
                0290 #endif /* ALLOW_SHAP_FILT */
                0291 
                0292       RETURN
                0293       END