! source file: /sfs/fs6/home-geomar/smomw258/UVic_ESCM/2.9/source/common/filtr.F
subroutine filtr (s, im, mm, n, iss)
!=======================================================================
! ===
! filter fourier analyses the arrays of various ===
! physical quantities, then truncates the series and ===
! resynthesizes the filtered quantities where: ===
! s =the string to be filtered ===
! im =the length of s ===
! mm =1 (cosine series, deriv at bndry pts=0) ===
! =2 ( sine series, bndry pts=0) ===
! =3 (full series, cyclic) ===
! n =number of waves to keep ===
! iss=0 (cant use fourier coefs from previous call) ===
! iss>0 (can use fourier coefs from previous call) ===
!=======================================================================
implicit none
!---------------------------------------------------------------------
! define global data
!---------------------------------------------------------------------
integer imtx2, ni, imtd2, lqmsum, lhsum, imtx4, imtx8, imtimt
integer imp1x2, im, mm, n, iss, imsave, i, imm1, imqc, nmax
integer nmaxp1, lcy, lh, lhm1, lqm, l2cy, lcym1, lcyp1, imx2
integer imx4, imx8, nprint, maxind, ncyc, maxndx, npwr, np, j
integer ioff1, ioff2, joff, ioff, jbase, ibase
real ssm, fimr, cc1, cc2, fnorm, ssum, fim, stemp, fact1
real fact2, genadj
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "scalar.h"
include "switch.h"
!---------------------------------------------------------------------
! define local data and dimension argument arrays
!---------------------------------------------------------------------
parameter (imtx2=imt*2,ni=imt)
parameter (imtd2=imt/2,lqmsum=imtd2*(imt-imtd2),lhsum=imt*imtp1/2)
parameter (imtx4=imt*4,imtx8=imt*8,imtimt=imt*imt)
parameter (imp1x2=imtp1*2)
integer icbase(imtp1), idbase(imtp1), ind(imtx8), indx(imtx8)
common /cfilt_i/ ind, idbase, icbase
common /cfilt_i/ imsave, jbase, ibase
real*8 cossav(lqmsum)
real temp(imtx4), denmsv(lhsum), cosnpi(imt)
real cof(imtx8), cosine(imtx8), ftarr(imtimt)
real denom(imtx4), s(imt), sprime(imt), circle(4)
common /cfilt_d/ cossav
common /cfilt_r/ denmsv, cosnpi, ftarr, circle
! data circle /0.,-1.,0.,1./
circle(1) = 0.0
circle(2) = -1.0
circle(3) = 0.0
circle(4) = 1.0
!---------------------------------------------------------------------
! begin executable code
!---------------------------------------------------------------------
if (im.lt.1 .or. mm.lt.1 .or. mm.gt.3 .or. n.lt.0 .or. iss.lt.0)
$ then
write (stdout,99) im, mm, n, iss
write (stderr,99) im, mm, n, iss
stop ' filtr 1'
endif
if (first) then
! this section sets up tables for filter; it must be called once
! per execution of ocean
! note: lqmsum is the sum of (im-1)/2 for im=1,imtp1
! lhsum is the sum of im-1 for im=1,imtp1
imsave = im
! assemble index array
do 100 i=1,imtx8
ind(i) = i
100 continue
! calculate and save all cosines which will be needed
ibase = 0
jbase = 0
do 200 im=1,imtp1
fimr = c1/float(im)
imm1 = im-1
if (imm1.eq.0) goto 181
do 180 i=1,imm1
denmsv(ibase+i) = c1/(c1-cos(pi*float(i)*fimr))
180 continue
181 continue
idbase(im) = ibase
ibase = ibase + imm1
imqc = (im-1)/2
if (imqc .eq. 0) goto 191
do 190 i=1,imqc
cossav(jbase+i) = cos(pi*float(i)*fimr)
190 continue
191 continue
icbase(im) = jbase
jbase = jbase + imqc
200 continue
! calculate adjustments for general fourier case if im=2*n
do 300 im=1,imt
cosnpi(im) = circle(mod(im-1,4)+1)
300 continue
im = imsave
endif
! calculate some useful constants
if (mm.eq.2 .and. n.eq.0) then
do 400 i=1,im
s(i) = c0
400 continue
goto 3201
endif
if (mm .eq. 1) then
nmax = n - 1
else
nmax = n
endif
nmaxp1 = nmax + 1
cc1 = p5*float(nmax) + p25
cc2 = float(nmax) + p5
if (mm .eq. 2) then
lcy = 2*(im + 1)
fnorm = c2/float(im + 1)
else
lcy = 2*im
fnorm = c2/float(im)
endif
lh = lcy/2
lhm1 = lh - 1
lqm = (lh - 1)/2
l2cy = 2*lcy
lcym1 = lcy - 1
lcyp1 = lcy + 1
imx2 = im*2
imx4 = im*4
imx8 = im*8
! average incoming array
ssum = c0
do 500 i=1,im
ssum = ssum + s(i)
500 continue
! mm = 1 derivative must be zero at boundaries (cosine)
! mm = 2 value must be zero at boundaries (sine)
! mm = 3 cyclic boundary conditions (general fourier series)
fim = float(im)
fimr = c1/fim
stemp = ssum*fimr
if (n.gt.1 .or. mm.ne.1) goto 601
do 600 i=1,im
s(i)=stemp
600 continue
go to 3201
601 continue
if (mm .ne. 2) then
do 700 i=1,im
s(i) = s(i) - stemp
700 continue
endif
if (iss .gt. 0) goto 2501
! assemble appropriate 1-cycle (2*pi) cosine array
! use stored 1/4 cycle to calculate first 1/2 cycle
jbase = icbase(lh)
do 800 i=1,lqm
cosine(i) = cossav(jbase+i)
800 continue
do 900 i=1,lqm
cosine(lh-i) = -cossav(jbase+i)
900 continue
! fill in cos(pi/2) if lh is even
if (2*(lqm+1) .eq. lh) cosine(lqm+1) = c0
! fill in cos(pi) in any case
cosine(lh) = -c1
! fill in rest of cycle
do 1000 i=1,lh
cosine(lh+i) = -cosine(i)
1000 continue
! assemble denominator array
ibase = idbase(lh)
do 1100 i=1,lhm1
denom(i) = p25*denmsv(ibase+i)
1100 continue
denom(lh) = 0.125
do 1200 i=1,lhm1
temp(i) = denom(lh-i)
1200 continue
do 1300 i=1,lhm1
denom(lh+i) = temp(i)
1300 continue
nprint = 0
denom(lcy) = c0
do 1400 i=lcyp1,imx4
denom(i) = denom(i-lcy)
1400 continue
! assemble appropriate subscript arrays
! calculate needed indices
if (mm.eq.3) then
fact1 = 2*nmax
fact2 = 2*nmaxp1
else
fact1 = nmax
fact2 = nmaxp1
endif
do 1500 i=1,imx4
indx(i) = ind(i)*fact1
1500 continue
do 1600 i=1,imx4
indx(imx4+i) = ind(i)*fact2
1600 continue
! calculate parameters for reducing indices
maxind = imx4*fact2
ncyc = (maxind-1)/lcy + 1
maxndx = lcy
if (maxndx .ge. maxind) goto 1801
do 1700 npwr=1,ncyc+2
maxndx = 2*maxndx
if (maxndx .ge. maxind) goto 1701
1700 continue
write (stdout,999)
write (stderr,999)
stop ' filtr 2'
1701 continue
do 1800 np=1,npwr
maxndx = maxndx/2
do 1790 i=1,imx8
if (indx(i) .gt. maxndx) indx(i) = indx(i) - maxndx
1790 continue
1800 continue
1801 continue
! gather coefficients
do 1900 j=1,imx8
cof(j) = cosine(indx(j))
1900 continue
! assemble transformation array which will filter s
if (mm.eq.1) then
! cosine transform
ioff1 = lcy
ioff2 = lcy + imx4
do 2000 j=1,im
joff = (j-1)*imt
do 1990 i=1,im
ftarr(joff+i) =
$ (cof(i-j+ioff1) - cof(i-j+ioff2)) *denom(i-j+ioff1) +
$ (cof(i + j - 1) - cof(imx4+i+j-1))*denom(i+j-1) - p5
1990 continue
2000 continue
do 2100 j=1,im
ftarr(j*imtp1-imt) = ftarr(j*imtp1-imt) + cc1
2100 continue
elseif (mm .eq. 2) then
! sine transform
ioff1 = lcy
ioff2 = lcy + imx4
do 2200 j=1,im
joff = (j-1)*imt
do 2190 i=1,im
ftarr(joff+i) =
$ (cof(i-j+ioff1) - cof(i-j+ioff2))*denom(i-j+ioff1) -
$ (cof(i + j) - cof(imx4+i+j)) *denom(i+j)
2190 continue
2200 continue
do 2300 j=1,im
ftarr(j*imtp1-imt) = ftarr(j*imtp1-imt) + cc1
2300 continue
elseif (mm.eq.3) then
! general fourier transform
if (2*n .eq. im) then
genadj = p5
else
genadj = c0
endif
ioff1 = lcy
ioff2 = lcy + imx4
do 2400 j=1,im
joff = (j-1)*imt
do 2390 i=1,im
ftarr(joff+i) = (c2*(cof(i-j+ioff1) - cof(i-j+ioff2)))
$ *denom(2*i-2*j+ioff1) - p5 - genadj*cosnpi(i)*cosnpi(j)
2390 continue
2400 continue
do 2500 j=1,im
ftarr(j*imtp1-imt) = ftarr(j*imtp1-imt) + cc2
2500 continue
endif
! filter s
2501 continue
do 2600 i=1,im
sprime(i) = c0
2600 continue
! note that ftarr(j,i)=ftarr(i,j), so following is legal
do 2700 i=1,im
ioff = (i-1)*imt
do 2690 j=1,im
sprime(j) = sprime(j) + s(i)*ftarr(ioff+j)
2690 continue
2700 continue
do 2800 i=1,im
sprime(i) = fnorm*sprime(i)
2800 continue
if (mm.eq.2) then
do 2900 i=1,im
s(i) = sprime(i)
2900 continue
goto 3201
endif
ssm = c0
do 3100 i=1,im
ssm = ssm + sprime(i)
3100 continue
ssm = (ssum-ssm)*fimr
do 3200 i=1,im
s(i) = ssm+sprime(i)
3200 continue
3201 continue
99 format (/' error => bad argument(s) in call to filtr'
$ /' im,mm,n,iss = ',4i10)
999 format (/' error => can not calculate parameters for reducing',
$ ' indices in filtr')
return
end