var/da/da_ffts/da_solve_poissoneqn_fct.inc

   1 subroutine da_solve_poissoneqn_fct(grid, xbx, del2b, b)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Solve Del**2 B = A for B with zero gradient boundary conditions.
   5    !
   6    ! Method:  1) Compute spectral del2b using double forward FCT.
   7    !          2) Calculate spectral b.
   8    !          3) Reform gridpt. b using inverse double FCT.
   9    !          4) Remove mean b (arbitrary constant).
  10    !--------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    type(domain),     intent(inout) :: grid
  15    type(xbx_type),   intent(in)    :: xbx     ! Header & non-gridded vars.
  16    real,             intent(in)    :: del2b(ims:ime,jms:jme,kms:kme)   ! Del**2 B.
  17    real,             intent(out)   :: b(ims:ime,jms:jme,kms:kme)       ! B
  18
  19    integer           :: vector_inc    ! Increment between FST data.
  20    integer           :: vector_jump   ! Jump between start of vectors.
  21    integer           :: vector_size   ! Of form 2**p 3**q 5**r for FSTs.
  22    integer           :: num_vectors   ! Number of FSTs to perform.
  23    integer           :: work_area     ! Dimension for FST routine.
  24    integer           :: idim          ! Size of 1st dimension for FST.
  25    integer           :: jdim          ! Size of 2nd dimension for FST.
  26    integer           :: i, j, k, n, ij     ! loop counter
  27    real, allocatable :: work_1d(:)     ! FFT work array
  28    real              :: global_mean(kts:kte)
  29    real              :: local_mean(kts:kte)
  30    real              :: rij
  31
  32    if (trace_use) call da_trace_entry("da_solve_poissoneqn_fct")
  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 calculation of gridpoint b level by level:
  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
  75    ! [2.3] Perform forward FFT:
  76
  77    do k = grid%xp%ktsx, grid%xp%ktex
  78       ij = 0
  79       do j=grid%xp%jtsx, grid%xp%jtex
  80          do i=ids, ide
  81             ij=ij+1
  82             work_1d(ij) = grid%xp%v1x(i,j,k)
  83          end do
  84
  85          do i=1, xbx%fft_pad_i
  86             ij=ij+1
  87             work_1d(ij) = 0.0
  88          end do
  89       end do
  90       call fft551(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
  95       ij = 0
  96       do j=grid%xp%jtsx, grid%xp%jtex
  97          do i=ids, ide
  98             ij=ij+1
  99             grid%xp%v1x(i,j,k) = work_1d(ij)
 100          end do
 101
 102          do n=1, xbx%fft_pad_i
 103             i=(n-1)*xbx%pad_inc + 1
 104             ij=ij+1
 105             grid%xp%v2x(i,j,k) = work_1d(ij)
 106          end do
 107       end do
 108    end do
 109
 110    !---------------------------------------------------------------------------
 111    ! [3.0] For each k-level, perform forward FFT in y direction, apply spectral
 112    !       Poisson equation, and then perform inverse FFT in y direction:
 113    !---------------------------------------------------------------------------
 114
 115    ! [3.1] Apply (i,j',k' -> i',j,k') transpose (v1x -> v1y).
 116
 117    call da_transpose_x2y (grid)
 118    call da_transpose_x2y_v2 (grid)
 119
 120    ! [3.2] Set up FFT parameters:
 121
 122    idim = grid%xp%itey - grid%xp%itsy + 1 + xbx%pad_num
 123    jdim = xbx%fft_jy
 124
 125    vector_inc  = idim
 126    vector_jump = 1
 127    vector_size = jdim - 1
 128    num_vectors = idim
 129
 130    work_area   = (vector_size+1)*num_vectors
 131
 132
 133    ! [2.3] Perform forward FFT in j:
 134
 135    do k = grid%xp%ktsy, grid%xp%ktey
 136       ij = 0
 137       do j=jds, jde
 138          do i=grid%xp%itsy, grid%xp%itey
 139             ij=ij+1
 140             work_1d(ij) = grid%xp%v1y(i,j,k)
 141          end do
 142
 143          do n=1, xbx%pad_num
 144             i=xbx%pad_loc(n)
 145             ij=ij+1
 146             work_1d(ij) = grid%xp%v2y(i,j,k)
 147          end do
 148       end do
 149
 150       do j=1, xbx%fft_pad_j
 151          do i=grid%xp%itsy, grid%xp%itey+xbx%pad_num
 152             ij=ij+1
 153             work_1d(ij) = 0.0
 154          end do
 155       end do
 156
 157       call fft551(Forward_FFT, vector_inc, vector_jump, &
 158                                    num_vectors, vector_size, &
 159                                    xbx % fft_factors_y, xbx % trig_functs_y, &
 160                                    work_1d(1), work_area)
 161
 162       !------------------------------------------------------------------------
 163       ! [4.0] Solve spectral Poisson equation:
 164       !------------------------------------------------------------------------
 165
 166       ij = 0
 167       do j=grid%xp%jds, xbx%fft_jy
 168          do i=grid%xp%itsy, grid%xp%itey
 169             ij=ij+1
 170             work_1d(ij) = xbx%fft_coeffs(i,j)*work_1d(ij)
 171          end do
 172          do n=1, xbx%pad_num
 173             i=xbx%pad_pos(n)
 174             ij=ij+1
 175             work_1d(ij) = xbx%fft_coeffs(i,j)*work_1d(ij)
 176          end do
 177       end do
 178
 179       ! [2.3] Reform gridpt. b using inverse double FST in i.
 180
 181       call fft551(Inverse_FFT, vector_inc, vector_jump, &
 182                                      num_vectors, vector_size, &
 183                                    xbx % fft_factors_y, xbx % trig_functs_y, &
 184                                    work_1d(1), work_area)
 185
 186       ij = 0
 187       do j=jds, jde
 188          do i=grid%xp%itsy, grid%xp%itey
 189             ij=ij+1
 190             grid%xp%v1y(i,j,k) = work_1d(ij)
 191          end do
 192
 193          do n=1, xbx%pad_num
 194             i=xbx%pad_loc(n)
 195             ij=ij+1
 196             grid%xp%v2y(i,j,k) = work_1d(ij)
 197          end do
 198       end do
 199    end do
 200
 201    !---------------------------------------------------------------------------
 202    ! Perform inverse FFT in x direction:
 203    !---------------------------------------------------------------------------
 204
 205    ! Apply (i',j,k' -> i,j',k') transpose (v1y -> v1x).
 206
 207    call da_transpose_y2x (grid)
 208    call da_transpose_y2x_v2 (grid)
 209
 210    ! Set up FFT parameters:
 211
 212    idim = xbx%fft_ix
 213    jdim = grid%xp%jtex-grid%xp%jtsx+1
 214
 215    vector_inc  = 1
 216    vector_jump = idim
 217    vector_size = idim - 1
 218
 219    num_vectors = jdim
 220
 221    work_area   = (vector_size+1)*num_vectors
 222
 223    do k = grid%xp%ktsx, grid%xp%ktex
 224       ij = 0
 225       do j=grid%xp%jtsx, grid%xp%jtex
 226          do i=ids, ide
 227             ij=ij+1
 228             work_1d(ij) = grid%xp%v1x(i,j,k)
 229          end do
 230
 231          do n=1, xbx%fft_pad_i
 232             i=(n-1)*xbx%pad_inc + 1
 233             ij=ij+1
 234             work_1d(ij) = grid%xp%v2x(i,j,k)
 235          end do
 236       end do
 237
 238       call fft551(Inverse_FFT, vector_inc, vector_jump, &
 239                                      num_vectors, vector_size, &
 240                                      xbx % fft_factors_x, xbx % trig_functs_x, &
 241                                      work_1d(1), work_area)
 242       ij = 0
 243       do j=grid%xp%jtsx, grid%xp%jtex
 244          do i=ids, ide
 245             ij=ij+1
 246             grid%xp%v1x(i,j,k) = work_1d(ij)
 247          end do
 248
 249          ij=ij+xbx%fft_pad_i
 250       end do
 251    end do
 252
 253    ! Apply (i,j',k') -> i',j',k) transpose to restore v1z.
 254
 255    call da_transpose_x2z (grid)
 256
 257    ! Remove mean b (set arbitrary constant to zero):
 258
 259    rij = 1.0/real((ide-ids+1)*(jde-jds+1))
 260
 261    do k=kts, kte
 262       local_mean(k) = sum(grid%xp%v1z(its:ite,jts:jte,k))*rij
 263    end do
 264
 265    call wrf_dm_sum_reals (local_mean, global_mean)
 266
 267    do k=kts,kte
 268       write (unit=stdout,fmt=*) 'TEST_COVERAGE_DA_Solve_PoissonEqn_FCT:  global_mean(', &
 269          k,') = ', global_mean(k)
 270    end do
 271
 272    ! [2.5] Write data array into b:
 273
 274    do k=kts, kte
 275       b(its:ite,jts:jte,k) = grid%xp%v1z(its:ite,jts:jte,k) - global_mean(k)
 276    end do
 277
 278    !---------------------------------------------------------------------------
 279    ! [5.0] Tidy up:
 280    !---------------------------------------------------------------------------
 281
 282    if (allocated(work_1d)) deallocate(work_1d)
 283
 284    if (trace_use) call da_trace_exit("da_solve_poissoneqn_fct")
 285
 286 end subroutine da_solve_poissoneqn_fct
 287
 288