subroutine da_test_spectral (max_wave, sizec, xbx, field) 3,11
!-----------------------------------------------------------------------
! Purpose: TBD
!-----------------------------------------------------------------------
implicit none
integer, intent(in) :: max_wave ! Maximum wavenumber.
integer, intent(in) :: sizec ! Size of complex cv.
type (xbx_type), intent(in) :: xbx ! For header & non-grid arrays.
real, intent(in) :: field(its:ite,jts:jte) ! Gridpoint field.
real :: field_out(its:ite,jts:jte)
real*8 :: power(0:max_wave) ! Power spectrum
real :: rcv(1:2*sizec) ! Spectral modes.
real :: rcv_out(1:2*sizec) ! Spectral modes.
integer :: m,mm, index_start, index_end
complex :: r_leg(1:jde)
complex :: ccv(1:sizec) ! Spectral modes.
complex :: ccv1(1:sizec) ! Spectral modes.
real :: den, num, xx
if (trace_use) call da_trace_entry
("da_test_spectral")
write (unit=stdout, fmt='(A)') &
' Test orthogonality of Associated Legendre Polynomials:'
! Initialise Power spectrum
power = 1.0
call da_setlegpol_test( jde, max_wave, xbx%alp_size, xbx%int_wgts, xbx%alp )
write(unit=stdout,fmt='(A)') &
' Test invertibility of spectral (Fourier, Legendre) transforms:'
! Gridpoint to spectral:
rcv = 0.0
call da_vv_to_v_spectral( xbx%ni, xbx%nj, max_wave, xbx%inc, xbx%lenr, xbx%lensav, xbx%lenwrk, &
xbx%alp_size, 2*sizec, xbx%alp, xbx%wsave, xbx%int_wgts, rcv, field)
field_out = 0.0
! Spectral to gridpoint:
call da_vtovv_spectral( max_wave, sizec, &
xbx % lenr, xbx % lenwrk, xbx % lensav, &
xbx % inc, xbx % alp_size, xbx % alp, &
xbx % wsave, power, rcv, field_out)
write(unit=stdout,fmt='(1x,a,e30.10)') &
'Domain-Averaged (Grid->Spectral->Grid) Error = ', &
sqrt( sum( ( field_out(1:xbx%ni,1:xbx%nj) - field(1:xbx%ni,1:xbx%nj) )**2 ) / &
sum( field(1:xbx%ni,1:xbx%nj)**2 ) )
rcv_out = 0.0
! Gridpoint to spectral (again):
call da_vv_to_v_spectral( xbx%ni, xbx%nj, max_wave, xbx%inc, xbx%lenr, xbx%lensav, xbx%lenwrk, &
xbx%alp_size, 2*sizec, xbx%alp, xbx%wsave, xbx%int_wgts, rcv_out, field_out)
rcv_out(1:2*sizec) = rcv_out(1:2*sizec) - rcv(1:2*sizec) ! Create difference for test diags.
write(unit=stdout,fmt='(1x,a,e30.10)') &
' Domain-Averaged (Spectral->Grid->Spectral) Error = ', &
sqrt( sum( rcv_out(1:2*sizec)**2 ) ) / sqrt( sum( rcv(1:2*sizec)**2 ) )
! Adjoint test for Spectral Transform
rcv_out = 0.0
call da_vtovv_spectral_adj( max_wave, sizec, &
xbx % lenr, xbx % lenwrk, xbx % lensav, &
xbx % inc, xbx % alp_size, xbx % alp, &
xbx % wsave, power, rcv_out, field_out)
write(unit=stdout,fmt='(A)') ' Adjoint test result for Spectral -> Grid : '
write(unit=stdout,fmt='(1x,a,e30.10)') &
' LHS ( LX.LX) = ',&
sum( field_out(1:xbx%ni,1:xbx%nj)*field_out(1:xbx%ni,1:xbx%nj) )
write(unit=stdout,fmt='(1x,a,e30.10)') &
' RHS ( X.L^TLX ) = ', sum( rcv(1:2*sizec)*rcv_out(1:2*sizec) )
! Inverse test for Legendre Transform
write(unit=stdout,fmt='(A)') ' Inverse and Adjoint Legendre test result:'
do m = 0, max_wave
index_start = m * ( max_wave + 1 - m ) + m * ( m + 1 ) / 2 + 1
index_end = index_start + max_wave - m
do mm = index_start, index_end
if (2*mm > 2*sizec) then
call da_error(__FILE__,__LINE__, &
(/"rcv_out index bounce"/))
end if
ccv(mm) = cmplx (rcv_out(2*mm-1), rcv_out(2*mm))
end do
r_leg = 0.0
call da_legtra_inv( jds, jde, jts, jte, max_wave, xbx%alp_size, m, xbx%alp, &
ccv(index_start:index_end), r_leg )
ccv1(index_start:index_end) = 0.0
call da_legtra ( xbx%nj, max_wave, xbx%alp_size, m, xbx%int_wgts, xbx%alp, r_leg, &
ccv1(index_start:index_end) )
ccv1(index_start:index_end) = ccv1(index_start:index_end) - &
ccv(index_start:index_end)
num = sum ( real(ccv1(index_start:index_end))*real(ccv1(index_start:index_end))+&
aimag(ccv1(index_start:index_end))* aimag(ccv1(index_start:index_end)) )
den = sum ( real(ccv(index_start:index_end))*real(ccv(index_start:index_end))+&
aimag(ccv(index_start:index_end))* aimag(ccv(index_start:index_end)) )
write(unit=stdout,fmt='(A,I4,A,E30.20)') &
'For zonal wave number',m,' difference ',sqrt(num)/sqrt(den)
xx = sum( real(r_leg(1:xbx%nj))* real(r_leg(1:xbx%nj))+ &
aimag(r_leg(1:xbx%nj))*aimag(r_leg(1:xbx%nj)) )
write(unit=stdout,fmt='(a,i5,a,e30.20)') 'For Wave = ',m,' LX.LX = ',xx
ccv1(index_start:index_end) = 0.0
call da_legtra_inv_adj( jds, jde, jts, jte, max_wave, xbx%alp_size, m, xbx%alp, &
ccv1(index_start:index_end), r_leg )
xx = sum( real(ccv(index_start:index_end))* &
real(ccv1(index_start:index_end)) +&
aimag(ccv(index_start:index_end))* &
aimag(ccv1(index_start:index_end)) )
write(unit=stdout,fmt='(a,i5,a,e30.20,/)') 'For Wave = ',m,' X.L^TLX = ',xx
end do
if (trace_use) call da_trace_exit("da_test_spectral")
end subroutine da_test_spectral