da_test_spectral.inc
References to this file elsewhere.
1 subroutine da_test_spectral( max_wave, sizec, xbx, field)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: TBD
5 !-----------------------------------------------------------------------
6
7 implicit none
8
9 integer, intent(in) :: max_wave ! Maximum wavenumber.
10 integer, intent(in) :: sizec ! Size of complex cv.
11 type (xbx_type),intent(in) :: xbx ! For header & non-grid arrays.
12
13 real, dimension(its:ite,jts:jte), intent(in) :: field ! Gridpoint field.
14
15 real, dimension(its:ite,jts:jte) :: field_out
16 real, dimension(0:max_wave) :: power ! Power spectrum
17 real, dimension(1:2*sizec) :: rcv ! Spectral modes.
18 real, dimension(1:2*sizec) :: rcv_out ! Spectral modes.
19
20 integer :: m,mm, index_start, index_end
21
22 complex :: r_leg(1:jde)
23 complex, dimension(1:sizec) :: ccv, ccv1 ! Spectral modes.
24 real :: den, num, xx
25
26 write (unit=stdout, fmt='(A)') &
27 ' Test orthogonality of Associated Legendre Polynomials:'
28
29 ! Initialise Power spectrum
30 power = 1.0
31 call da_setlegpol_test( jde, max_wave, xbx%alp_size, xbx%int_wgts, xbx%alp )
32
33 write(unit=stdout,fmt='(A)') &
34 ' Test invertibility of spectral (Fourier, Legendre) transforms:'
35
36 ! Gridpoint to spectral:
37 rcv = 0.
38 call da_vv_to_v_spectral( xbx%ni, xbx%nj, max_wave, xbx%inc, xbx%lenr, xbx%lensav, xbx%lenwrk, &
39 xbx%alp_size, 2*sizec, xbx%alp, xbx%wsave, xbx%int_wgts, rcv, field)
40
41 field_out = 0.
42 ! Spectral to gridpoint:
43 call da_vtovv_spectral( max_wave, sizec, &
44 xbx % lenr, xbx % lenwrk, xbx % lensav, &
45 xbx % inc, xbx % alp_size, xbx % alp, &
46 xbx % wsave, power, rcv, field_out)
47
48 write(unit=stdout,fmt='(1x,a,e30.10)') &
49 'Domain-Averaged (Grid->Spectral->Grid) Error = ', &
50 sqrt( sum( ( field_out(1:xbx%ni,1:xbx%nj) - field(1:xbx%ni,1:xbx%nj) )**2 ) / &
51 sum( field(1:xbx%ni,1:xbx%nj)**2 ) )
52 rcv_out = 0.
53
54 ! Gridpoint to spectral (again):
55 call da_vv_to_v_spectral( xbx%ni, xbx%nj, max_wave, xbx%inc, xbx%lenr, xbx%lensav, xbx%lenwrk, &
56 xbx%alp_size, 2*sizec, xbx%alp, xbx%wsave, xbx%int_wgts, rcv_out, field_out)
57
58 rcv_out(1:2*sizec) = rcv_out(1:2*sizec) - rcv(1:2*sizec) ! Create difference for test diags.
59
60 write(unit=stdout,fmt='(1x,a,e30.10)') &
61 ' Domain-Averaged (Spectral->Grid->Spectral) Error = ', &
62 sqrt( sum( rcv_out(1:2*sizec)**2 ) ) / sqrt( sum( rcv(1:2*sizec)**2 ) )
63
64 ! Adjoint test for Spectral Transform
65 rcv_out = 0.0
66 call da_vtovv_spectral_adj( max_wave, sizec, &
67 xbx % lenr, xbx % lenwrk, xbx % lensav, &
68 xbx % inc, xbx % alp_size, xbx % alp, &
69 xbx % wsave, power, rcv_out, field_out)
70
71 write(unit=stdout,fmt='(A)') ' Adjoint test result for Spectral -> Grid : '
72 write(unit=stdout,fmt='(1x,a,e30.10)') &
73 ' LHS ( LX.LX) = ',&
74 sum( field_out(1:xbx%ni,1:xbx%nj)*field_out(1:xbx%ni,1:xbx%nj) )
75 write(unit=stdout,fmt='(1x,a,e30.10)') &
76 ' RHS ( X.L^TLX ) = ', sum( rcv(1:2*sizec)*rcv_out(1:2*sizec) )
77
78 ! Inverse test for Legendre Transform
79
80 write(unit=stdout,fmt='(A)') ' Inverse and Adjoint Legendre test result:'
81
82 do m = 0, max_wave
83 index_start = m * ( max_wave + 1 - m ) + m * ( m + 1 ) / 2 + 1
84 index_end = index_start + max_wave - m
85
86 do mm = index_start, index_end
87 if (2*mm > 2*sizec) then
88 call da_error(__FILE__,__LINE__, &
89 (/"rcv_out index bounce"/))
90 end if
91 ccv(mm) = cmplx (rcv_out(2*mm-1), rcv_out(2*mm))
92 end do
93 r_leg = 0.0
94 call da_legtra_inv( jds, jde, jts, jte, max_wave, xbx%alp_size, m, xbx%alp, &
95 ccv(index_start:index_end), r_leg )
96
97 ccv1(index_start:index_end) = 0.
98 call da_legtra ( xbx%nj, max_wave, xbx%alp_size, m, xbx%int_wgts, xbx%alp, r_leg, &
99 ccv1(index_start:index_end) )
100 ccv1(index_start:index_end) = ccv1(index_start:index_end) - &
101 ccv(index_start:index_end)
102 num = sum ( real(ccv1(index_start:index_end))*real(ccv1(index_start:index_end))+&
103 aimag(ccv1(index_start:index_end))* aimag(ccv1(index_start:index_end)) )
104 den = sum ( real(ccv(index_start:index_end))*real(ccv(index_start:index_end))+&
105 aimag(ccv(index_start:index_end))* aimag(ccv(index_start:index_end)) )
106 write(unit=stdout,fmt='(A,I4,A,E30.20)') &
107 'For zonal wave number',m,' difference ',sqrt(num)/sqrt(den)
108
109 xx = sum( real(r_leg(1:xbx%nj))* real(r_leg(1:xbx%nj))+ &
110 aimag(r_leg(1:xbx%nj))*aimag(r_leg(1:xbx%nj)) )
111 write(unit=stdout,fmt='(a,i5,a,e30.20)') 'For Wave = ',m,' LX.LX = ',xx
112
113 ccv1(index_start:index_end) = 0.
114 call da_legtra_inv_adj( jds, jde, jts, jte, max_wave, xbx%alp_size, m, xbx%alp, &
115 ccv1(index_start:index_end), r_leg )
116 xx = sum( real(ccv(index_start:index_end))* &
117 real(ccv1(index_start:index_end)) +&
118 aimag(ccv(index_start:index_end))* &
119 aimag(ccv1(index_start:index_end)) )
120 write(unit=stdout,fmt='(a,i5,a,e30.20,/)') 'For Wave = ',m,' X.L^TLX = ',xx
121 end do
122
123 end subroutine da_test_spectral
124
125