Back to home page

MITgcm

 
 

    


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

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
b75758c66e Jean*0001 #include "SHAP_FILT_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: SHAP_FILT_UV_S2C
                0005 C     !INTERFACE:
                0006       SUBROUTINE SHAP_FILT_UV_S2C(
                0007      U           uFld, vFld, tmpFldU, tmpFldV,
                0008      I           kSize, myTime, myThid )
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | S/R SHAP_FILT_UV_S2C
                0012 C     | o Applies Shapiro filter to velocity field (u & v).
                0013 C     | o use filtering function "S2" = [1 - (d_xx+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:
                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 
                0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     == Routine arguments
                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_S2C
                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)
                0043       _RL     myTime
                0044       INTEGER myThid
                0045 
                0046 #ifdef ALLOW_SHAP_FILT
                0047 
                0048 C     !LOCAL VARIABLES: 
                0049 C     == Local variables ==
                0050       INTEGER bi,bj,k,i,j,N
                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
                0054       _RL noSlipFact
                0055 CEOP
                0056 
                0057       noSlipFact = Shap_noSlip*2. _d 0
                0058 
                0059       IF (nShapUV.gt.0) THEN
                0060 
                0061         DO bj=myByLo(myThid),myByHi(myThid)
                0062          DO bi=myBxLo(myThid),myBxHi(myThid)
                0063           DO k=1,kSize
                0064            DO j=1-OLy,sNy+OLy
                0065             DO i=1-OLx,sNx+OLx
                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---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0078 
                0079 C    [d_xx+d_yy]^n tmpFld
                0080 
                0081        DO N=1,nShapUV
                0082 
                0083         IF (kSize.EQ.Nr) THEN
                0084           CALL EXCH_UV_XYZ_RL(tmpFldU,tmpFldV,.TRUE.,myThid)
                0085         ELSE
                0086           CALL EXCH_UV_XY_RL(tmpFldU,tmpFldV,.TRUE.,myThid)
                0087         ENDIF
                0088 
                0089 
                0090         DO bj=myByLo(myThid),myByHi(myThid)
                0091          DO bi=myBxLo(myThid),myBxHi(myThid)
                0092           DO k=1,kSize
                0093 
                0094 C    Uxx + Uyy :
                0095            DO j=1,sNy
                0096             DO i=1,sNx+1
                0097              maskZj=_maskS(i-1, j ,k,bi,bj)*_maskS( i , j ,k,bi,bj)
                0098              maskZp=_maskS(i-1,j+1,k,bi,bj)*_maskS( i ,j+1,k,bi,bj)
                0099              tmpGrdU(i,j) = 
                0100      &         tmpFldU(i-1,j,k,bi,bj) + tmpFldU(i+1,j,k,bi,bj)
                0101      &             - 2.*tmpFldU(i,j,k,bi,bj) 
                0102      &       +(tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i, j ,k,bi,bj))*maskZp
                0103      &       -(tmpFldU(i, j ,k,bi,bj)-tmpFldU(i,j-1,k,bi,bj))*maskZj
                0104      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
                0105             ENDDO
                0106            ENDDO
                0107 
                0108            IF (useCubedSphereExchange) THEN
                0109             DO i=1,sNx+1,sNx
                0110             j=1
                0111              maskZj=_maskS(i-1, j ,k,bi,bj)*_maskS( i , j ,k,bi,bj)
                0112              maskZp=_maskS(i-1,j+1,k,bi,bj)*_maskS( i ,j+1,k,bi,bj)
                0113              tmpGrdU(i,j) = 
                0114      &         tmpFldU(i-1,j,k,bi,bj) + tmpFldU(i+1,j,k,bi,bj)
                0115      &             - 2.*tmpFldU(i,j,k,bi,bj) 
                0116      &       +(tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i, j ,k,bi,bj))*maskZp
                0117      &       -(tmpFldU(i,j,k,bi,bj)-0*tmpFldU(i,j-1,k,bi,bj))*maskZj
                0118      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
                0119             j=sNy
                0120              maskZj=_maskS(i-1, j ,k,bi,bj)*_maskS( i , j ,k,bi,bj)
                0121              maskZp=_maskS(i-1,j+1,k,bi,bj)*_maskS( i ,j+1,k,bi,bj)
                0122              tmpGrdU(i,j) = 
                0123      &         tmpFldU(i-1,j,k,bi,bj) + tmpFldU(i+1,j,k,bi,bj)
                0124      &             - 2.*tmpFldU(i,j,k,bi,bj) 
                0125      &       +(0*tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i,j,k,bi,bj))*maskZp
                0126      &       -(tmpFldU(i, j ,k,bi,bj)-tmpFldU(i,j-1,k,bi,bj))*maskZj
                0127      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
                0128             ENDDO
                0129            ENDIF
                0130 
                0131            DO j=1,sNy
                0132             DO i=1,sNx+1
                0133              tmpFldU(i,j,k,bi,bj) = 
                0134      &        -0.125 _d 0*tmpGrdU(i,j)*_maskW(i,j,k,bi,bj)
                0135             ENDDO
                0136            ENDDO
                0137 
                0138 C    Vyy + Vxx :
                0139            DO j=1,sNy+1
                0140             DO i=1,sNx
                0141              maskZj=_maskW( i ,j-1,k,bi,bj)*_maskW( i , j ,k,bi,bj)
                0142              maskZp=_maskW(i+1,j-1,k,bi,bj)*_maskW(i+1, j ,k,bi,bj)
                0143              tmpGrdV(i,j) =
                0144      &         tmpFldV(i,j-1,k,bi,bj) + tmpFldV(i,j+1,k,bi,bj)
                0145      &             - 2.*tmpFldV(i,j,k,bi,bj)
                0146      &       +(tmpFldV(i+1,j,k,bi,bj)-tmpFldV( i ,j,k,bi,bj))*maskZp
                0147      &       -(tmpFldV( i ,j,k,bi,bj)-tmpFldV(i-1,j,k,bi,bj))*maskZj
                0148      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
                0149             ENDDO
                0150            ENDDO
                0151 
                0152            IF (useCubedSphereExchange) THEN
                0153             DO j=1,sNy+1,sNy
                0154             i=1
                0155              maskZj=_maskW( i ,j-1,k,bi,bj)*_maskW( i , j ,k,bi,bj)
                0156              maskZp=_maskW(i+1,j-1,k,bi,bj)*_maskW(i+1, j ,k,bi,bj)
                0157              tmpGrdV(i,j) =
                0158      &         tmpFldV(i,j-1,k,bi,bj) + tmpFldV(i,j+1,k,bi,bj)
                0159      &             - 2.*tmpFldV(i,j,k,bi,bj)
                0160      &       +(tmpFldV(i+1,j,k,bi,bj)-tmpFldV( i ,j,k,bi,bj))*maskZp
                0161      &       -(tmpFldV(i,j,k,bi,bj)-0*tmpFldV(i-1,j,k,bi,bj))*maskZj
                0162      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
                0163             i=sNx
                0164              maskZj=_maskW( i ,j-1,k,bi,bj)*_maskW( i , j ,k,bi,bj)
                0165              maskZp=_maskW(i+1,j-1,k,bi,bj)*_maskW(i+1, j ,k,bi,bj)
                0166              tmpGrdV(i,j) =
                0167      &         tmpFldV(i,j-1,k,bi,bj) + tmpFldV(i,j+1,k,bi,bj)
                0168      &             - 2.*tmpFldV(i,j,k,bi,bj)
                0169      &       +(0*tmpFldV(i+1,j,k,bi,bj)-tmpFldV(i,j,k,bi,bj))*maskZp
                0170      &       -(tmpFldV( i ,j,k,bi,bj)-tmpFldV(i-1,j,k,bi,bj))*maskZj
                0171      &       -noSlipFact*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
                0172             ENDDO
                0173            ENDIF
                0174 
                0175            DO j=1,sNy+1
                0176             DO i=1,sNx
                0177              tmpFldV(i,j,k,bi,bj) =
                0178      &        -0.125 _d 0*tmpGrdV(i,j)*_maskS(i,j,k,bi,bj)
                0179             ENDDO
                0180            ENDDO
                0181 
                0182           ENDDO
                0183          ENDDO
                0184         ENDDO
                0185 
                0186        ENDDO
                0187 
                0188 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0189 
                0190 C      F <-  [1 - (d_xx+d_yy)^n *deltaT/tau].F
                0191 
                0192        DO bj=myByLo(myThid),myByHi(myThid)
                0193         DO bi=myBxLo(myThid),myBxHi(myThid)
                0194          DO k=1,kSize
                0195           DO j=1,sNy
                0196            DO i=1,sNx+1
                0197             uFld(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
6824f48020 Jean*0198      &                  -tmpFldU(i,j,k,bi,bj)*deltaTmom/Shap_uvtau
                0199             tmpFldU(i,j,k,bi,bj)= -tmpFldU(i,j,k,bi,bj)/Shap_uvtau
b75758c66e Jean*0200            ENDDO
                0201           ENDDO
                0202           DO j=1,sNy+1
                0203            DO i=1,sNx
                0204             vFld(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
6824f48020 Jean*0205      &                  -tmpFldV(i,j,k,bi,bj)*deltaTmom/Shap_uvtau
                0206             tmpFldV(i,j,k,bi,bj)= -tmpFldV(i,j,k,bi,bj)/Shap_uvtau
b75758c66e Jean*0207            ENDDO
                0208           ENDDO
                0209          ENDDO
                0210         ENDDO
                0211        ENDDO
                0212 
                0213         IF (kSize.EQ.Nr) THEN
                0214           CALL EXCH_UV_XYZ_RL(uFld,vFld,.TRUE.,myThid)
                0215         ELSEIF (kSize.EQ.1) THEN
                0216           CALL EXCH_UV_XY_RL(uFld,vFld,.TRUE.,myThid)
                0217         ELSE
                0218           STOP 'S/R SHAP_FILT_UV_S2C: kSize is wrong'
                0219         ENDIF    
                0220 
                0221       ENDIF
                0222 #endif /* ALLOW_SHAP_FILT */
                0223 
                0224       RETURN
                0225       END