init_kt_grids_box Subroutine

  subroutine init_kt_grids_box(kt_grids_box_config_in)
!CMR, 14/10/2013: 
! New namelist variables: n0, rhostar_box. 
! If n0 and rhostar_box defined, set ky(1) using toroidal mode number.

    use file_utils, only: error_unit
    use theta_grid, only: shat, drhodpsi
    use constants, only: pi
    use mp, only: mp_abort, proc0
    implicit none
    type(kt_grids_box_config_type), intent(in), optional :: kt_grids_box_config_in

    if (initialized) return
    initialized = .true.
    
    call read_parameters_box(kt_grids_box_config_in)

    if (ny==0 .and. naky==0) call mp_abort("ERROR: ny==0 .and. naky==0", .true.) 
    if (nx==0 .and. ntheta0==0) call mp_abort("ERROR: nx==0 .and. ntheta0==0", .true.) 

    if (rhostar_box .gt. 0.0 .and. n0 .gt. 0) y0=1.0/(n0*rhostar_box*drhodpsi)

    if (gryfx) then
      !WHEN RUNNING IN GRYFX, ONLY EVOLVE ky=0 MODES.
      naky = 1
      ny = 1

      !NEED TO ACCOUNT FOR SQRT(2) DIFFERENCE BETWEEN rho_GS2 and rho_GryfX, so
      !change x0 and y0
      y0 = y0/sqrt(2.)  
      !this still needs to be renormalized even though only ky=0
      !is running because sometimes x0 is set from y0 and jtwist later.
      x0 = x0/sqrt(2.) 
      !this is in case x0 is not set from y0 and jtwist later.
    end if


    if (y0 < 0) y0 = -1./y0

    !EGH This line does not affect any existing gs2 runs but is
    !here because gryfx uses jtwist < 0 to signal using the
    ! default
    if (jtwist < 0) jtwist = max(int(2.0*pi*shat + 0.5),1)

    if (ly == 0.) ly = 2.0*pi*y0
    if (naky == 0) naky = (ny-1)/3 + 1
    if (ntheta0 == 0) ntheta0 = 2*((nx-1)/3) + 1
    if (rtwist == 0.) rtwist = real(jtwist)

    if (mod(ntheta0, 2) /= 1) then
       call mp_abort("ERROR: ntheta0 must be an odd number in box mode", .true.)
    end if

    ! Now we make sure that we set ny and nx for given 
    ! choices of naky and ntheta0. If e.g. both ny and naky
    ! are set and they are not consistent with each other
    ! raise an error.

    if (ny == 0) then
       ! If ny hasn't been set the determine it from naky
       ny = (naky - 1)*  3 + 1
       if (proc0) write (error_unit(), '("INFO: ny (",I0,") set from naky (",I0,").")') ny, naky
    else
       ! If both naky and ny are set then check that the resulting padding is at least
       ! as much as required for dealiasing
       if (naky < (ny-1)/3 + 1) then ! Excess padding
          if (proc0) then
             write (error_unit(), '("INFO: Both ny (",I0,") and naky (",I0,") have been set by the user.")') ny, naky
             write (error_unit(), '("      these values lead to excess padding (",I0,".) for the FFTs.")') ny - ((naky - 1)*  3 + 1)
             write (error_unit(), '("      This may be desirable if the values chosen increase the number of parallelization sweetspots,")')
             write (error_unit(), '("      but this also makes the resolution in ny larger than is strictly necessary.")')
          end if

       else if (naky > (ny-1)/3 + 1) then ! Insufficient padding
          if (proc0) then
             write (error_unit(), '("ERROR: ny (",I0,") and naky (",I0,") have been set by the user.")') ny, naky
             write (error_unit(), '("       but these values lead to insufficient padding (",I0,") for the FFTS.")') ny - ((naky - 1)*  3 + 1)
             write (error_unit(), '("       The zero-padding must satisfy the one thirds rule to avoid aliasing, which requires naky <= (ny-1)/3 + 1.")')
             write (error_unit(), '("  Please do one of the following:")')
             write (error_unit(), '("       1. Increase ny or decrease naky to ensure naky <= (ny-1)/3 + 1")')
             write (error_unit(), '("       2. Set just one of ny and naky, GS2 will then set the other appropriately.")')
          end if

          call mp_abort("ERROR: naky and ny both set resulting in insufficient padding. See error file for more details.", .true.)
       end if
    end if

    if (nx == 0) then
       ! If nx hasn't been set the determine it from ntheta0
       nx = ((ntheta0 - 1) /  2) * 3 + 1
       if (proc0) write (error_unit(), '("INFO: nx (",I0,") set from ntheta0 (",I0,").")') nx, ntheta0
    else
       ! If both ntheta0 and nx are set then check that the resulting padding is at least
       ! as much as required for dealiasing
       if (ntheta0 < 2*((nx-1)/3) + 1) then ! Excess padding
          if (proc0) then
             write (error_unit(), '("INFO: Both nx (",I0,") and ntheta0 (",I0,") have been set by the user.")') nx, ntheta0
             write (error_unit(), '("      these values lead to excess padding (",I0,".) for the FFTs.")') nx - (((ntheta0 - 1) /  2) * 3 + 1)
             write (error_unit(), '("      This may be desirable if the values chosen increase the number of parallelization sweetspots,")')
             write (error_unit(), '("      but this also makes the resolution in ny larger than is strictly necessary.")')
          end if

       else if (ntheta0 > 2*((nx-1)/3) + 1) then ! Insufficient padding
          if (proc0) then
             write (error_unit(), '("ERROR: nx (",I0,") and ntheta0 (",I0,") have been set by the user.")') nx, ntheta0
             write (error_unit(), '("       but these values lead to insufficient padding (",I0,") for the FFTS.")') nx - (((ntheta0 - 1) /  2) * 3 + 1)
             write (error_unit(), '("       The zero-padding must satisfy the two thirds rule to avoid aliasing, which requires ntheta0 <= 2*((nx-1)/3)+1.")')
             write (error_unit(), '("  Please do one of the following:")')
             write (error_unit(), '("       1. Increase nx or decrease ntheta0 to ensure ntheta0 <= 2*((nx-1)/3)+1")')
             write (error_unit(), '("       2. Set just one of nx and ntheta0, GS2 will then set the other appropriately.")')
          end if

          call mp_abort("ERROR: ntheta0 and nx both set resulting in insufficient padding. See error file for more details.", .true.)
       end if
    end if

    naky_private = naky
    ntheta0_private = ntheta0
    nx_private = nx
    ny_private = ny
  end subroutine init_kt_grids_box