/[MITgcm]/MITgcm/pkg/zonal_filt/zonal_filter.F
ViewVC logotype

Diff of /MITgcm/pkg/zonal_filt/zonal_filter.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by adcroft, Wed Jan 24 16:56:07 2001 UTC revision 1.2 by adcroft, Fri Feb 2 21:36:30 2001 UTC
# Line 0  Line 1 
1    C $Header$
2    
3    #include "ZONAL_FILT_OPTIONS.h"
4    
5          SUBROUTINE ZONAL_FILTER(
6         U           field, fieldMask,
7         I           jMin, jMax, kMin, kMax, bi, bj, gridLoc, myThid )
8    C     /==========================================================\
9    C     | S/R ZONAL_FILTER                                         |
10    C     | o Apply FFT filter to a latitude circle.                 |
11    C     \==========================================================/
12          IMPLICIT NONE
13    
14    C     == Global data ==
15    #include "SIZE.h"
16    #include "ZONAL_FILT.h"
17    #include "FFTPACK.h"
18    
19    C     == Routine arguments ==
20    C     jMin - Range of points to filter
21    C     jMax
22    C     kMin
23    C     kMax
24    C     bi
25    C     bj
26    C     myThid  - Thread number of this instance of FILTER_LATCIRC_FFT_APPLY
27    C     field   - Field to filter
28    C     gridLoc - Orientation (U or V) of field.
29          INTEGER myThid
30          INTEGER gridLoc
31          Real*8 field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
32          Real*8 fieldMask(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
33          INTEGER jMin, jMax, kMin, kMax, bi, bj
34    
35    #ifdef ALLOW_ZONAL_FILT
36    
37    C     == Local data ==
38          Real*8 phi(Nx)
39          Real*8 phiMask(Nx)
40          Real*8 avPhi
41          INTEGER I, J, K
42    
43          DO k=kMin, kMax
44           DO j=jMin, jMax
45    
46    C     o Copy zonal line of field into local workspace
47            DO i=1,sNx
48             phi(I) = field(i,j,k,bi,bj)
49             phiMask(I) = fieldMask(i,j,k,bi,bj)
50            ENDDO
51    
52    C Interpolate through land
53            CALL ZONAL_FILT_PRESMOOTH( phiMask,phi,avPhi,sNx,myThid )
54    
55    C     o Forward transform (using specific FFT package)
56    C       CALL R8FFTF( Nx, phi, FFTPACKWS(1,bj) )
57            CALL R8FFTF1( Nx, phi,
58         &    FFTPACKWS1(1,bj), FFTPACKWS2(1,bj),FFTPACKWS3(1,bj) )
59    
60    C     o Apply amplitude filter and normalize
61            IF (gridLoc .EQ. 1) THEN
62             DO i=1, Nx
63              phi(i)=phi(i)*ampFactor(i,j,bi,bj)/float(Nx)
64             ENDDO
65            ELSEIF (gridLoc .EQ. 2) THEN
66             DO i=1, Nx
67              phi(i)=phi(i)*ampFactorV(i,j,bi,bj)/float(Nx)
68             ENDDO
69            ELSE
70             WRITE(0,*) 'Error: gridLoc = ',gridLoc
71             STOP 'Error: gridLoc has illegal value'
72            ENDIF
73    
74    C     o Backward transform (using specific FFT package)
75    C       CALL R8FFTB( Nx, phi, FFTPACKWS(1,bj) )
76            CALL R8FFTB1( Nx, phi,
77         &    FFTPACKWS1(1,bj), FFTPACKWS2(1,bj),FFTPACKWS3(1,bj) )
78    
79    C De-interpolate through land
80            CALL ZONAL_FILT_POSTSMOOTH(phiMask,phi,avPhi,sNx,myThid)
81    
82    C       o Do periodic wrap around by hand
83            DO i=1-OLx,0
84             field(i,j,k,bi,bj) = phi(sNx+i)
85            ENDDO
86            DO i=1,sNx
87             field(i,j,k,bi,bj) = phi(I)
88            ENDDO
89            DO i=sNx+1,sNx+OLx
90             field(i,j,k,bi,bj) = phi(i-sNx)
91            ENDDO
92    
93           ENDDO
94          ENDDO
95    
96    #endif /* ALLOW_ZONAL_FILT */
97    
98          RETURN
99          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22