var/da/da_spectral/da_vtovv_spectral_adj.inc

   1 subroutine da_vtovv_spectral_adj(max_wavenumber, sizec, lenr, lenwrk, lensav, inc, &
   2                                     alp_size, alp, wsave, power, rcv, field)
   3
   4    !----------------------------------------------------------------------
   5    ! Purpose: Performs Adjoint of spectral to grid transformation on a sphere.
   6    !----------------------------------------------------------------------
   7
   8    implicit none
   9
  10    integer, intent(in) :: max_wavenumber             ! 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_wavenumber)    ! Power spectrum
  20    real, intent(out)   :: rcv(1:2*sizec)             ! Spectral modes.
  21    real, intent(in)    :: field(its:ite,jts:jte)     ! Gridpoint field.
  22
  23    integer             :: j, m, n                    ! Loop counters.
  24    integer             :: index_start                ! Position markers in cv.
  25 #ifdef DM_PARALLEL
  26    integer             :: index_end                  ! Position markers in cv.
  27 #endif
  28    integer             :: index_r, index_c           ! Array index for complex v_fou.
  29
  30    real                :: r_fou(1:lenr)              ! FFT array.
  31    complex             :: v_fou(its:ite,0:max_wavenumber)! Intermediate Fourier state.
  32    complex             :: ccv(1:sizec)               ! Spectral modes.
  33    complex             :: ccv_local(1:sizec)         ! Spectral modes.
  34
  35    integer              :: l, js, je           ! Loop counters.
  36    integer              :: index_m, index_j    ! Markers.
  37    complex              :: sum_legtra          ! Summation scalars.
  38
  39    integer              :: jc, iequator, temp
  40
  41 #ifdef FFTPACK
  42    real                :: work(1:lenwrk)             ! FFT work array.
  43 #endif
  44
  45    if (trace_use) call da_trace_entry("da_vtovv_spectral_adj")
  46
  47    !----------------------------------------------------------------------
  48    ! [1] Perform Adjoint of inverse Fourier decomposition in E-W direction:
  49    !----------------------------------------------------------------------
  50
  51    v_fou = 0.0
  52    do j = jts, jte
  53       r_fou(its:ite) = field(its:ite,j)
  54 #ifdef FFTPACK
  55       call rfft1f(ide, inc, r_fou, lenr, wsave, lensav, work, lenwrk, ierr)
  56 #else
  57       call da_error(__FILE__,__LINE__,(/"Must compile with FFTPACK"/))
  58 #endif
  59
  60       !----------------------------------------------------------------------
  61       ! Adjust the output for adjoint test
  62       !----------------------------------------------------------------------
  63       r_fou      =  real(ite)/2.0 * r_fou
  64       r_fou(its) =  r_fou(its)   * 2.0
  65
  66       ! if(.not. odd_longitudes) r_fou(ite) = 2.0*r_fou(ite)
  67       ! make r_fou(ide) zero as there is no power computed for this wavenumber
  68       r_fou(ite) = 0.0
  69
  70       v_fou(j,0) = CMPLX(r_fou(its), r_fou(ite))
  71
  72       do m = 1, max_wavenumber
  73          index_r = 2 * m
  74          index_c = 2 * m + 1
  75          v_fou(j,m)  = v_fou(j,m) + cmplx(r_fou(index_r),r_fou(index_c))
  76       end do
  77    end do
  78
  79    !----------------------------------------------------------------------
  80    ! [2] Perform adjoint of inverse Legendre decomposition in N-S direction:
  81    !----------------------------------------------------------------------
  82
  83    ccv_local(:) = 0.0
  84
  85    do m = 0, max_wavenumber
  86       index_start = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 +1
  87
  88       index_m     = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
  89
  90       jc = (jde-jds+1)/2
  91
  92       iequator = mod(jde-jds+1, 2)
  93
  94       js = max(jts, jc+iequator+1)
  95       je = min(jc+iequator, jte)
  96
  97       temp = (max_wavenumber + 1) * (max_wavenumber + 2) / 2
  98
  99       do l = m, max_wavenumber
 100          sum_legtra = da_zero_complex
 101
 102          if (mod(l+m,2) == 1) then
 103             do j = js, jte
 104                index_j = (jds+jde - j - 1) * temp
 105                sum_legtra = sum_legtra - v_fou(j,m) * alp(index_j + index_m + l)
 106             end do
 107          else
 108             do j = js, jte
 109                index_j = (jds+jde - j - 1) * temp
 110                sum_legtra = sum_legtra + v_fou(j,m) * alp(index_j + index_m + l)
 111             end do
 112          end if
 113
 114          do j = jts, je
 115             index_j = (j - 1) * temp
 116             sum_legtra = sum_legtra + v_fou(j,m) * alp(index_j + index_m + l)
 117          end do
 118
 119          ccv_local(index_start+l-m) = sum_legtra
 120       end do
 121    end do
 122
 123 #ifdef DM_PARALLEL
 124    index_start = 1
 125    index_end   = max_wavenumber + &
 126       max_wavenumber * (max_wavenumber + 1) / 2 + 1
 127
 128    n = index_end - index_start + 1
 129    call mpi_allreduce(ccv_local(index_start:index_end), &
 130                       ccv(index_start:index_end), n, true_mpi_complex, mpi_sum, &
 131                       comm, ierr)
 132 #else
 133    ccv(:) = ccv_local(:)
 134 #endif
 135
 136    !----------------------------------------------------------------------
 137    ! [2] Apply Power spectrum
 138    !-------------------------------------------------------------------------
 139
 140    if (.not. test_transforms) call da_apply_power(power, max_wavenumber, ccv, sizec)
 141
 142    do n = 1, sizec
 143       rcv(2*n - 1) = real (ccv(n))
 144       rcv(2*n    ) = aimag(ccv(n))
 145    end do
 146
 147    if (trace_use) call da_trace_exit("da_vtovv_spectral_adj")
 148
 149 end subroutine da_vtovv_spectral_adj
 150
 151