subroutine frfhvo(p1,iv) !$$$ subprogram documentation block ! . . . . ! subprogram: frfhvo performs vertical smoothing of fields ! prgmmr: derber org: np23 date: 2004-05-13 ! ! abstract: performs vertical smoothing of fields ! ! program history log: ! 2004-05-13 derber, document ! 2005-01-22 parrish - make use of balmod and rename variables ! 2005-07-14 wu - add max bound to l2 ! 2008-04-11 safford - rm unsed vars ! ! input argument list: ! p1 - input field to be smoothed (lat2,lon2,nsig) ! iv - location in alv for smoothing coefficients ! - iv = 1 streamfunction ! - iv = 2 velocity potential ! - iv = 3 temperature ! - iv = 4 specific humidity ! - iv = 5 ozone ! - iv = 6 cloud condensate ! ! output argument list ! p1 - output field after vertical smoothing ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind use constants, only: ione,zero,one use gridmod, only: nsig,regional,lat2,lon2 use balmod, only: rllat1,llmax use berror, only: alv,be,ndeg implicit none ! lat2 = number of latitudes (lat2) ! lon2 = number of longitudes (lon2) ! nsig = number of model levels ! ndeg = degree of smoothing (ndeg=4) integer(i_kind) ,intent(in ) :: iv real(r_kind),dimension(lat2,lon2,nsig),intent(inout) :: p1 integer(i_kind) kmod2,j,i,k,kr,ki,i1,l,l2 real(r_kind) dl1,dl2,alvl,alvr,alvi,alvi1,alvr1 real(r_kind) gaki,gakr,dekr,deki,bekr,beki real(r_kind),dimension(lat2,lon2,nsig):: p2 real(r_kind),dimension(lat2,lon2,ndeg):: ga,de real(r_kind) newr, newi, oldr, oldi real(r_kind) newr0, newi0, oldr0, oldi0 real(r_kind) newr1, newi1, oldr1, oldi1 kmod2=mod(ndeg,2_i_kind) ! Zero output array do k=1,nsig do j=1,lon2 do i=1,lat2 p2(i,j,k)=p1(i,j,k) p1(i,j,k)=zero end do end do end do ! Zero local work arrays do k=1,ndeg do j=1,lon2 do i=1,lat2 ga(i,j,k)=zero de(i,j,k)=zero end do end do end do if (regional)then ! Regional mode, odd degree if (kmod2 == ione) then ! advancing filter: do i=1,nsig do j=1,lon2 do k=1,lat2 l=int(rllat1(k,j)) l2=min0(l+ione,llmax) dl2=rllat1(k,j)-float(l) dl1=one-dl2 alvl=dl1*alv(l,1,i,iv)+dl2*alv(l2,1,i,iv) ga(k,j,1)=alvl*ga(k,j,1)+be(1)*p2(k,j,i) end do enddo do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+ga(k,j,1) end do enddo ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 l=int(rllat1(k,j)) l2=min0(l+ione,llmax) dl2=rllat1(k,j)-float(l) dl1=one-dl2 alvr=dl1*alv(l,kr,i,iv)+dl2*alv(l2,kr,i,iv) alvi=dl1*alv(l,ki,i,iv)+dl2*alv(l2,ki,i,iv) gakr=ga(k,j,kr) gaki=ga(k,j,ki) ga(k,j,kr)=alvr*gakr-alvi*gaki+bekr*p2(k,j,i) ga(k,j,ki)=alvi*gakr+alvr*gaki+beki*p2(k,j,i) end do enddo do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+ga(k,j,kr) end do enddo enddo enddo ! backing filter: do i=nsig,1,-1 do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+de(k,j,1) end do enddo do j=1,lon2 do k=1,lat2 l=int(rllat1(k,j)) l2=min0(l+ione,llmax) dl2=rllat1(k,j)-float(l) dl1=one-dl2 alvl=dl1*alv(l,1,i,iv)+dl2*alv(l2,1,i,iv) de(k,j,1)=alvl*(de(k,j,1)+be(1)*p2(k,j,i)) end do enddo ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+de(k,j,kr) end do enddo do j=1,lon2 do k=1,lat2 l=int(rllat1(k,j)) l2=min0(l+ione,llmax) dl2=rllat1(k,j)-float(l) dl1=one-dl2 alvr=dl1*alv(l,kr,i,iv)+dl2*alv(l2,kr,i,iv) alvi=dl1*alv(l,ki,i,iv)+dl2*alv(l2,ki,i,iv) dekr=de(k,j,kr)+bekr*p2(k,j,i) deki=de(k,j,ki)+beki*p2(k,j,i) de(k,j,kr)=alvr*dekr-alvi*deki de(k,j,ki)=alvi*dekr+alvr*deki end do enddo enddo enddo ! Regional mode, even degree else ! advancing filter: do i=1,nsig-ione,2 i1=i+ione ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 l=int(rllat1(k,j)) l2=min0(l+ione,llmax) dl2=rllat1(k,j)-float(l) dl1=one-dl2 alvr=dl1*alv(l,kr,i,iv)+dl2*alv(l2,kr,i,iv) alvi=dl1*alv(l,ki,i,iv)+dl2*alv(l2,ki,i,iv) alvr1=dl1*alv(l,kr,i1,iv)+dl2*alv(l2,kr,i1,iv) alvi1=dl1*alv(l,ki,i1,iv)+dl2*alv(l2,ki,i1,iv) gakr=ga(k,j,kr) gaki=ga(k,j,ki) ga(k,j,kr)=alvr*gakr-alvi*gaki+bekr*p2(k,j,i) ga(k,j,ki)=alvi*gakr+alvr*gaki+beki*p2(k,j,i) p1(k,j,i)=p1(k,j,i)+ga(k,j,kr) gakr=ga(k,j,kr) gaki=ga(k,j,ki) ga(k,j,kr)=alvr1*gakr-alvi1*gaki+bekr*p2(k,j,i1) ga(k,j,ki)=alvi1*gakr+alvr1*gaki+beki*p2(k,j,i1) p1(k,j,i1)=p1(k,j,i1)+ga(k,j,kr) end do enddo enddo enddo if (mod(nsig,2_i_kind)==ione) then i=nsig ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 l=int(rllat1(k,j)) l2=min0(l+ione,llmax) dl2=rllat1(k,j)-float(l) dl1=one-dl2 alvr=dl1*alv(l,kr,i,iv)+dl2*alv(l2,kr,i,iv) alvi=dl1*alv(l,ki,i,iv)+dl2*alv(l2,ki,i,iv) gakr=ga(k,j,kr) gaki=ga(k,j,ki) ga(k,j,kr)=alvr*gakr-alvi*gaki+bekr*p2(k,j,i) ga(k,j,ki)=alvi*gakr+alvr*gaki+beki*p2(k,j,i) p1(k,j,i)=p1(k,j,i)+ga(k,j,kr) end do enddo enddo endif ! backing filter: do i=nsig,2,-2 i1=i-ione ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 l=int(rllat1(k,j)) l2=min0(l+ione,llmax) dl2=rllat1(k,j)-float(l) dl1=one-dl2 alvr=dl1*alv(l,kr,i,iv)+dl2*alv(l2,kr,i,iv) alvi=dl1*alv(l,ki,i,iv)+dl2*alv(l2,ki,i,iv) alvr1=dl1*alv(l,kr,i1,iv)+dl2*alv(l2,kr,i1,iv) alvi1=dl1*alv(l,ki,i1,iv)+dl2*alv(l2,ki,i1,iv) p1(k,j,i)=p1(k,j,i)+de(k,j,kr) dekr=de(k,j,kr)+bekr*p2(k,j,i) deki=de(k,j,ki)+beki*p2(k,j,i) de(k,j,kr)=alvr*dekr-alvi*deki de(k,j,ki)=alvi*dekr+alvr*deki p1(k,j,i1)=p1(k,j,i1)+de(k,j,kr) dekr=de(k,j,kr)+bekr*p2(k,j,i1) deki=de(k,j,ki)+beki*p2(k,j,i1) de(k,j,kr)=alvr1*dekr-alvi1*deki de(k,j,ki)=alvi1*dekr+alvr1*deki end do enddo enddo enddo if (mod(nsig,2_i_kind)==ione) then i=ione ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 l=int(rllat1(k,j)) l2=min0(l+ione,llmax) dl2=rllat1(k,j)-float(l) dl1=one-dl2 alvr=dl1*alv(l,kr,i,iv)+dl2*alv(l2,kr,i,iv) alvi=dl1*alv(l,ki,i,iv)+dl2*alv(l2,ki,i,iv) p1(k,j,i)=p1(k,j,i)+de(k,j,kr) dekr=de(k,j,kr)+bekr*p2(k,j,i) deki=de(k,j,ki)+beki*p2(k,j,i) de(k,j,kr)=alvr*dekr-alvi*deki de(k,j,ki)=alvi*dekr+alvr*deki end do enddo enddo endif endif ! Global branch else ! Global mode, odd degree if (kmod2 == ione) then ! advancing filter: do i=1,nsig do j=1,lon2 do k=1,lat2 ga(k,j,1)=alv(k,1,i,iv)*ga(k,j,1)+be(1)*p2(k,j,i) end do enddo do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+ga(k,j,1) end do enddo ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 gakr=ga(k,j,kr) gaki=ga(k,j,ki) ga(k,j,kr)=alv(k,kr,i,iv)*gakr-alv(k,ki,i,iv)*gaki+bekr*p2(k,j,i) ga(k,j,ki)=alv(k,ki,i,iv)*gakr+alv(k,kr,i,iv)*gaki+beki*p2(k,j,i) end do enddo do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+ga(k,j,kr) end do enddo enddo enddo ! backing filter: do i=nsig,1,-1 do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+de(k,j,1) end do enddo do j=1,lon2 do k=1,lat2 de(k,j,1)=alv(k,1,i,iv)*(de(k,j,1)+be(1)*p2(k,j,i)) end do enddo ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+de(k,j,kr) end do enddo do j=1,lon2 do k=1,lat2 dekr=de(k,j,kr)+bekr*p2(k,j,i) deki=de(k,j,ki)+beki*p2(k,j,i) de(k,j,kr)=alv(k,kr,i,iv)*dekr-alv(k,ki,i,iv)*deki de(k,j,ki)=alv(k,ki,i,iv)*dekr+alv(k,kr,i,iv)*deki end do enddo enddo enddo ! Global branch, even degree else ! advancing filter: do i=1,nsig-ione,2 i1=i+ione ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2-ione, 2 oldr0=ga(k ,j,kr) oldr1=ga(k+ione,j,kr) oldi0=ga(k ,j,ki) oldi1=ga(k+ione,j,ki) newr0=alv(k ,kr,i,iv)*oldr0-alv(k ,ki,i,iv)*oldi0+bekr*p2(k ,j,i) newr1=alv(k+ione,kr,i,iv)*oldr1-alv(k+ione,ki,i,iv)*oldi1+bekr*p2(k+ione,j,i) newi0=alv(k ,ki,i,iv)*oldr0+alv(k ,kr,i,iv)*oldi0+beki*p2(k ,j,i) newi1=alv(k+ione,ki,i,iv)*oldr1+alv(k+ione,kr,i,iv)*oldi1+beki*p2(k+ione,j,i) p1(k ,j,i)=p1(k ,j,i)+newr0 p1(k+ione,j,i)=p1(k+ione,j,i)+newr1 ga(k ,j,kr)=alv(k ,kr,i1,iv)*newr0-alv(k ,ki,i1,iv)*newi0+bekr*p2(k ,j,i1) ga(k+ione,j,kr)=alv(k+ione,kr,i1,iv)*newr1-alv(k+ione,ki,i1,iv)*newi1+bekr*p2(k+ione,j,i1) ga(k ,j,ki)=alv(k ,ki,i1,iv)*newr0+alv(k ,kr,i1,iv)*newi0+beki*p2(k ,j,i1) ga(k+ione,j,ki)=alv(k+ione,ki,i1,iv)*newr1+alv(k+ione,kr,i1,iv)*newi1+beki*p2(k+ione,j,i1) p1(k ,j,i1)=p1(k ,j,i1)+ga(k ,j,kr) p1(k+ione,j,i1)=p1(k+ione,j,i1)+ga(k+ione,j,kr) end do do k=k,lat2 oldr=ga(k,j,kr) oldi=ga(k,j,ki) newr=alv(k,kr,i,iv)*oldr-alv(k,ki,i,iv)*oldi+bekr*p2(k,j,i) newi=alv(k,ki,i,iv)*oldr+alv(k,kr,i,iv)*oldi+beki*p2(k,j,i) p1(k,j,i)=p1(k,j,i)+newr ga(k,j,kr)=alv(k,kr,i1,iv)*newr-alv(k,ki,i1,iv)*newi+bekr*p2(k,j,i1) ga(k,j,ki)=alv(k,ki,i1,iv)*newr+alv(k,kr,i1,iv)*newi+beki*p2(k,j,i1) p1(k,j,i1)=p1(k,j,i1)+ga(k,j,kr) end do enddo enddo enddo if (mod(nsig,2_i_kind)==ione) then i=nsig ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 gakr=ga(k,j,kr) gaki=ga(k,j,ki) ga(k,j,kr)=alv(k,kr,i,iv)*gakr-alv(k,ki,i,iv)*gaki+bekr*p2(k,j,i) ga(k,j,ki)=alv(k,ki,i,iv)*gakr+alv(k,kr,i,iv)*gaki+beki*p2(k,j,i) p1(k,j,i)=p1(k,j,i)+ga(k,j,kr) end do enddo enddo endif ! backing filter: do i=nsig,2,-2 i1=i-ione ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2-ione, 2 p1(k ,j,i)=p1(k ,j,i)+de(k ,j,kr) p1(k+ione,j,i)=p1(k+ione,j,i)+de(k+ione,j,kr) oldr0=de(k ,j,kr)+bekr*p2(k ,j,i) oldr1=de(k+ione,j,kr)+bekr*p2(k+ione,j,i) oldi0=de(k ,j,ki)+beki*p2(k ,j,i) oldi1=de(k+ione,j,ki)+beki*p2(k+ione,j,i) newr0=alv(k ,kr,i,iv)*oldr0-alv(k ,ki,i,iv)*oldi0 newr1=alv(k+ione,kr,i,iv)*oldr1-alv(k+ione,ki,i,iv)*oldi1 newi0=alv(k ,ki,i,iv)*oldr0+alv(k ,kr,i,iv)*oldi0 newi1=alv(k+ione,ki,i,iv)*oldr1+alv(k+ione,kr,i,iv)*oldi1 p1(k ,j,i1)=p1(k ,j,i1)+newr0 p1(k+ione,j,i1)=p1(k+ione,j,i1)+newr1 newr0=newr0+bekr*p2(k ,j,i1) newr1=newr1+bekr*p2(k+ione,j,i1) newi0=newi0+beki*p2(k ,j,i1) newi1=newi1+beki*p2(k+ione,j,i1) de(k ,j,kr)=alv(k ,kr,i1,iv)*newr0-alv(k ,ki,i1,iv)*newi0 de(k+ione,j,kr)=alv(k+ione,kr,i1,iv)*newr1-alv(k+ione,ki,i1,iv)*newi1 de(k ,j,ki)=alv(k ,ki,i1,iv)*newr0+alv(k ,kr,i1,iv)*newi0 de(k+ione,j,ki)=alv(k+ione,ki,i1,iv)*newr1+alv(k+ione,kr,i1,iv)*newi1 end do do k=k,lat2 p1(k,j,i)=p1(k,j,i)+de(k,j,kr) oldr=de(k,j,kr)+bekr*p2(k,j,i) oldi=de(k,j,ki)+beki*p2(k,j,i) newr=alv(k,kr,i,iv)*oldr-alv(k,ki,i,iv)*oldi newi=alv(k,ki,i,iv)*oldr+alv(k,kr,i,iv)*oldi p1(k,j,i1)=p1(k,j,i1)+newr newr=newr+bekr*p2(k,j,i1) newi=newi+beki*p2(k,j,i1) de(k,j,kr)=alv(k,kr,i1,iv)*newr-alv(k,ki,i1,iv)*newi de(k,j,ki)=alv(k,ki,i1,iv)*newr+alv(k,kr,i1,iv)*newi end do enddo enddo enddo if (mod(nsig,2_i_kind)==ione) then i=ione ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components bekr=be(kr) beki=be(ki) do j=1,lon2 do k=1,lat2 p1(k,j,i)=p1(k,j,i)+de(k,j,kr) dekr=de(k,j,kr)+bekr*p2(k,j,i) deki=de(k,j,ki)+beki*p2(k,j,i) de(k,j,kr)=alv(k,kr,i,iv)*dekr-alv(k,ki,i,iv)*deki de(k,j,ki)=alv(k,ki,i,iv)*dekr+alv(k,kr,i,iv)*deki end do enddo enddo endif endif end if return end subroutine frfhvo subroutine smoothzo(vx,samp,rate,iv,jx,dsv) !$$$ subprogram documentation block ! . . . . ! subprogram: smoothzo initializes and renormalizes vertical smoothing coefs. ! prgmmr: derber org: np22 date: 2004-05-13 ! ! abstract: initializes and renormalizes vertical smoothing coefficients ! initializes dssv and alv ! ! program history log: ! 2004-05-13 derber, document ! 2004-11-30 treadon - add longitude dimension to variance array dssv ! 2010-03-01 zhu - decide the location in alv and dsv based on anavinfo ! - add dsv in the interface, rm dssv in berror ! ! input argument list: ! vx - vertical smoothing scales ! samp - parameter for smoothing ! rate - parameter for smoothing ! iv - location in alv and dssv for smoothing coefficients ! jx - latitude index ! ! output argument list ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind use constants, only: zero use gridmod, only: nsig,lon2 use berror, only: alv,ndeg implicit none integer(i_kind) ,intent(in ) :: jx,iv real(r_kind),dimension(nsig),intent(in ) :: vx real(r_kind) ,intent(in ) :: samp real(r_kind),dimension(ndeg),intent(in ) :: rate real(r_kind),dimension(lon2,nsig),intent(out):: dsv integer(i_kind) i,k,m real(r_kind),dimension(nsig):: dss real(r_kind),dimension(nsig,nsig):: p1 real(r_kind),dimension(nsig,ndeg):: al call rfdparv(vx,rate,al,nsig,ndeg) do m=1,ndeg do k=1,nsig alv(jx,m,k,iv)=al(k,m) end do end do p1=zero do k=1,nsig dss(k)=sqrt(samp*vx(k)) p1(k,k)=dss(k) end do call rfhvo(p1,nsig,nsig,al) call rfhvo(p1,nsig,nsig,al) do k=1,nsig do i=1,lon2 dsv(i,k)=sqrt(dss(k)/p1(k,k)) end do end do return end subroutine smoothzo subroutine rfhvo(p1,nc,n,al) !$$$ subprogram documentation block ! . . . . ! subprogram: rfhvo performs vertical smoothing for renormalization ! prgmmr: derber org: np22 date: 2004-05-13 ! ! abstract: performs vertical smoothing of identity matrix for use with ! renormalization ! ! program history log: ! 2004-05-13 derber, document ! ! input argument list: ! p1 - input field to be smoothed (nc,n) ! nc - first array dimension for p1 ! n - second array dimension for p1 ! al - smoothing coefficients ! ! output argument list ! p1 - output field after vertical smoothing ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind use constants, only: ione,zero use berror, only: be,ndeg implicit none integer(i_kind) ,intent(in ) :: n,nc real(r_kind),dimension(n,ndeg),intent(in ) :: al real(r_kind),dimension(nc,n) ,intent(inout) :: p1 integer(i_kind) i,j,kmod2,kr,ki real(r_kind) gaki,dekr,deki,gakr real(r_kind),dimension(nc,n):: p2 real(r_kind),dimension(nc,ndeg):: ga,de kmod2=mod(ndeg,2_i_kind) ! Zero local work arrays. do j=1,ndeg do i=1,nc ga(i,j)=zero de(i,j)=zero end do end do do j=1,n do i=1,nc p2(i,j)=p1(i,j) p1(i,j)=zero end do end do if (kmod2 == ione) then ! advancing filter: do i=1,n do j=1,nc ga(j,1)=al(i,1)*ga(j,1)+be(1)*p2(j,i) enddo do j=1,nc p1(j,i)=p1(j,i)+ga(j,1) enddo ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components do j=1,nc gakr=ga(j,kr) gaki=ga(j,ki) ga(j,kr)=al(i,kr)*gakr-al(i,ki)*gaki+be(kr)*p2(j,i) ga(j,ki)=al(i,ki)*gakr+al(i,kr)*gaki+be(ki)*p2(j,i) enddo do j=1,nc p1(j,i)=p1(j,i)+ga(j,kr) enddo enddo enddo ! backing filter: do i=n,1,-1 do j=1,nc p1(j,i)=p1(j,i)+de(j,1) enddo do j=1,nc de(j,1)=al(i,1)*(de(j,1)+be(1)*p2(j,i)) enddo ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components do j=1,nc p1(j,i)=p1(j,i)+de(j,kr) enddo do j=1,nc dekr=de(j,kr)+be(kr)*p2(j,i) deki=de(j,ki)+be(ki)*p2(j,i) de(j,kr)=al(i,kr)*dekr-al(i,ki)*deki de(j,ki)=al(i,ki)*dekr+al(i,kr)*deki enddo enddo enddo else ! advancing filter: do i=1,n ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components do j=1,nc gakr=ga(j,kr) gaki=ga(j,ki) ga(j,kr)=al(i,kr)*gakr-al(i,ki)*gaki+be(kr)*p2(j,i) ga(j,ki)=al(i,ki)*gakr+al(i,kr)*gaki+be(ki)*p2(j,i) enddo do j=1,nc p1(j,i)=p1(j,i)+ga(j,kr) enddo enddo enddo ! backing filter: do i=n,1,-1 ! treat remaining complex roots: do kr=kmod2+ione,ndeg,2 ! <-- index of "real" components ki=kr+ione ! <-- index of "imag" components do j=1,nc p1(j,i)=p1(j,i)+de(j,kr) enddo do j=1,nc dekr=de(j,kr)+be(kr)*p2(j,i) deki=de(j,ki)+be(ki)*p2(j,i) de(j,kr)=al(i,kr)*dekr-al(i,ki)*deki de(j,ki)=al(i,ki)*dekr+al(i,kr)*deki enddo enddo enddo endif return end subroutine rfhvo