m_shapiro.F90 1.91 KB
``````module m_shapiro
contains
subroutine shapiro(order,y,idim)
integer, intent(in)    :: order
integer, intent(in)    :: idim
real,    intent(inout) :: y(idim)

integer, parameter :: shdim=16
real sh(0:shdim)

call shapiro_ini(order,sh,shdim)
call shfilt(sh,shdim,order,y,idim)
end subroutine

subroutine shapiro_ini(n,sh,shdim)
integer, intent(in) :: shdim
real,    intent(out):: sh(0:shdim)
integer, intent(in) :: n

integer i,j
real f2n,fj,f2nj,ff

if((n == 1).OR.(n == 2).OR.(n == 4).OR.(n == 8))then
ff=2.0**(2*n)
f2n=1

do j=1,2*n
f2n=f2n*float(j)
enddo

fj=1
sh(0)=((-1.0)**(n-1))/ff

do j=1,n
f2nj=1
do i=1,2*n-j
f2nj=f2nj*float(i)
enddo
fj=fj*float(j)
sh(j)=(-1.0)**(n+1-j)*f2n/(ff*fj*f2nj)
enddo
else
write(*,*)'error in shfact.  n=',n
stop
endif
end subroutine shapiro_ini

subroutine shfilt(sh,shdim,ish,y,yn)
integer, intent(in) :: shdim
real,    intent(in) :: sh(0:shdim)
integer, intent(in) :: ish
integer, intent(in) :: yn
real,    intent(out):: y(yn)
integer i,j,m

real, allocatable :: z(:)

allocate(z(yn))

z=y

do m=2,ish
do j=0,ish-1
mm=1-(m-ish+j)
if(mm.GE.0)then
z(m)=z(m)+sh(j)*(2.0*y(1)-y(1+mm)+y(m+ish-j))
else
z(m)=z(m)+sh(j)*(y(m+ish-j)+y(m-ish+j))
endif
enddo
z(m)=(sh(ish))*y(m)+z(m)
enddo

do m=ish+1,yn-ish
do j=0,ish-1
z(m)=z(m)+sh(j)*(y(m+ish-j)+y(m-ish+j))
enddo
z(m)=(sh(ish))*y(m)+z(m)
enddo

do m=yn-ish+1,yn-1
do j=0,ish-1
mm=(m+ish-j)-yn
if(mm.GE.0)then
z(m)=z(m)+sh(j)*(2.0*y(yn)-y(yn-mm)+y(m-ish+j))
else
z(m)=z(m)+sh(j)*(y(m+ish-j)+y(m-ish+j))
endif
enddo
z(m)=(sh(ish))*y(m)+z(m)
enddo

y=z

deallocate(z)

end subroutine shfilt

end module
``````