var/da/da_spectral/da_vtovv_spectral.inc

   1 subroutine da_vtovv_spectral(max_wave, sizec, lenr, lenwrk, lensav, inc, &
   2                                 alp_size, alp, wsave, power, rcv, field)
   3
   4    !-------------------------------------------------------------------------
   5    ! Purpose: Performs spectral to gridpoint transformation on a sphere.
   6    !----------------------------------------------------------------------
   7
   8    implicit none
   9
  10    integer, intent(in) :: max_wave                 ! Max total wavenumber.
  11    integer, intent(in) :: sizec                    ! Size of packed spectral array.
  12    integer, intent(in) :: lenr                     ! FFT info.
  13    integer, intent(in) :: lenwrk                   ! FFT info.
  14    integer, intent(in) :: lensav                   ! FFT info.
  15    integer, intent(in) :: inc                      ! FFT info.
  16    integer, intent(in) :: alp_size                 ! Size of alp array.
  17    real, intent(in)    :: alp(1:alp_size)          ! Associated Legendre Polynomials.
  18    real, intent(in)    :: wsave(1:lensav)          ! Primes for FFT.
  19    real*8, intent(in)  :: power(0:max_wave)        ! Power spectrum
  20    real, intent(in)    :: rcv(1:2*sizec)           ! Spectral modes.
  21    real, intent(out)   :: field(its:ite,jts:jte)   ! Gridpoint field.
  22
  23    integer :: j, l,m, n                ! Loop counters.
  24    integer :: index_start              ! Position markers in cv.
  25    integer :: index_r, index_c         ! Array index for complex v_fou.
  26
  27    real    :: r_fou(1:lenr)            ! FFT array.
  28    complex :: v_fou(its:ite,0:max_wave)! Intermediate Fourier state.
  29    complex :: ccv(1:sizec)             ! Spectral modes.
  30    integer              :: index_m, index_j
  31    integer              :: jc, js, je, iequator
  32    complex              :: sum_legtra          ! Summation scalars.
  33
  34 #ifdef FFTPACK
  35    real    :: work(1:lenwrk)           ! FFT work array.
  36 #endif
  37
  38
  39    if (trace_use) call da_trace_entry("da_vtovv_spectral")
  40
  41    !----------------------------------------------------------------------------
  42    ! [1] Create complex array from read array:
  43    !----------------------------------------------------------------------------
  44
  45    v_fou = 0.0
  46    do n = 1, sizec
  47       ccv(n)  = CMPLX(rcv(2*n-1), rcv(2*n))
  48    end do
  49
  50    !----------------------------------------------------------------------------
  51    ! [2]  Apply power spectrum
  52    !----------------------------------------------------------------------------
  53
  54    if (.not. test_transforms) call da_apply_power(power, max_wave, ccv, sizec)
  55
  56    !----------------------------------------------------------------------------
  57    ! [3] Perform inverse Legendre decomposition in N-S direction:
  58    !----------------------------------------------------------------------------
  59
  60    do m = 0, max_wave
  61       index_start = m * (max_wave + 1 - m) + m * (m + 1) / 2 + 1
  62
  63       index_m = m * (max_wave + 1 - m) + m * (m + 1) / 2 + 1 - m
  64
  65       jc = (jde-jds+1)/2
  66
  67       iequator = mod(jde-jds+1, 2)
  68
  69       je = min(jc+iequator, jte)
  70
  71       do j = jts, je
  72          index_j = (j - 1) * (max_wave + 1) * (max_wave + 2) / 2
  73
  74          v_fou(j,m) = sum(ccv(index_start:index_start-m+max_wave) * &
  75             alp(index_j+index_m+m:index_j+index_m+max_wave))
  76       end do
  77
  78       js = max(jts, jc+iequator+1)
  79
  80       do j = js, jte
  81          index_j = (jds+jde - j - 1) * (max_wave + 1) * (max_wave + 2) / 2
  82
  83          sum_legtra = da_zero_complex
  84          do l = m, max_wave
  85             ! Calculate second quadrant values:
  86             if(mod(l+m,2) == 1) then
  87                sum_legtra = sum_legtra - ccv(index_start-m+l) * alp(index_j + index_m + l)
  88             else
  89                sum_legtra = sum_legtra + ccv(index_start-m+l) * alp(index_j + index_m + l)
  90             end if
  91          end do
  92          v_fou(j,m) = sum_legtra
  93       end do
  94    end do
  95
  96    !----------------------------------------------------------------------------
  97    ! [4] Perform inverse Fourier decomposition in E-W direction:
  98    !----------------------------------------------------------------------------
  99
 100    do j = jts, jte
 101       r_fou(its) =  real(v_fou(j,0))  ! R(m=0) is real.
 102       ! r_fou(ite) = aimag(v_fou(j,0)) ! R(m=NI/2) is real, but packed in imag m = 0)
 103       ! make r_fou(ide) zero as there is no power computed corresponding to this wavenumber
 104       r_fou(ite) = 0.0
 105
 106       do m = 1, max_wave
 107          index_r = 2 * m
 108          index_c = 2 * m + 1
 109          r_fou(index_r) = real(v_fou(j,m))
 110          r_fou(index_c) = aimag(v_fou(j,m))
 111       end do
 112
 113 #ifdef FFTPACK
 114       call rfft1b(ide, inc, r_fou, lenr, wsave, lensav, work, lenwrk, ierr)
 115 #else
 116       call da_error(__FILE__,__LINE__,(/"Must compile with FFTPACK"/))
 117 #endif
 118       field(its:ite,j) = r_fou(its:ite)
 119    end do
 120
 121    if (trace_use) call da_trace_exit("da_vtovv_spectral")
 122
 123 end subroutine da_vtovv_spectral
 124
 125