pkg/zonal_filt/zonal_filter.F

42c525bfb4 Alis*0001 #include "ZONAL_FILT_OPTIONS.h"
                0002
e4b65e705c Jean*0003 CBOP 0
                0004 C     !ROUTINE: ZONAL_FILTER
                0005
                0006 C     !INTERFACE:
bc49b7aac9 Jean*0007       SUBROUTINE ZONAL_FILTER(
                0008      U           field,
                0009      I           fieldMask,
e4b65e705c Jean*0010      I           jMin, jMax, kSize, bi, bj, gridLoc, myThid )
bc49b7aac9 Jean*0011
e4b65e705c Jean*0012 C     !DESCRIPTION:
bc49b7aac9 Jean*0013 C     *==========================================================*
                0014 C     | S/R ZONAL_FILTER
                0015 C     | o Apply FFT filter to a latitude circle.
                0016 C     *==========================================================*
                0017
e4b65e705c Jean*0018 C     !USES:
42c525bfb4 Alis*0019       IMPLICIT NONE
                0020 C     == Global data ==
                0021 #include "SIZE.h"
d4c5f8fe62 Jean*0022 #include "EEPARAMS.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
42c525bfb4 Alis*0025 #include "ZONAL_FILT.h"
                0026 #include "FFTPACK.h"
                0027
e4b65e705c Jean*0028 C     !INPUT/OUTPUT PARAMETERS:
                0029 C     jMin      :: Range of points to filter
42c525bfb4 Alis*0030 C     jMax
e4b65e705c Jean*0031 C     kSize     :: Number of levels to filter
                0032 C     bi, bj    :: tile indices
                0033 C     field     :: Field to filter
                0034 C     fieldMask :: mask corresponding to field to filter
                0035 C     gridLoc   :: Position on the grid (U or V) of field.
                0036 C     myThid    :: my Thread Id number
                0037       INTEGER kSize
                0038       _RL     field    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize)
                0039       _RS     fieldMask(1-Olx:sNx+Olx,1-Oly:sNy+Oly,kSize)
                0040       INTEGER jMin, jMax, bi, bj
bc49b7aac9 Jean*0041       INTEGER gridLoc
                0042       INTEGER myThid
e4b65e705c Jean*0043 CEOP
42c525bfb4 Alis*0044
                0045 #ifdef ALLOW_ZONAL_FILT
                0046
e4b65e705c Jean*0047 C     !LOCAL VARIABLES:
42c525bfb4 Alis*0048       Real*8 phi(Nx)
                0049       Real*8 phiMask(Nx)
                0050       Real*8 avPhi
e4b65e705c Jean*0051       INTEGER i, j, k
42c525bfb4 Alis*0052
e4b65e705c Jean*0053       DO k= 1, kSize
42c525bfb4 Alis*0054        DO j=jMin, jMax
d4c5f8fe62 Jean*0055         IF ( (gridLoc.EQ.1 .AND.ABS(yC(1,j,bi,bj)).GE.zonal_filt_lat)
                0056      &   .OR.(gridLoc.EQ.2 .AND.ABS(yG(2,j,bi,bj)).GE.zonal_filt_lat)
                0057      &   .OR. zonal_filt_mode2dx.EQ.2 ) THEN
42c525bfb4 Alis*0058
                0059 C     o Copy zonal line of field into local workspace
                0060         DO i=1,sNx
e4b65e705c Jean*0061          phi(i) = field(i,j,k)
                0062          phiMask(i) = fieldMask(i,j,k)
42c525bfb4 Alis*0063         ENDDO
                0064
                0065 C Interpolate through land
                0066         CALL ZONAL_FILT_PRESMOOTH( phiMask,phi,avPhi,sNx,myThid )
                0067
                0068 C     o Forward transform (using specific FFT package)
                0069 C       CALL R8FFTF( Nx, phi, FFTPACKWS(1,bj) )
                0070         CALL R8FFTF1( Nx, phi,
                0071      &    FFTPACKWS1(1,bj), FFTPACKWS2(1,bj),FFTPACKWS3(1,bj) )
                0072
                0073 C     o Apply amplitude filter and normalize
                0074         IF (gridLoc .EQ. 1) THEN
                0075          DO i=1, Nx
                0076           phi(i)=phi(i)*ampFactor(i,j,bi,bj)/float(Nx)
                0077          ENDDO
                0078         ELSEIF (gridLoc .EQ. 2) THEN
                0079          DO i=1, Nx
                0080           phi(i)=phi(i)*ampFactorV(i,j,bi,bj)/float(Nx)
                0081          ENDDO
                0082         ELSE
2cfc9d59a2 Patr*0083          WRITE(*,*) 'Error: gridLoc = ',gridLoc
42c525bfb4 Alis*0084          STOP 'Error: gridLoc has illegal value'
                0085         ENDIF
                0086
                0087 C     o Backward transform (using specific FFT package)
                0088 C       CALL R8FFTB( Nx, phi, FFTPACKWS(1,bj) )
                0089         CALL R8FFTB1( Nx, phi,
                0090      &    FFTPACKWS1(1,bj), FFTPACKWS2(1,bj),FFTPACKWS3(1,bj) )
                0091
                0092 C De-interpolate through land
                0093         CALL ZONAL_FILT_POSTSMOOTH(phiMask,phi,avPhi,sNx,myThid)
                0094
                0095 C       o Do periodic wrap around by hand
                0096         DO i=1-OLx,0
e4b65e705c Jean*0097          field(i,j,k) = phi(sNx+i)
42c525bfb4 Alis*0098         ENDDO
                0099         DO i=1,sNx
e4b65e705c Jean*0100          field(i,j,k) = phi(i)
42c525bfb4 Alis*0101         ENDDO
                0102         DO i=sNx+1,sNx+OLx
e4b65e705c Jean*0103          field(i,j,k) = phi(i-sNx)
42c525bfb4 Alis*0104         ENDDO
                0105
d4c5f8fe62 Jean*0106         ENDIF
42c525bfb4 Alis*0107        ENDDO
                0108       ENDDO
                0109
                0110 #endif /* ALLOW_ZONAL_FILT */
                0111
                0112       RETURN
                0113       END