var/da/da_spectral/da_vv_to_v_spectral.inc

   1 subroutine da_vv_to_v_spectral(ni, nj, max_wavenumber, inc, lenr, lensav, lenwrk, &
   2                                 alp_size, r_cvsize, alp, wsave, int_wgts, rcv, field)
   3
   4    !-------------------------------------------------------------------------
   5    ! Purpose: Performs gridpoint to spectral transformation on a sphere.
   6    ! Note: Routine works for both regular and Gaussian latitude (latitude
   7    ! integration weights contain necessary info).
   8    !-------------------------------------------------------------------------
   9
  10    implicit none
  11
  12    integer, intent(in) :: ni               ! Number of longitudes.
  13    integer, intent(in) :: nj               ! Number of latitudes.
  14    integer, intent(in) :: r_cvsize         ! Size of real control cv-array.
  15    integer, intent(in) :: max_wavenumber   ! Smallest scale required (ni/2 - 1).
  16    integer, intent(in) :: inc              ! Jump between elements of vector in array.
  17    integer, intent(in) :: lenr             ! FFT array dimension (at least inc*(n-1)+1).
  18    real, intent(in)    :: field(1:ni,1:nj) ! Gridpoint field.
  19    real, intent(out)   :: rcv(1:r_cvsize)  ! Spectral modes.
  20    integer, intent(in) :: lensav           ! wsave dimension (n+int(log(real(ni)))+4).
  21    integer, intent(in) :: lenwrk           ! Dimension of work array.
  22    integer, intent(in) :: alp_size         ! Size of ALP vector.
  23    real, intent(in)    :: alp(1:alp_size)  ! Associated Legendre Polynomials.
  24    real, intent(in)    :: wsave(1:lensav)  ! Primes for FFT.
  25    real, intent(in)    :: int_wgts(1:nj)   ! Legendre integration weights.
  26
  27    integer             :: i, j, m, mm            ! Loop counters.
  28    integer             :: sizec                  ! Size of complex cv-array
  29    integer             :: index_r, index_c       ! Array index for complex v_fou
  30    integer             :: index_start, index_end ! Position markers in cv.
  31    real                :: r_fou(1:lenr)          ! FFT array.
  32    logical             :: odd_longitudes
  33    complex             :: v_fou(1:nj,0:max_wavenumber)! Intermediate Fourier state.
  34    complex             :: r_leg(1:nj)                 ! Intermediate Fourier state.
  35    complex, allocatable:: ccv(:)                      ! Spectral modes.
  36
  37 #ifdef FFTPACK
  38    real                :: work(1:lenwrk)         ! FFT work array.
  39 #endif
  40
  41 !   if (trace_use) call da_trace_entry("da_vv_to_v_spectral")
  42
  43    sizec = int(0.5 * r_cvsize)
  44    allocate (ccv(1:sizec))
  45
  46    if ((ni+1) / 2 == ni/2 + 1) then  ! Odd number of longitudes:
  47       odd_longitudes = .true.
  48    else                                ! Even number of longitudes:
  49       odd_longitudes = .false.
  50    end if
  51
  52    !-------------------------------------------------------------------------
  53    ! [1] Perform Adjoint of inverse Fourier decomposition in E-W direction:
  54    !-------------------------------------------------------------------------
  55
  56    if ((ni+1) / 2 == ni/2 + 1) then  ! Odd number of longitudes:
  57       odd_longitudes = .true.
  58    else                                ! Even number of longitudes:
  59       odd_longitudes = .false.
  60    end if
  61
  62    ! [1] Perform Fourier decomposition in E-W direction:
  63
  64    do j = 1, nj
  65       r_fou(1:ni) = field(1:ni,j)
  66 #ifdef FFTPACK
  67       call rfft1f(ni, inc, r_fou, lenr, wsave, lensav, work, lenwrk, ierr)
  68 #else
  69       call da_error(__FILE__,__LINE__,(/"Needs to be compiled with FFTPACK"/))
  70 #endif
  71
  72       if (odd_longitudes) then
  73          v_fou(j,0) = CMPLX(r_fou(1), 0.0) ! m = 0 is real.
  74       else
  75          ! m = 0 is real, but pack R(NI/2) in imag m = 0:
  76          v_fou(j,0) = CMPLX(r_fou(1), r_fou(ni))
  77       end if
  78
  79       do m = 1, max_wavenumber
  80          index_r = 2 * m
  81          index_c = 2 * m + 1
  82          v_fou(j,m) = CMPLX(r_fou(index_r), r_fou(index_c)) ! 2.0 * Fourier mode.
  83       end do
  84    end do
  85
  86    ! [2] Perform Legendre decomposition in N-S direction:
  87
  88    do m = 0, max_wavenumber
  89       index_start = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1
  90       index_end   = index_start + max_wavenumber - m
  91       r_leg(1:nj) = v_fou(1:nj,m)
  92       call da_legtra (nj, max_wavenumber, alp_size, m, int_wgts, alp, r_leg, &
  93                         ccv(index_start:index_end))
  94    end do
  95
  96    do i=1,sizec
  97       mm = 2*i - 1
  98       rcv(mm ) = real (ccv(i))
  99       rcv(mm+1) = aimag(ccv(i))
 100    end do
 101    deallocate (ccv)
 102
 103 !   if (trace_use) call da_trace_exit("da_vv_to_v_spectral")
 104
 105 end subroutine da_vv_to_v_spectral
 106
 107