! source file: /Users/oschlies/UVIC/master/source/mom/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, 1) === ! n =number of waves to keep === ! iss=0 (cant use fourier coefs from previous call) === ! iss>0 (can use fourier coefs from previous call) === ! based on code by: M. Cox === ! === !======================================================================= !--------------------------------------------------------------------- ! define global data !--------------------------------------------------------------------- include "param.h" include "ndcon.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) ! cossav must remain full precision if most of filter is made half-p real cossav dimension icbase(imtp1),idbase(imtp1),ind(imtx8),temp(imtx4) dimension cossav(lqmsum),denmsv(lhsum),cosnpi(imt) dimension circle(4) dimension indx(imtx8),cof(imtx8) dimension cosine(imtx8),ftarr(imtimt) dimension denom(imtx4) dimension s(imt),sprime(imt) common /cfilti/ ind, idbase, icbase common /cfilti/ imsave, jbase, ibase common /cfiltr/ denmsv, cossav, cosnpi, ftarr common /cfilt1/ 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 1 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