pkg/shap_filt/shap_filt_uv_s2.F

aea29c8517 Alis*0001 #include "SHAP_FILT_OPTIONS.h"
113b21cd45 Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
aea29c8517 Alis*0005
15c70d7cd1 Jean*0006 CBOP
                0007 C     !ROUTINE: SHAP_FILT_UV_S2
                0008 C     !INTERFACE:
                0009       SUBROUTINE SHAP_FILT_UV_S2(
                0010      U           uFld, vFld, tmpFldU, tmpFldV,
                0011      I           kSize, myTime, myThid )
                0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
                0014 C     | S/R SHAP_FILT_UV_S2
                0015 C     | o Applies Shapiro filter to velocity field (u & v).
                0016 C     | o use filtering function "S2" = [1 - (d_xx+d_yy)^n]
                0017 C     | o Options for computational filter (no grid spacing)
                0018 C     |   or physical space filter (with grid spacing) or both.
                0019 C     *==========================================================*
                0020 C     \ev
                0021
                0022 C     !USES:
aea29c8517 Alis*0023       IMPLICIT NONE
                0024
                0025 C     == Global variables ===
                0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
                0029 #include "GRID.h"
                0030 #include "SHAP_FILT.h"
c424ee7cc7 Jean*0031 c#ifdef ALLOW_EXCH2
                0032 c#include "W2_EXCH2_SIZE.h"
                0033 c#include "W2_EXCH2_TOPOLOGY.h"
                0034 c#endif /* ALLOW_EXCH2 */
aea29c8517 Alis*0035
15c70d7cd1 Jean*0036 C     !INPUT/OUTPUT PARAMETERS:
aea29c8517 Alis*0037 C     == Routine arguments
15c70d7cd1 Jean*0038 C     uFld :: velocity field (U component) on which filter applies
                0039 C     vFld :: velocity field (V component) on which filter applies
                0040 C     tmpFldU :: working temporary array
                0041 C     tmpFldV :: working temporary array
                0042 C     kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
                0043 C     myTime :: Current time in simulation
                0044 C     myThid :: Thread number for this instance of SHAP_FILT_UV_S2
                0045       INTEGER kSize
                0046       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
                0047       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
8edc6fec00 Jean*0048       _RL tmpFldU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
                0049       _RL tmpFldV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
aea29c8517 Alis*0050       _RL     myTime
                0051       INTEGER myThid
                0052
                0053 #ifdef ALLOW_SHAP_FILT
15c70d7cd1 Jean*0054
e64f43ca56 Jean*0055 C------
1891130b05 Jean*0056 C  Combine computational Filter of Div & Vorticity
                0057 C   and Physical Filter of U,V field
e64f43ca56 Jean*0058 C   nShapUVPhys = 0  ==> use only computational Filter
1891130b05 Jean*0059 C   nShapUVPhys = 1  ==> compute Div & Vort. with  Grid factors,
e64f43ca56 Jean*0060 C                        Filter Div & Vort. Numerically (power nShapUV-1)
1891130b05 Jean*0061 C                        and return filtered U.V in physical space
e64f43ca56 Jean*0062 C   nShapUVPhys = nShapUV  ==> Filter in Physical space only (power nShapUV)
1891130b05 Jean*0063 C------
aea29c8517 Alis*0064
1891130b05 Jean*0065 C     !LOCAL VARIABLES:
aea29c8517 Alis*0066 C     == Local variables ==
8edc6fec00 Jean*0067       INTEGER bi,bj,k,i,j,n
aea29c8517 Alis*0068       _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8edc6fec00 Jean*0069       _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
aea29c8517 Alis*0070       _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0071       _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113b21cd45 Mart*0072 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0073       INTEGER kkey
113b21cd45 Mart*0074 #endif
15c70d7cd1 Jean*0075 CEOP
aea29c8517 Alis*0076
113b21cd45 Mart*0077 #ifdef ALLOW_AUTODIFF_TAMC
                0078 CADJ INIT loctape_shapfilt_bibj_k = COMMON, nSx*nSy*Nr
                0079 #endif
e64f43ca56 Jean*0080       IF (nShapUV.GT.0 .AND. Shap_uvtau.GT.0.) THEN
aea29c8517 Alis*0081
e63e6d8a75 Jean*0082         IF (useCubedSphereExchange) THEN
                0083 C-      need to initialise hDiv for MOM_VI_DEL2UV(call FILL_CS_CORNER_TR_RL)
41b8eb6e45 Jean*0084           DO j=1-OLy,sNy+OLy
                0085            DO i=1-OLx,sNx+OLx
e63e6d8a75 Jean*0086             hDiv(i,j) = 0. _d 0
                0087            ENDDO
                0088           ENDDO
                0089         ENDIF
                0090
aea29c8517 Alis*0091         DO bj=myByLo(myThid),myByHi(myThid)
                0092          DO bi=myBxLo(myThid),myBxHi(myThid)
d4399095cb Jean*0093           DO k=1,kSize
41b8eb6e45 Jean*0094            DO j=1-OLy,sNy+OLy
                0095             DO i=1-OLx,sNx+OLx
aea29c8517 Alis*0096              tmpFldU(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
                0097      &                *_maskW(i,j,k,bi,bj)
                0098              tmpFldV(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
                0099      &                *_maskS(i,j,k,bi,bj)
                0100             ENDDO
                0101            ENDDO
                0102           ENDDO
                0103          ENDDO
                0104         ENDDO
                0105
e64f43ca56 Jean*0106 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
aea29c8517 Alis*0107
d4399095cb Jean*0108 C--   [d_xx+d_yy]^n tmpFld
aea29c8517 Alis*0109
8edc6fec00 Jean*0110        DO n=1,nShapUV
aea29c8517 Alis*0111
8edc6fec00 Jean*0112         IF ( MOD(n,2).EQ.1 .OR. Shap_alwaysExchUV ) THEN
e63e6d8a75 Jean*0113           CALL EXCH_UV_3D_RL( tmpFldU,tmpFldV, .TRUE., kSize, myThid )
15c70d7cd1 Jean*0114         ENDIF
aea29c8517 Alis*0115
                0116         DO bj=myByLo(myThid),myByHi(myThid)
                0117          DO bi=myBxLo(myThid),myBxHi(myThid)
d4399095cb Jean*0118           DO k=1,kSize
aea29c8517 Alis*0119
113b21cd45 Mart*0120 #ifdef ALLOW_AUTODIFF_TAMC
                0121 C     This store uses a wrongly-sized tape and stores the same field
                0122 C     nShapUV times, but since this field does not change here, we
                0123 C     can do it like this and suppress the recomputation warning
20dee61641 Mart*0124            kkey = k + (bi-1 + (bj-1)*nSx) * Nr
113b21cd45 Mart*0125 CADJ STORE hFacZ = loctape_shapfilt_bibj_k, key = kkey
                0126 #endif
8edc6fec00 Jean*0127            IF ( n.LE.nShapUVPhys .OR.
                0128      &          n.GT.nShapUV-nShapUVPhys )
                0129      &     CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
                0130
d4399095cb Jean*0131 C-    [d_xx+d_yy] tmpFld
8edc6fec00 Jean*0132          IF (n.LE.nShapUVPhys) THEN
999b38430d Alis*0133            CALL MOM_CALC_HDIV(bi,bj,k,2,
e64f43ca56 Jean*0134      I                    tmpFldU(1-OLx,1-OLy,k,bi,bj),
                0135      I                    tmpFldV(1-OLx,1-OLy,k,bi,bj),
                0136      &                    hDiv,myThid)
d4399095cb Jean*0137 #ifdef USE_SHAP_CALC_VORTICITY
d88f4c0b16 Jean*0138            CALL SHAP_FILT_RELVORT3(bi,bj,k,
e64f43ca56 Jean*0139      I                    tmpFldU(1-OLx,1-OLy,k,bi,bj),
                0140      I                    tmpFldV(1-OLx,1-OLy,k,bi,bj),
                0141      &                    hFacZ,vort3,myThid)
d4399095cb Jean*0142 #else
                0143            CALL MOM_CALC_RELVORT3(bi,bj,k,
                0144      I                    tmpFldU(1-OLx,1-OLy,k,bi,bj),
                0145      I                    tmpFldV(1-OLx,1-OLy,k,bi,bj),
                0146      &                    hFacZ,vort3,myThid)
                0147 #endif
e64f43ca56 Jean*0148          ELSE
d4399095cb Jean*0149 C-    replace Physical calc Div & Vort by computational one :
41b8eb6e45 Jean*0150            DO j=1-OLy,sNy+OLy-1
                0151             DO i=1-OLx,sNx+OLx-1
8edc6fec00 Jean*0152              hDiv(i,j)=(tmpFldU(i+1,j,k,bi,bj)-tmpFldU(i,j,k,bi,bj))
                0153      &                +(tmpFldV(i,j+1,k,bi,bj)-tmpFldV(i,j,k,bi,bj))
aea29c8517 Alis*0154             ENDDO
                0155            ENDDO
8edc6fec00 Jean*0156            CALL SHAP_FILT_COMPUTVORT(
                0157      I                    tmpFldU(1-OLx,1-OLy,k,bi,bj),
                0158      I                    tmpFldV(1-OLx,1-OLy,k,bi,bj),
                0159      O                    vort3,
                0160      I                    k,bi,bj,myThid)
e64f43ca56 Jean*0161          ENDIF
                0162
                0163 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0164
8edc6fec00 Jean*0165          IF (n.GT.nShapUV-nShapUVPhys) THEN
ae409e69d3 Jean*0166            IF (Shap_uvLength.LT.0.) THEN
41b8eb6e45 Jean*0167              DO j=1-OLy,sNy+OLy-1
                0168               DO i=1-OLx,sNx+OLx-1
ae409e69d3 Jean*0169                 hDiv(i,j) = hDiv(i,j) * rA(i,j,bi,bj)
e63e6d8a75 Jean*0170               ENDDO
                0171              ENDDO
41b8eb6e45 Jean*0172              DO j=2-OLy,sNy+OLy
                0173               DO i=2-OLx,sNx+OLx
ae409e69d3 Jean*0174                 vort3(i,j)= vort3(i,j)*rAz(i,j,bi,bj)
                0175               ENDDO
                0176              ENDDO
                0177            ENDIF
e64f43ca56 Jean*0178            CALL MOM_VI_DEL2UV(
                0179      I                    bi,bj,k,hDiv,vort3,hFacZ,
                0180      O                    tmpFldU(1-OLx,1-OLy,k,bi,bj),
                0181      O                    tmpFldV(1-OLx,1-OLy,k,bi,bj),
8edc6fec00 Jean*0182      I                    myThid)
ae409e69d3 Jean*0183            IF (Shap_uvLength.LT.0.) THEN
41b8eb6e45 Jean*0184             DO j=2-OLy,sNy+OLy-1
                0185              DO i=2-OLx,sNx+OLx-1
ae409e69d3 Jean*0186               tmpFldU(i,j,k,bi,bj) = -0.125 _d 0*tmpFldU(i,j,k,bi,bj)
                0187      &                                          *_maskW(i,j,k,bi,bj)
                0188               tmpFldV(i,j,k,bi,bj) = -0.125 _d 0*tmpFldV(i,j,k,bi,bj)
                0189      &                                          *_maskS(i,j,k,bi,bj)
                0190              ENDDO
                0191             ENDDO
                0192            ELSEIF (Shap_uvLength.EQ.0.) THEN
41b8eb6e45 Jean*0193             DO j=2-OLy,sNy+OLy-1
                0194              DO i=2-OLx,sNx+OLx-1
d4399095cb Jean*0195               tmpFldU(i,j,k,bi,bj) = -0.125 _d 0*tmpFldU(i,j,k,bi,bj)
                0196      &                            *rAw(i,j,bi,bj)*_maskW(i,j,k,bi,bj)
                0197               tmpFldV(i,j,k,bi,bj) = -0.125 _d 0*tmpFldV(i,j,k,bi,bj)
                0198      &                            *rAs(i,j,bi,bj)*_maskS(i,j,k,bi,bj)
e64f43ca56 Jean*0199              ENDDO
                0200             ENDDO
                0201            ELSE
41b8eb6e45 Jean*0202             DO j=2-OLy,sNy+OLy-1
                0203              DO i=2-OLx,sNx+OLx-1
d4399095cb Jean*0204               tmpFldU(i,j,k,bi,bj) = -0.125 _d 0*tmpFldU(i,j,k,bi,bj)
                0205      &               *Shap_uvLength*Shap_uvLength*_maskW(i,j,k,bi,bj)
                0206               tmpFldV(i,j,k,bi,bj) = -0.125 _d 0*tmpFldV(i,j,k,bi,bj)
                0207      &               *Shap_uvLength*Shap_uvLength*_maskS(i,j,k,bi,bj)
e64f43ca56 Jean*0208              ENDDO
                0209             ENDDO
8edc6fec00 Jean*0210            ENDIF
e64f43ca56 Jean*0211          ELSE
8edc6fec00 Jean*0212 C-       replace derivatives in physical space of Div & Vort by computational ones:
41b8eb6e45 Jean*0213 #ifndef ALLOW_AUTODIFF
8edc6fec00 Jean*0214            IF ( .NOT.Shap_alwaysExchUV
                0215      &          .AND. useCubedSphereExchange ) THEN
                0216 C          to compute d/dx(hDiv), fill corners with appropriate values:
93e3461d85 Jean*0217              CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
1891130b05 Jean*0218      &                                  hDiv, bi,bj, myThid )
8edc6fec00 Jean*0219            ENDIF
                0220 #endif
41b8eb6e45 Jean*0221            DO j=2-OLy,sNy+OLy-1
                0222             DO i=2-OLx,sNx+OLx-1
d4399095cb Jean*0223              tmpFldU(i,j,k,bi,bj) = -0.125 _d 0*
8edc6fec00 Jean*0224      &                   ( (hDiv(i,j) - hDiv(i-1,j))
                0225      &                    -(vort3(i,j+1)-vort3(i,j))
e64f43ca56 Jean*0226      &                   )*maskW(i,j,k,bi,bj)
aea29c8517 Alis*0227             ENDDO
                0228            ENDDO
41b8eb6e45 Jean*0229 #ifndef ALLOW_AUTODIFF
8edc6fec00 Jean*0230            IF ( .NOT.Shap_alwaysExchUV
                0231      &          .AND. useCubedSphereExchange ) THEN
                0232 C          to compute d/dy(hDiv), fill corners with appropriate values:
93e3461d85 Jean*0233              CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
1891130b05 Jean*0234      &                                  hDiv, bi,bj, myThid )
8edc6fec00 Jean*0235            ENDIF
                0236 #endif
41b8eb6e45 Jean*0237            DO j=2-OLy,sNy+OLy-1
                0238             DO i=2-OLx,sNx+OLx-1
d4399095cb Jean*0239              tmpFldV(i,j,k,bi,bj) = -0.125 _d 0*
8edc6fec00 Jean*0240      &                   ( (hDiv(i,j) - hDiv(i,j-1))
                0241      &                    +(vort3(i+1,j)-vort3(i,j))
e64f43ca56 Jean*0242      &                   )*maskS(i,j,k,bi,bj)
aea29c8517 Alis*0243             ENDDO
                0244            ENDDO
                0245
e64f43ca56 Jean*0246          ENDIF
aea29c8517 Alis*0247
d4399095cb Jean*0248 C end loops  k / bi / bj
aea29c8517 Alis*0249           ENDDO
                0250          ENDDO
                0251         ENDDO
8edc6fec00 Jean*0252 C end loop n=1,nShapUV
aea29c8517 Alis*0253        ENDDO
                0254
e64f43ca56 Jean*0255 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0256
15c70d7cd1 Jean*0257 C      F <-  [1 - (d_xx+d_yy)^n *deltaT/tau].F
aea29c8517 Alis*0258        DO bj=myByLo(myThid),myByHi(myThid)
                0259         DO bi=myBxLo(myThid),myBxHi(myThid)
d4399095cb Jean*0260          DO k=1,kSize
8edc6fec00 Jean*0261           DO j=1,sNy
                0262            DO i=1,sNx+1
aea29c8517 Alis*0263             uFld(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
41b8eb6e45 Jean*0264      &                  -tmpFldU(i,j,k,bi,bj)*deltaTMom/Shap_uvtau
6824f48020 Jean*0265             tmpFldU(i,j,k,bi,bj)= -tmpFldU(i,j,k,bi,bj)/Shap_uvtau
e64f43ca56 Jean*0266            ENDDO
                0267           ENDDO
d4399095cb Jean*0268           DO j=1,sNy+1
                0269            DO i=1,sNx
aea29c8517 Alis*0270             vFld(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
41b8eb6e45 Jean*0271      &                  -tmpFldV(i,j,k,bi,bj)*deltaTMom/Shap_uvtau
6824f48020 Jean*0272             tmpFldV(i,j,k,bi,bj)= -tmpFldV(i,j,k,bi,bj)/Shap_uvtau
aea29c8517 Alis*0273            ENDDO
                0274           ENDDO
                0275          ENDDO
                0276         ENDDO
                0277        ENDDO
                0278
e63e6d8a75 Jean*0279        IF ( ( OLx.LE.2 .OR. OLy.LE.2 ) .AND.
                0280      &       MOD(nShapUV,2).EQ.0 .AND. .NOT.Shap_alwaysExchUV )
                0281      &  CALL EXCH_UV_3D_RL( uFld, vFld, .TRUE., kSize, myThid )
aea29c8517 Alis*0282
                0283       ENDIF
                0284 #endif /* ALLOW_SHAP_FILT */
                0285
                0286       RETURN
                0287       END