pkg/shap_filt/shap_filt_tracer_s2.F

aea29c8517 Alis*0001 #include "SHAP_FILT_OPTIONS.h"
ae409e69d3 Jean*0002
15c70d7cd1 Jean*0003 CBOP
                0004 C     !ROUTINE: SHAP_FILT_TRACER_S2
                0005 C     !INTERFACE:
                0006       SUBROUTINE SHAP_FILT_TRACER_S2(
                0007      U           field, tmpFld,
e93c2e7dac Jean*0008      I           nShapTr, exchInOut, kSize,
                0009      I           myTime, myIter, myThid )
15c70d7cd1 Jean*0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
                0012 C     | S/R SHAP_FILT_TRACER_S2
                0013 C     | o Applies Shapiro filter to 2D field (cell center).
                0014 C     | o use filtering function "S2" = [1 - (d_xx+d_yy)^n]
                0015 C     | o Options for computational filter (no grid spacing)
                0016 C     |   or physical space filter (with grid spacing) or both.
                0017 C     *==========================================================*
                0018 C     \ev
ae409e69d3 Jean*0019
15c70d7cd1 Jean*0020 C     !USES:
aea29c8517 Alis*0021       IMPLICIT NONE
ae409e69d3 Jean*0022
aea29c8517 Alis*0023 C     == Global variables ===
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
                0028 #include "SHAP_FILT.h"
                0029
15c70d7cd1 Jean*0030 C     !INPUT/OUTPUT PARAMETERS:
aea29c8517 Alis*0031 C     == Routine arguments
e93c2e7dac Jean*0032 C     field     :: cell-centered 2D field on which filter applies
                0033 C     tmpFld    :: working temporary array
                0034 C     nShapTr   :: (total) power of the filter for this tracer
                0035 C     exchInOut :: apply Exchanges to fill overlap region:
                0036 C            = 0 : do not apply Exch on neither input nor output field
                0037 C            = 1 : apply Exch on input field
                0038 C                  (needed if input field has invalid overlap)
                0039 C            = 2 : apply Exch on output field (after the filter)
                0040 C            = 3 : apply Exch on both input & output field
                0041 C     kSize     :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
                0042 C     myTime    :: Current time in simulation
                0043 C     myIter    :: Current iteration number in simulation
                0044 C     myThid    :: Thread number for this instance of SHAP_FILT_TRACER_S2
8a5457c6b4 Jean*0045       INTEGER kSize
15c70d7cd1 Jean*0046       _RL field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
                0047       _RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
8a5457c6b4 Jean*0048       INTEGER nShapTr, exchInOut
aea29c8517 Alis*0049       _RL     myTime
e93c2e7dac Jean*0050       INTEGER myIter
aea29c8517 Alis*0051       INTEGER myThid
ae409e69d3 Jean*0052
aea29c8517 Alis*0053 #ifdef ALLOW_SHAP_FILT
                0054
15c70d7cd1 Jean*0055 C     !LOCAL VARIABLES:
aea29c8517 Alis*0056 C     == Local variables ==
8edc6fec00 Jean*0057       INTEGER nShapComput
                0058       INTEGER bi,bj,k,i,j,n
aea29c8517 Alis*0059       _RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8edc6fec00 Jean*0060       _RL tmpFdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL tmpFdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
15c70d7cd1 Jean*0062 CEOP
aea29c8517 Alis*0063
e93c2e7dac Jean*0064       IF ( exchInOut.LT.0 .OR. exchInOut.GT.3 ) THEN
                0065          STOP 'S/R SHAP_FILT_TRACER_S2: exchInOut is wrong'
                0066       ENDIF
                0067
ae409e69d3 Jean*0068       IF (nShapTr.GT.0) THEN
e64f43ca56 Jean*0069 C-------
                0070 C  Apply computational filter ^(nShap-nShapPhys) without grid factor
8edc6fec00 Jean*0071 C  then apply Physical filter ^nShapPhys  with grid factors
e64f43ca56 Jean*0072 C-------
8edc6fec00 Jean*0073         nShapComput = nShapTr - nShapTrPhys
aea29c8517 Alis*0074
                0075         DO bj=myByLo(myThid),myByHi(myThid)
                0076          DO bi=myBxLo(myThid),myBxHi(myThid)
8edc6fec00 Jean*0077           DO k=1,kSize
41b8eb6e45 Jean*0078            DO j=1-OLy,sNy+OLy
                0079             DO i=1-OLx,sNx+OLx
aea29c8517 Alis*0080              tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
                0081             ENDDO
                0082            ENDDO
                0083           ENDDO
                0084          ENDDO
                0085         ENDDO
                0086
8edc6fec00 Jean*0087 C      ( d_xx +d_yy )^n tmpFld
aea29c8517 Alis*0088
8edc6fec00 Jean*0089        DO n=1,nShapTr
aea29c8517 Alis*0090
e93c2e7dac Jean*0091         IF ( ( MOD(n,2).EQ.1 .OR. Shap_alwaysExchTr ) .AND.
                0092      &       ( n.GE.2 .OR. MOD(exchInOut,2).EQ.1 )  ) THEN
8edc6fec00 Jean*0093          IF (kSize.EQ.Nr) THEN
7163a40534 Jean*0094           _EXCH_XYZ_RL( tmpFld, myThid )
8edc6fec00 Jean*0095          ELSEIF (kSize.EQ.1) THEN
7163a40534 Jean*0096           _EXCH_XY_RL( tmpFld, myThid )
8edc6fec00 Jean*0097          ELSE
                0098           STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
                0099          ENDIF
15c70d7cd1 Jean*0100         ENDIF
ae409e69d3 Jean*0101
aea29c8517 Alis*0102         DO bj=myByLo(myThid),myByHi(myThid)
                0103          DO bi=myBxLo(myThid),myBxHi(myThid)
8edc6fec00 Jean*0104           DO k=1,kSize
aea29c8517 Alis*0105
8edc6fec00 Jean*0106 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
aea29c8517 Alis*0107
8edc6fec00 Jean*0108 C--        Calculate gradient in X direction:
41b8eb6e45 Jean*0109 #ifndef ALLOW_AUTODIFF
8edc6fec00 Jean*0110            IF ( .NOT.Shap_alwaysExchTr
                0111      &          .AND. useCubedSphereExchange ) THEN
                0112 C          to compute d/dx(tmpFld), fill corners with appropriate values:
93e3461d85 Jean*0113              CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
41b8eb6e45 Jean*0114      &                                 tmpFld(1-OLx,1-OLy,k,bi,bj),
8edc6fec00 Jean*0115      &                                 bi,bj, myThid )
                0116            ENDIF
e13cb171e9 Jean*0117 #endif
8edc6fec00 Jean*0118            IF ( n.LE.nShapComput ) THEN
                0119 C-         Computational space: del_i
                0120              DO j=0,sNy+1
                0121               DO i=0,sNx+2
e93c2e7dac Jean*0122                tmpFdx(i,j) =
8edc6fec00 Jean*0123      &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
                0124      &                     *_maskW(i,j,k,bi,bj)
                0125               ENDDO
                0126              ENDDO
                0127            ELSE
                0128 C-         Physical space: grad_x
                0129              DO j=0,sNy+1
                0130               DO i=0,sNx+2
e93c2e7dac Jean*0131                tmpFdx(i,j) =
8edc6fec00 Jean*0132      &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
                0133      &                     *_hFacW(i,j,k,bi,bj)
                0134      &                     *dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)
                0135               ENDDO
                0136              ENDDO
                0137            ENDIF
e64f43ca56 Jean*0138
8edc6fec00 Jean*0139 C--        Calculate gradient in Y direction:
41b8eb6e45 Jean*0140 #ifndef ALLOW_AUTODIFF
8edc6fec00 Jean*0141            IF ( .NOT.Shap_alwaysExchTr
                0142      &          .AND. useCubedSphereExchange ) THEN
                0143 C          to compute d/dy(tmpFld), fill corners with appropriate values:
93e3461d85 Jean*0144              CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
41b8eb6e45 Jean*0145      &                                 tmpFld(1-OLx,1-OLy,k,bi,bj),
8edc6fec00 Jean*0146      &                                 bi,bj, myThid )
                0147            ENDIF
e13cb171e9 Jean*0148 #endif
8edc6fec00 Jean*0149            IF ( n.LE.nShapComput ) THEN
                0150 C-         Computational space: del_j
                0151              DO j=0,sNy+2
                0152               DO i=0,sNx+1
e93c2e7dac Jean*0153                tmpFdy(i,j) =
8edc6fec00 Jean*0154      &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
                0155      &                     *_maskS(i,j,k,bi,bj)
                0156               ENDDO
                0157              ENDDO
                0158            ELSE
                0159 C-         Physical space: grad_y
                0160              DO j=0,sNy+2
                0161               DO i=0,sNx+1
e93c2e7dac Jean*0162                tmpFdy(i,j) =
8edc6fec00 Jean*0163      &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
                0164      &                     *_hFacS(i,j,k,bi,bj)
                0165      &                     *dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)
                0166               ENDDO
                0167              ENDDO
                0168            ENDIF
aea29c8517 Alis*0169
8edc6fec00 Jean*0170 C--        Calculate (d_xx + d_yy) tmpFld :
                0171            DO j=0,sNy+1
                0172              DO i=0,sNx+1
                0173                tmpGrd(i,j) = ( tmpFdx(i+1,j) - tmpFdx(i,j) )
                0174      &                     + ( tmpFdy(i,j+1) - tmpFdy(i,j) )
                0175              ENDDO
e64f43ca56 Jean*0176            ENDDO
                0177
e93c2e7dac Jean*0178 C--        Computational space Filter
8edc6fec00 Jean*0179            IF ( n.LE.nShapComput ) THEN
                0180              DO j=0,sNy+1
                0181               DO i=0,sNx+1
                0182                tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
                0183               ENDDO
                0184              ENDDO
e93c2e7dac Jean*0185 C--        Physical space Filter
8edc6fec00 Jean*0186            ELSEIF (Shap_TrLength.LE.0.) THEN
                0187              DO j=0,sNy+1
                0188               DO i=0,sNx+1
                0189                tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
                0190      &             *recip_hFacC(i,j,k,bi,bj)
                0191               ENDDO
e64f43ca56 Jean*0192              ENDDO
                0193            ELSE
8edc6fec00 Jean*0194              DO j=0,sNy+1
                0195               DO i=0,sNx+1
                0196                tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
                0197      &             *recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
                0198      &             *Shap_TrLength*Shap_TrLength
                0199               ENDDO
e64f43ca56 Jean*0200              ENDDO
                0201            ENDIF
8edc6fec00 Jean*0202 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
e64f43ca56 Jean*0203
8edc6fec00 Jean*0204 C      end k,bi,bj loop:
e64f43ca56 Jean*0205           ENDDO
                0206          ENDDO
                0207         ENDDO
8edc6fec00 Jean*0208 C      end loop n=1,nShapTr
aea29c8517 Alis*0209        ENDDO
                0210
15c70d7cd1 Jean*0211 C      F <-  [1 - (d_xx+d_yy)^n *deltaT/tau].F
aea29c8517 Alis*0212        DO bj=myByLo(myThid),myByHi(myThid)
                0213         DO bi=myBxLo(myThid),myBxHi(myThid)
8edc6fec00 Jean*0214          DO k=1,kSize
                0215           DO j=1,sNy
                0216            DO i=1,sNx
e64f43ca56 Jean*0217             field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
6824f48020 Jean*0218      &            -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
                0219             tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
aea29c8517 Alis*0220            ENDDO
                0221           ENDDO
                0222          ENDDO
                0223         ENDDO
                0224        ENDDO
                0225
e93c2e7dac Jean*0226        IF ( exchInOut.GE.2 ) THEN
                0227         IF (kSize.EQ.Nr) THEN
7163a40534 Jean*0228           _EXCH_XYZ_RL( field, myThid )
e93c2e7dac Jean*0229         ELSEIF (kSize.EQ.1) THEN
7163a40534 Jean*0230           _EXCH_XY_RL( field, myThid )
e93c2e7dac Jean*0231         ELSE
                0232           STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
                0233         ENDIF
                0234        ENDIF
aea29c8517 Alis*0235
                0236       ENDIF
                0237 #endif /* ALLOW_SHAP_FILT */
                0238
                0239       RETURN
                0240       END