var/da/da_spectral/da_test_spectral.inc

   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    real,            intent(in) :: field(its:ite,jts:jte)    ! Gridpoint field.
  13
  14    real    :: field_out(its:ite,jts:jte)
  15    real*8  :: power(0:max_wave)   ! Power spectrum
  16    real    :: rcv(1:2*sizec)      ! Spectral modes.
  17    real    :: rcv_out(1:2*sizec)  ! Spectral modes.
  18    integer :: m,mm, index_start, index_end
  19    complex :: r_leg(1:jde)
  20    complex :: ccv(1:sizec)     ! Spectral modes.
  21    complex :: ccv1(1:sizec)    ! Spectral modes.
  22    real    :: den, num, xx
  23
  24    if (trace_use) call da_trace_entry("da_test_spectral")
  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.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.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.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.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.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    if (trace_use) call da_trace_exit("da_test_spectral")
 124
 125 end subroutine da_test_spectral
 126
 127