var/da/da_setup_structures/da_setup_runconstants.inc

   1 subroutine da_setup_runconstants(grid, xbx)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Define constants used.
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (domain),   intent(inout) :: grid
  10    type (xbx_type), intent(inout) :: xbx     ! Header & non-gridded vars.
  11
  12    integer                  :: n             ! Loop counter.
  13
  14    integer                  :: fft_size_i    ! Efficient FFT 1st dimension.
  15    integer                  :: fft_size_j    ! Efficient FFT 2nd dimension.
  16    logical                  :: found_magic   ! true when efficient FFT found
  17    logical                  :: need_pad      ! True when need pad_i > 0
  18
  19    integer                  :: i, j
  20    integer                  :: nx            ! nx + 1 = ix + pad_i
  21    integer                  :: ny            ! ny + 1 = jy + pad_j
  22    real                     :: const         ! Multiplicative constants.
  23    real                     :: coeff_nx      ! Multiplicative coefficients.
  24    real                     :: coeff_ny      ! Multiplicative coefficients.
  25    real                     :: cos_coeff_nx  ! Multiplicative coefficients.
  26    real                     :: cos_coeff_ny  ! Multiplicative coefficients.
  27
  28    if (trace_use) call da_trace_entry("da_setup_runconstants")
  29
  30    !---------------------------------------------------------------------------
  31    ! [1.0] Calculate padded grid-size required for efficient FFTs and factors:
  32    !---------------------------------------------------------------------------
  33
  34    ! [1.1] In x direction:
  35
  36    nx = ide - ids
  37
  38    allocate(xbx % fft_factors_x(num_fft_factors))
  39
  40    do n = nx, nx+nrange
  41       call da_find_fft_factors(n, found_magic, xbx % fft_factors_x)
  42
  43       if (found_magic) then
  44          if (mod(n, 2) == 0) then
  45             fft_size_i = n
  46             xbx % fft_pad_i = n - nx
  47             exit
  48          end if
  49       end if
  50    end do
  51
  52    if (.NOT. found_magic) then
  53      call da_error(__FILE__,__LINE__, &
  54        (/"No FFT factor found: invalid e_we value"/))
  55    end if
  56
  57    allocate(xbx % trig_functs_x(1:3*fft_size_i))
  58
  59    call da_find_fft_trig_funcs(fft_size_i, xbx % trig_functs_x(:))
  60
  61    ! [1.2] In y direction:
  62
  63    ny = jde - jds
  64
  65    allocate(xbx % fft_factors_y(num_fft_factors))
  66
  67    do n = ny, ny+nrange
  68       call da_find_fft_factors(n, found_magic, xbx % fft_factors_y)
  69
  70       if (found_magic) then
  71          if (mod(n, 2) == 0) then
  72             fft_size_j = n
  73             xbx % fft_pad_j = n - ny
  74             exit
  75          end if
  76       end if
  77    end do
  78
  79    if (.NOT. found_magic) then
  80      call da_error(__FILE__,__LINE__, &
  81        (/"No FFT factor found: invalid e_sn value"/))
  82    end if
  83
  84    allocate(xbx % trig_functs_y(1:3*fft_size_j))
  85
  86    call da_find_fft_trig_funcs(fft_size_j, xbx % trig_functs_y(:))
  87
  88    !-----Multiplicative coefficent for solution of spectral Poission eqn:
  89
  90    !mslee.Bgrid
  91    ! const = -0.5 * grid%xb%ds * grid%xb%ds
  92    const = -2.0 * grid%xb%ds * grid%xb%ds
  93
  94    nx = ide - ids + xbx%fft_pad_i
  95    ny = jde - jds + xbx%fft_pad_j
  96    ! YRG: A-grid:
  97    coeff_nx = 2.0 * pi / real(nx)
  98    coeff_ny = 2.0 * pi / real(ny)
  99
 100    ! YRG: B-grid::
 101    ! coeff_nx = pi / real(nx)
 102    ! coeff_ny = pi / real(ny)
 103
 104    xbx%fft_ix = nx + 1
 105    xbx%fft_jy = ny + 1
 106
 107    allocate(xbx%fft_coeffs(1:xbx%fft_ix,1:xbx%fft_jy))
 108
 109    do j = 2, ny
 110       cos_coeff_ny = COS(coeff_ny * real(j-1))
 111       do i = 2, nx
 112          cos_coeff_nx = COS(coeff_nx * real(i-1))
 113          !mslee.Bgrid
 114          ! xbx%fft_coeffs(i,j) = const / (1.0 - cos_coeff_nx * cos_coeff_ny)
 115          if (cos_coeff_nx.eq.1.and.cos_coeff_ny.eq.1) then
 116             xbx%fft_coeffs(i,j) = 0.0
 117          else
 118             xbx%fft_coeffs(i,j) = const / (2.0 - cos_coeff_nx - cos_coeff_ny)
 119          end if
 120       end do
 121    end do
 122
 123    ! Set to zero all coefficients which are multiplied by sin(0.0) in FST:
 124
 125    !mslee      xbx%fft_coeffs(1,1:xbx%fft_jy)  = 0.0
 126    !mslee      xbx%fft_coeffs(xbx%fft_ix,1:xbx%fft_jy) = 0.0
 127    !mslee      xbx%fft_coeffs(1:xbx%fft_ix,1)  = 0.0
 128    !mslee      xbx%fft_coeffs(1:xbx%fft_ix,xbx%fft_jy) = 0.0
 129    !mslee. we need to check the following
 130
 131    xbx%fft_adjoint_factor = 4.0 / real(nx * ny)
 132
 133    !---------------------------------------------------------------------------
 134
 135    ! Calculate i increment for distributing pad region across processors.
 136
 137    need_pad = .true.
 138    if (xbx%fft_pad_i <= 0) then
 139       need_pad = .false.
 140    else if (xbx%fft_pad_i > (ide-ids+1)) then
 141       write(unit=message(1), fmt='(a)') &
 142        'FFT xbx%fft_pad_i is too large!'
 143       write(unit=message(2), fmt='(2x,a,i4)') &
 144          '(ide-ids+1) = ', (ide-ids+1), &
 145          'xbx%fft_pad_i = ', xbx%fft_pad_i
 146       call da_error (__FILE__,__LINE__,message(1:2))
 147    else
 148       xbx%pad_inc = (ide-ids+1)/xbx%fft_pad_i
 149    end if
 150
 151    ! Calculate number of pad vectors in x to be done on this processor.
 152    ! Need to save xbx%pad_num, xbx%pad_inc, and xbx%pad_loc
 153
 154    xbx%pad_num = 0
 155    if (need_pad) then
 156       do n=1, xbx%fft_pad_i
 157          i = (n-1)*xbx%pad_inc + 1
 158          if ((i >= grid%xp%itsy) .and. (i <= grid%xp%itey)) then
 159             xbx%pad_num = xbx%pad_num + 1
 160          end if
 161       end do
 162
 163       if (xbx%pad_num > 0) then
 164          allocate(xbx%pad_loc(1:xbx%pad_num))
 165          allocate(xbx%pad_pos(1:xbx%pad_num))
 166       end if
 167
 168       xbx%pad_num = 0
 169       do n=1, xbx%fft_pad_i
 170          i = (n-1)*xbx%pad_inc + 1
 171          if ((i >= grid%xp%itsy) .and. (i <= grid%xp%itey)) then
 172             xbx%pad_num = xbx%pad_num + 1
 173             xbx%pad_loc(xbx%pad_num) = i
 174             xbx%pad_pos(xbx%pad_num) = grid%xp%ide + n
 175          end if
 176       end do
 177    end if
 178
 179    !---------------------------------------------------------------------------
 180
 181    if (global) then
 182       ! Set up Spectral transform constants:
 183       xbx%inc    = 1
 184       xbx%ni     = ide-ids+1
 185       xbx%nj     = jde-jds+1
 186       xbx%lenr   = xbx%inc * (xbx%ni - 1) + 1
 187       xbx%lensav = xbx%ni + int(log(real(xbx%ni))) + 4
 188       xbx%lenwrk = xbx%ni
 189       xbx%max_wavenumber = xbx%ni/2 -1
 190       xbx%alp_size = (xbx%nj+ 1)*(xbx%max_wavenumber+1)*(xbx%max_wavenumber+2)/4
 191
 192       if (print_detail_spectral) then
 193          write (unit=stdout,fmt='(a)') "Spectral transform constants"
 194          write (unit=stdout, fmt='(a, i8)') &
 195             '   inc            =', xbx%inc   , &
 196             '   ni             =', xbx%ni    , &
 197             '   nj             =', xbx%nj    , &
 198             '   lenr           =', xbx%lenr  , &
 199             '   lensav         =', xbx%lensav, &
 200             '   lenwrk         =', xbx%lenwrk, &
 201             '   max_wavenumber =', xbx%max_wavenumber, &
 202             '   alp_size       =', xbx%alp_size
 203          write (unit=stdout,fmt='(a)') " "
 204       end if
 205
 206       allocate(xbx%coslat(1:xbx%nj))
 207       allocate(xbx%sinlat(1:xbx%nj))
 208       allocate(xbx%coslon(1:xbx%ni))
 209       allocate(xbx%sinlon(1:xbx%ni))
 210       allocate(xbx%int_wgts(1:xbx%nj))      ! Interpolation weights
 211       allocate(xbx%alp(1:xbx%alp_size))
 212       allocate(xbx%wsave(1:xbx%lensav))
 213
 214       call da_initialize_h(xbx%ni, xbx%nj, xbx%max_wavenumber, &
 215                            xbx%lensav, xbx%alp_size,           &
 216                            xbx%wsave, grid%xb%lon, xbx%sinlon,      &
 217                            xbx%coslon, grid%xb%lat, xbx%sinlat,     &
 218                            xbx%coslat, xbx%int_wgts, xbx%alp)
 219
 220
 221    end if
 222
 223    if (trace_use) call da_trace_exit("da_setup_runconstants")
 224
 225 end subroutine da_setup_runconstants
 226
 227