27.1.1 fourfil

Due to the convergence of meridians on a sphere, longitudinal grid resolution decreases toward zero as the poles are approached. For global domains this may severely limit the length of the time step due to the CFL constraint given in Section 14.4.4. The instability can be removed by filtering (Bryan, Manabe, and Pacanowski, 1975 and Takacs, Balgovind, 1983) the smallest scales out of the solution since they impose the most severe restriction on the time step.

Filtering is not to be encouraged and is not recommended unless absolutely necessary. When necessary, it should be restricted to the northernmost and southernmost polar latitudes in the domain. Typically, the filtering works on strips of latitude defined by the researcher. In the southern hemisphere, the reference latitude is given by filter_reflat_s and in the northern hemisphere, the reference latitude is given by filter_reflat_n. By default, these values are set to $70^\circ$ N and $70^\circ$ S latitude for the test case resolution. With higher resolution, these reference latitudes can and should be pushed further poleward. To save computation over the land mass of Antarctica, filtering is turned off southward of variable rjfrst which is defaulted to $81^\circ$ S. If needed, these variables can be set through namelist. Refer to Section 14.4.6.

Option fourfil does a Fourier smoothing of prognostic variables in the longitudinal direction by truncating wavenumbers larger than a critical one given by

$$N= im\frac{\cstj}{\cos{filter\_reflat\_s}}$$ (27.1)

where N is the number of waves to retain, im is the length of a strip of ocean data along the latitude given by $\cstj$, and $filter\_reflat\_s$ is a reference latitude row in the southern hemisphere ( $filter\_reflat\_n$ would be used for the northern hemisphere.). If data were being filtered on the latitude of U-cells, then $\cstj$ would be replaced by $\csuj$. Because of geometry and topography break up each latitude circle, each filtering latitude is composed of strips defined by starting and stopping longitudes for ocean data and each strip is therefore a function of latitude and depth. There are a lot of different sized strips which rules out the use of an FFT. Also, FFT's are inherently square whereas for filtering purposes, fewer than im wavenumbers are generally needed. Instead, a fourier transform is used which has been optimized in prehistoric times. As such, it is difficult to understand^27.1 exactly how this fourier decomposition and systhesis works by looking at the code since no reference was left by the original designer (Mike Cox who is now deceased). In principle the filter is doing the following.

Any variable $\alpha_i$ where i=1,im can be decomposed and re-constructed in terms of a discrete Fourier transform given by

$$\tilde{\alpha}_i = A_\circ + \sum_{\ell=1}^{N/2}w_\ell A_\ell\cos... ...2\pi i \ell}{im} + \sum_{\ell=1}^{N/2-1}w_\ell B_\ell\sin\frac{2\pi i \ell}{im}$$ (27.2)

where $\tilde{\alpha}_i$ is the filtered value of $\alpha$composed of N wavenumbers ($N \le im$) given by Equation (27.1) and the even and odd Fourier coefficients $A_\ell$and $B_\ell$ are

$$A_\circ \frac{1}{im} \sum_{i=1}^{im} \alpha_i$$ = (27.3)

$$A_\ell \frac{2}{im} \sum_{i=1}^{im} \alpha_i \cos\frac{2\pi i \ell}{im} \;\;\;for \;\ell=1,2,\cdots \leq \frac{N}{2}-1$$ = (27.4)

$$\frac{1}{im} \sum_{i=1}^{im} \alpha_i \cos(i \pi)$$ A_im/2 = (27.5)

$$B_\ell \frac{2}{im} \sum_{i=1}^{im} \alpha_i \sin\frac{2\pi i \ell}{im} \;\;\;for \;\ell=1,2,\cdots \leq \frac{N}{2}-1$$ = (27.6)

The window $w_\ell$ is given by the boxcar function

$$w_\ell = 1$$ (27.7)

which can lead to oscillations when wavenumbers greater than N are truncated. These oscillations are known as the Gibbs effect which is particularly pronounced if energy exists at the truncation wavenumber N. The Gibbs effect can increase variance and this is a drawback. The Gibbs effect can be minimized by multiplying the wavenumbers by a window. A typical cosine window is

$$w_\ell = \cos((\pi/2)\ell/N).$$ (27.8)

Tracers are filtered using even Fourier coefficients (cosines) because the boundary condition at edges of the tracer strips is no-flux. Odd Fourier coefficients (sines) are used for velocity components because of the no-slip boundary condition at the edges of the velocity strips. Vorticity ( ztd_i,jrow) is filtered with even Fourier coefficients^27.2. Strips which completely encircle the earth are filtered using the full series of sines and cosines using a cyclic boundary condition. Any vector which is filtered must first remove the wavenumber which is due to the changing direction of a latitude circle. This is wavenumber one for all latitudes. For instance, the x-component of constant cross polar flow will exhibit a wavenumber one influence. After the filtering operation, this wavenumber one due to the changing direction of the latitude circle must be put back.

Apart from the Gibbs problem, another drawback of fourier filtering is that two strips of different length on the same latitude circle will contain different wavenumbers after filtering. The reason is that truncation wavenumber is a function of strip length and not total distance around the latitude circle. If the two strips were on adjacent latitudes, then taking the divergence of horizontal velocity might introduce ficticious vertical velocities (because of differing scales on the strips. Note that option firfil does not have this problem.). Also, fourier filtering can be a real time burner and FFT's are of little help because of the arbitrary strip widths and the need to enforce proper boundary conditions at the side walls (i.e. zero filling of land will not work). It should be noted that filtering a function with a fourier transform and the cosine window given above yields results very similar to the much faster finite impulse filter described below. Also, the finite impulse filter does not induce a ficticious vertical velocity because the filtering is not a function of strip size.