check_kt_grids_box Subroutine

  subroutine check_kt_grids_box(report_unit)
    use theta_grid, only: shat, is_effectively_zero_shear
    use constants, only: twopi, pi
    use warning_helpers, only: is_zero, not_exactly_equal
    implicit none
    integer, intent(in) :: report_unit
    real :: lx
    integer :: nx, ny, naky, ntheta0
    integer :: naky_calc, ntheta0_calc

    naky=naky_private
    ntheta0=ntheta0_private
    nx=nx_private
    ny=ny_private

    if (not_exactly_equal(y0, 2.)) then
       if (abs(twopi*y0 - ly) > 1.0e-7) then
          write (report_unit, *)
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('You cannot specify both ly and y0.')")
          write (report_unit, fmt="(' ly=',e12.4,'  2.0*pi*y0=',2e12.4)") ly
          write (report_unit, fmt="('THIS IS AN ERROR.')")
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *)
       end if
    end if

    write (report_unit, *)
    write (report_unit, fmt="('A rectangular simulation domain has been selected.')")
    if (rhostar_box > 0.0 .and. n0 > 0) write (report_unit, fmt="('The flux-tube size corresponds to toroidal mode number, n0=',i8/T44,' at rhostar_box=',1pe12.4)") n0,rhostar_box
    write (report_unit, *)
    write (report_unit, fmt="('The domain is ',f10.4,' rho in the y direction.')") ly


    if (is_effectively_zero_shear()) then
       if (is_zero(x0)) then
          if (rtwist > 0) then
             write (report_unit, fmt="('At theta=0, the domain has Lx = ',f10.5)")  ly*rtwist
          else
             write (report_unit, fmt="('At theta=0, the domain has Lx = ',f10.5)") -ly/rtwist
          end if
       else
          if (x0 > 0) then
             write (report_unit, fmt="('At theta=0, the domain has Lx = ',f10.5)")  twopi*x0
          else
             write (report_unit, fmt="('At theta=0, the domain has Lx = ',f10.5)") -twopi/x0
          end if
       end if
    else
       lx = ly * jtwist / (twopi * abs(shat))
       write (report_unit, fmt="('At theta=0, the domain is ',f10.4,' rho in the x direction.')") lx
       write (report_unit,*) ly, rtwist, jtwist, pi, shat
    end if

    write (report_unit, *)
    write (report_unit, fmt="('The nonlinear terms will be evaluated on a grid with ',&
         & i4,' points in x and ',i4,' points in y.')") nx, ny
    write (report_unit, *)
    naky_calc = (ny-1)/3+1
    ntheta0_calc = 2*((nx-1)/3)+1
    write (report_unit, fmt="('After de-aliasing, there will be ',i4,'  ky >= 0 modes and ',i4,' kx modes.')") naky_calc, ntheta0_calc
    write (report_unit, fmt="('The modes with ky < 0 are determined by the reality condition.')")

    if ( naky_calc < naky ) then
       write (report_unit, *)
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('ERROR : The requested number of ky >= 0 modes is ',I0,' which exceeds that available after de-aliasing.')") naky
       write (report_unit, fmt="('THIS IS AN ERROR.')")
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *)
    else if ( naky_calc > naky ) then
       write (report_unit, *)
       write (report_unit, fmt="('The requested number of ky >= 0 modes is ',I0,' which is less than available after de-aliasing.')") naky
       write (report_unit, fmt="('This corresponds to extra padding in the y grid of ',I0,' points.')") ny - ((naky - 1)*  3 + 1)
       write (report_unit, *)
    end if

    if ( ntheta0_calc < ntheta0 ) then
       write (report_unit, *)
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, fmt="('ERROR : The requested number of kx modes is ',I0,' which exceeds that available after de-aliasing.')") ntheta0
       write (report_unit, fmt="('THIS IS AN ERROR.')")
       write (report_unit, fmt="('################# WARNING #######################')")
       write (report_unit, *)
    else if ( ntheta0_calc > ntheta0 ) then
       write (report_unit, *)
       write (report_unit, fmt="('The requested number of kx modes is ',I0,' which is less than available after de-aliasing.')") ntheta0
       write (report_unit, fmt="('This corresponds to extra padding in the x grid of ',I0,' points.')") nx - (((ntheta0 - 1) /  2) * 3 + 1)
       write (report_unit, *)
    end if

  end subroutine check_kt_grids_box