var/da/da_ffts/da_solve_poissoneqn_fst.inc

   1 subroutine da_solve_poissoneqn_fst(grid, xbx, del2b, b)
   2
   3    !--------------------------------------------------------------------------
   4    !  Purpose: Solve Del**2 B = A for B with zero field boundary conditions.
   5    !
   6    !  Method:  1) Compute spectral del2b using double forward FST.
   7    !           2) Calculate spectral b.
   8    !           3) Reform gridpt. b using inverse double FST.
   9    !           Note no mean removed (would make b.ne.0 at boundaries) so
  10    !           arbitrary constant still present.
  11    !--------------------------------------------------------------------------
  12
  13    implicit none
  14
  15    type (domain),   intent(inout) :: grid
  16    type (xbx_type), intent(in)    :: xbx           ! Header & non-gridded vars.
  17    real,            intent(in)    :: del2b(ims:ime,jms:jme,kms:kme)  ! Del**2 B.
  18    real,            intent(out)   :: b(ims:ime,jms:jme,kms:kme)  ! B.
  19
  20    integer           :: vector_inc    ! Increment between FST data.
  21    integer           :: vector_jump   ! Jump between start of vectors.
  22    integer           :: vector_size   ! Of form 2**p 3**q 5**r for FSTs.
  23    integer           :: num_vectors   ! Number of FSTs to perform.
  24    integer           :: work_area     ! Dimension for FST routine.
  25    integer           :: idim          ! Size of 1st dimension for FST.
  26    integer           :: jdim          ! Size of 2nd dimension for FST.
  27
  28    integer           :: i, j, k, n, ij     ! loop counter
  29
  30    real, allocatable :: work_1d(:)     ! FFT work array
  31
  32    if (trace_use) call da_trace_entry("da_solve_poissoneqn_fst")
  33
  34    !------------------------------------------------------------------------------
  35    ! [1.0] Initialise:
  36    !------------------------------------------------------------------------------
  37
  38    ! Calculate work space needed.
  39
  40    n = max(xbx%fft_ix*(grid%xp%jtex-grid%xp%jtsx+1), &
  41            xbx%fft_jy*(grid%xp%itey-grid%xp%itsy+1+xbx%pad_num))
  42
  43    ! Allocate work arrays.
  44    allocate(work_1d(1:n))
  45
  46    ! Copy del2b for transpose.
  47
  48    grid%xp%v1z(its:ite,jts:jte,kts:kte) = del2b(its:ite,jts:jte,kts:kte)
  49
  50    if (ite == ide) grid%xp%v1z(ite,jts:jte,kts:kte) = 0.0
  51    if (jte == jde) grid%xp%v1z(its:ite,jte,kts:kte) = 0.0
  52
  53    !---------------------------------------------------------------------------
  54    ! [2.0] Perform forward FFT in x direction:
  55    !---------------------------------------------------------------------------
  56
  57    ! [2.1] Apply (i',j',k -> i,j',k') transpose (v1z -> v1x).
  58
  59    call da_transpose_z2x (grid)
  60
  61    ! [2.2] Set up FFT parameters:
  62
  63    idim = xbx%fft_ix
  64    jdim = grid%xp%jtex-grid%xp%jtsx+1
  65
  66    vector_inc  = 1
  67    vector_jump = idim
  68    vector_size = idim - 1
  69
  70    num_vectors = jdim
  71
  72    work_area   = (vector_size+1)*num_vectors
  73
  74    ! [2.3] Perform forward FFT:
  75
  76    do k = grid%xp%ktsx, grid%xp%ktex
  77       ij = 0
  78       do j=grid%xp%jtsx, grid%xp%jtex
  79          do i=ids, ide
  80             ij=ij+1
  81             work_1d(ij) = grid%xp%v1x(i,j,k)
  82          end do
  83
  84          do i=1, xbx%fft_pad_i
  85             ij=ij+1
  86             work_1d(ij) = 0.0
  87          end do
  88       end do
  89
  90       call fft661(Forward_FFT, vector_inc, vector_jump, &
  91          num_vectors, vector_size, &
  92          xbx % fft_factors_x, xbx % trig_functs_x, &
  93          work_1d(1), work_area)
  94       ij = 0
  95       do j=grid%xp%jtsx, grid%xp%jtex
  96          do i=ids, ide
  97             ij=ij+1
  98             grid%xp%v1x(i,j,k) = work_1d(ij)
  99          end do
 100
 101          do n=1, xbx%fft_pad_i
 102             i=(n-1)*xbx%pad_inc + 1
 103             ij=ij+1
 104             grid%xp%v2x(i,j,k) = work_1d(ij)
 105          end do
 106       end do
 107    end do
 108
 109    !---------------------------------------------------------------------------
 110    ! [3.0] For each k-level, perform forward FFT in y direction, apply spectral
 111    !        Poisson equation, and then perform inverse FFT in y direction:
 112    !---------------------------------------------------------------------------
 113
 114    ! [3.1] Apply (i,j',k' -> i',j,k') transpose (v1x -> v1y).
 115
 116    call da_transpose_x2y (grid)
 117    call da_transpose_x2y_v2 (grid)
 118
 119    ! [3.2] Set up FFT parameters:
 120
 121    idim = grid%xp%itey - grid%xp%itsy + 1 + xbx%pad_num
 122    jdim = xbx%fft_jy
 123
 124    vector_inc  = idim
 125    vector_jump = 1
 126    vector_size = jdim - 1
 127    num_vectors = idim
 128
 129    work_area   = (vector_size+1)*num_vectors
 130
 131    ! [2.3] Perform forward FFT in j:
 132
 133    do k = grid%xp%ktsy, grid%xp%ktey
 134       ij = 0
 135       do j=jds,jde
 136          do i=grid%xp%itsy, grid%xp%itey
 137             ij=ij+1
 138             work_1d(ij) = grid%xp%v1y(i,j,k)
 139          end do
 140
 141          do n=1, xbx%pad_num
 142             i=xbx%pad_loc(n)
 143             ij=ij+1
 144             work_1d(ij) = grid%xp%v2y(i,j,k)
 145          end do
 146       end do
 147
 148       do j=1, xbx%fft_pad_j
 149          do i=grid%xp%itsy, grid%xp%itey+xbx%pad_num
 150             ij=ij+1
 151             work_1d(ij) = 0.0
 152          end do
 153       end do
 154
 155       call fft661(Forward_FFT, vector_inc, vector_jump, &
 156          num_vectors, vector_size, &
 157          xbx % fft_factors_y, xbx % trig_functs_y, &
 158          work_1d(1), work_area)
 159
 160       !------------------------------------------------------------------------
 161       ! [4.0] Solve spectral Poisson equation:
 162       !------------------------------------------------------------------------
 163
 164       ij = 0
 165       do j=jds, xbx%fft_jy
 166          do i=grid%xp%itsy, grid%xp%itey
 167             ij=ij+1
 168             work_1d(ij) = xbx%fft_coeffs(i,j)*work_1d(ij)
 169          end do
 170          do n=1, xbx%pad_num
 171             i=xbx%pad_pos(n)
 172             ij=ij+1
 173             work_1d(ij) = xbx%fft_coeffs(i,j)*work_1d(ij)
 174          end do
 175       end do
 176
 177       ! [2.3] Reform gridpt. b using inverse double FST in i.
 178
 179       call fft661(Inverse_FFT, vector_inc, vector_jump, &
 180          num_vectors, vector_size, &
 181          xbx % fft_factors_y, xbx % trig_functs_y, &
 182          work_1d(1), work_area)
 183
 184       ij = 0
 185       do j=jds, jde
 186          do i=grid%xp%itsy, grid%xp%itey
 187             ij=ij+1
 188             grid%xp%v1y(i,j,k) = work_1d(ij)
 189          end do
 190
 191          do n=1, xbx%pad_num
 192             i=xbx%pad_loc(n)
 193             ij=ij+1
 194             grid%xp%v2y(i,j,k) = work_1d(ij)
 195          end do
 196       end do
 197    end do
 198
 199    !---------------------------------------------------------------------------
 200    ! Perform inverse FFT in x direction:
 201    !---------------------------------------------------------------------------
 202
 203    ! Apply (i',j,k' -> i,j',k') transpose (v1y -> v1x).
 204
 205    call da_transpose_y2x (grid)
 206    call da_transpose_y2x_v2 (grid)
 207
 208    ! Set up FFT parameters:
 209
 210    idim = xbx%fft_ix
 211    jdim = grid%xp%jtex-grid%xp%jtsx+1
 212
 213    vector_inc  = 1
 214    vector_jump = idim
 215    vector_size = idim - 1
 216
 217    num_vectors = jdim
 218
 219    work_area   = (vector_size+1)*num_vectors
 220
 221    do k = grid%xp%ktsx, grid%xp%ktex
 222       ij = 0
 223       do j=grid%xp%jtsx, grid%xp%jtex
 224          do i=grid%xp%ids, grid%xp%ide
 225             ij=ij+1
 226             work_1d(ij) = grid%xp%v1x(i,j,k)
 227          end do
 228
 229          do n=1, xbx%fft_pad_i
 230             i=(n-1)*xbx%pad_inc + 1
 231             ij=ij+1
 232             work_1d(ij) = grid%xp%v2x(i,j,k)
 233          end do
 234       end do
 235
 236       call fft661(Inverse_FFT, vector_inc, vector_jump, &
 237          num_vectors, vector_size, &
 238          xbx % fft_factors_x, xbx % trig_functs_x, &
 239          work_1d(1), work_area)
 240       ij = 0
 241       do j=grid%xp%jtsx, grid%xp%jtex
 242          do i=ids, ide
 243             ij=ij+1
 244             grid%xp%v1x(i,j,k) = work_1d(ij)
 245          end do
 246
 247          ij=ij+xbx%fft_pad_i
 248       end do
 249    end do
 250
 251    ! Apply (i,j',k') -> i',j',k) transpose to restore v1z.
 252
 253    call da_transpose_x2z (grid)
 254
 255    !---------------------------------------------------------------------------
 256    ! [5.0] Tidy up:
 257    !---------------------------------------------------------------------------
 258
 259    if (allocated(work_1d)) deallocate(work_1d)
 260
 261    ! [2.5] Write data array into b:
 262
 263    b(its:ite,jts:jte,kts:kte) = grid%xp%v1z(its:ite,jts:jte,kts:kte)
 264
 265    if (trace_use) call da_trace_exit("da_solve_poissoneqn_fst")
 266
 267 end subroutine da_solve_poissoneqn_fst
 268
 269