check_theta_grid_eik Subroutine

public subroutine check_theta_grid_eik(report_unit, dbetadrho)

FIXME : Add documentation

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: report_unit
real, intent(in) :: dbetadrho

Contents

Source Code


Source Code

  subroutine check_theta_grid_eik(report_unit, dbetadrho)
    use theta_grid_params, only: eps
    use geometry, only: surf, dp_mult, alpha_input, beta_prime_input, invLp_input
    use geometry, only: bishop, irho, eqfile, eqnormfile
    use geometry, only: idfit_eq, gen_eq, efit_eq, ppl_eq, local_eq, dfit_eq, chs_eq
    implicit none
    integer, intent(in) :: report_unit
    real, intent(in) :: dbetadrho
    real :: beta_prime_new, shat, rhoc
    beta_prime_new = eikcoefs_results%dbetadrho
    shat = eikcoefs_results%shat
    rhoc = eikcoefs_results%rhoc
    call checklogic_theta_grid_eik(report_unit)
    write (report_unit, *)
    if (local_eq) then
       write (report_unit, fmt="('The following local equilibrium model has been selected:')")
       select case (surf%geoType)
          case (0)
             write (report_unit, fmt="('Miller')")
          case (1)
             write (report_unit, fmt="('Global')")
          case (2)
             write (report_unit, fmt="('GeneralizedEllipticity')")
          case (3)
             write (report_unit, fmt="('FourierSeries')")
          case default
             write (report_unit, fmt="('ERROR: invalid analytic geometry specification')")
       end select
       if (surf%rmaj == 1.0) then
          write (report_unit, &
               & fmt="('Scale lengths are normalized to the major radius, R')")
       else
          write (report_unit, fmt="('The aspect ratio R/a = ',f7.4)") surf%rmaj
       end if
       if (surf%rmaj /= surf%r_geo) then
          write (report_unit, *)
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('R_geo is not equal to Rmaj.')")
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('################# WARNING #######################')")
       end if
       if (irho /= 2) then
          write (report_unit, *)
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('You have selected irho = ',i2)") irho
          write (report_unit, fmt="('For local equilibria, irho=2 is required.')")
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('################# WARNING #######################')")
       end if
       write (report_unit, *)
       write (report_unit, fmt="('The safety factor q =      ',f7.4)") surf%q
       eps = rhoc / surf%R_geo
       write (report_unit, fmt="('and epsilon == r/R =       ',f7.4)") eps
       write (report_unit, *)
       if (eps > epsilon(0.0)) then
          write (report_unit, fmt="('Trapped particles are included.')")
       else
          write (report_unit, fmt="('Trapped particles are neglected.')")
       end if
       write (report_unit, *)
       write (report_unit, fmt="('B_poloidal is determined by:')")
       if (surf%geoType == 1) then
          write (report_unit, *)
          write (report_unit, fmt="('    minor radius of shaped surface: aSurf = ',f7.4)") surf%aSurf
       elseif (surf%geoType == 2 .or. surf%geoType == 3) then
          write (report_unit, *)
          write (report_unit, fmt="('    first shaping mode number: mMode =  ',i2)") surf%mMode
          write (report_unit, fmt="('    second shaping mode number: nMode = ',i2)") surf%nMode
       end if
       write (report_unit, *)
       write (report_unit, fmt="('    deltam =       ',f7.4)") surf%delm
       write (report_unit, fmt="('    deltan =        ',f7.4)") surf%deln
       if (surf%geoType /= 1) then
          write (report_unit, *)
          write (report_unit, fmt="('  & gradient: d deltam /d rho =   ',f7.4)") surf%delmp
          write (report_unit, fmt="('  & gradient: d deltan /d rho = ',f7.4)") surf%delnp
       end if
       write (report_unit, *)
       write (report_unit, fmt="('    thetam =       ',f7.4)") surf%thm
       write (report_unit, fmt="('    thetan =        ',f7.4)") surf%thn
       write (report_unit, *)
       write (report_unit, fmt="('    shift =        ',f7.4)") surf%sHorz
       write (report_unit, fmt="('    shiftVert =    ',f7.4)") surf%sVert

       write (report_unit, *)
       write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") shat
       write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
       if (abs(shat) <= 1.e-5) then
          write (report_unit, fmt="('This is effectively zero; periodic boundary conditions are assumed.')")
       end if
       select case (bishop)
       case (3)
          write (report_unit, fmt="('The normalized inverse pressure gradient scale length = ',f8.4)") invLp_input
       case (4)
          write (report_unit, fmt="('The beta gradient d beta / d rho = ',f8.4)") beta_prime_input
          if (beta_prime_input > epsilon(0.0)) then
             write (report_unit, *)
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('beta_prime > 0.')")
             write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, *)
          end if
          if (abs(beta_prime_input - dbetadrho) > 1.e-2*max(abs(beta_prime_input),abs(dbetadrho))) then
             write (report_unit, *)
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('beta_prime_input is not consistent with beta and Lp.')")
             write (report_unit, fmt="('beta_prime_input=',f6.4,' dbeta_drho (from species)=',f6.4)") beta_prime_input, dbetadrho
             write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, *)
          end if
       case (5)
          write (report_unit, fmt="('The alpha parameter (R beta_prime q**2) = ',f8.4)") alpha_input
          !              write (*,*) alpha_input, dbetadrho, qinp, Rmaj
          if (abs(alpha_input + dbetadrho*surf%q**2*surf%rmaj) > 1.e-2) then
             write (report_unit, *)
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('alpha is not consistent with beta, q, and Lp.')")
             write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, *)
          end if
       case default
          write (report_unit, *)
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('You have selected bishop = ',i2)") bishop
          write (report_unit, fmt="('For local equilibria, bishop = 4 is recommended.')")
          if (bishop == 1) then
             write (report_unit, fmt="('For d beta / d rho = 0, bishop = 1 is ok.')")
             write (report_unit, fmt="('Otherwise, ')")
          end if
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *)
       end select
    end if
    if (.not. local_eq) then
       if (gen_eq) then
          write (report_unit, *)
          write (report_unit, fmt="('Equilibrium information obtained from NetCDF file:')")
          write (report_unit, fmt="(a)") trim(eqfile)
       end if
       if (ppl_eq) then
          write (report_unit, *)
          write (report_unit, fmt="('Equilibrium information obtained from NetCDF file:')")
          write (report_unit, fmt="(a)") trim(eqfile)
       end if
       if (dfit_eq) then
          write (report_unit, *)
          write (report_unit, fmt="('Dipole equilibrium information obtained from file:')")
          write (report_unit, fmt="(a)") trim(eqfile)
          write (report_unit, fmt="(a)") trim(eqnormfile)
       end if
       if (idfit_eq) then
          write (report_unit, *)
          write (report_unit, fmt="('Dipole equilibrium information obtained from file:')")
          write (report_unit, fmt="(a)") trim(eqfile)
       end if
       if (efit_eq) then
          write (report_unit, *)
          write (report_unit, fmt="('Equilibrium information obtained from eqdsk:')")
          write (report_unit, fmt="(a)") trim(eqfile)
       end if
       if (chs_eq) then
          write (report_unit, *)
          write (report_unit, fmt="('Equilibrium information obtained from chease:')")
          write (report_unit, fmt="(a)") trim(eqfile)
       end if
       select case (bishop)
       case (1)
          write (report_unit, *)
          write (report_unit, fmt="('You have set bishop=1, so dp/drho and s_hat will be found from the equilibrium file.')")
          write (report_unit, *)
       case (3)
          write (report_unit, *)
          write (report_unit, fmt="('You have set bishop=3.')")
          write (report_unit, *)
          write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") surf%shat
          write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
          if (abs(shat) <= 1.e-5) then
             write (report_unit, fmt="('This is effectively zero; periodic boundary conditions are assumed.')")
          end if
          write (report_unit, fmt="('The normalized inverse pressure gradient scale length = ',f8.4)") invLp_input
       case (4)
          write (report_unit, *)
          write (report_unit, fmt="('You have set bishop=4.')")
          write (report_unit, *)
          write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") surf%shat
          write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
          if (abs(shat) <= 1.e-5) then
             write (report_unit, fmt="('This is effectively zero; periodic boundary conditions are assumed.')")
          end if
          write (report_unit, fmt="('The beta gradient d beta / d rho = ',f8.4)") beta_prime_input
          if (beta_prime_input > epsilon(0.0)) then
             write (report_unit, *)
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('beta_prime > 0.')")
             write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, *)
          end if
          if (abs(beta_prime_input - dbetadrho) > 1.e-2*abs(dbetadrho)) then
             write (report_unit, *)
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('beta_prime_input is not consistent with beta and Lp.')")
             write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, *)
          end if
       case (5)
          write (report_unit, *)
          write (report_unit, fmt="('You have set bishop=5.')")
          write (report_unit, *)
          write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") surf%shat
          write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
          if (abs(shat) <= 1.e-5) then
             write (report_unit, fmt="('This is effectively zero; periodic boundary conditions are assumed.')")
          end if
          write (report_unit, fmt="('The alpha parameter (R beta_prime q**2) = ',f8.4)") alpha_input
          write (*,*) alpha_input, dbetadrho, surf%q, surf%rmaj
          if (abs(alpha_input + dbetadrho*surf%q**2*surf%rmaj) > 1.e-2) then
             write (report_unit, *)
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('alpha is not consistent with beta, q, and Lp.')")
             write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, *)
          end if
       case (6)
          write (report_unit, *)
          write (report_unit, fmt="('You have set bishop=6.')")
          write (report_unit, *)
          write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") surf%shat
          write (report_unit, fmt="('This value is set by s_hat_input in the theta_grid_eik_knobs namelist.')")
          if (abs(shat) <= 1.e-5) then
             write (report_unit, fmt="('This is effectively zero; periodic boundary conditions are assumed.')")
          end if
          write (report_unit, fmt="('The value of dp/drho will be found from the equilibrium file.')")
       case (7)
          write (report_unit, *)
          write (report_unit, fmt="('You have set bishop=7.')")
          write (report_unit, fmt="('The value of s_hat will be found from the equilibrium file.')")
          write (report_unit, fmt="('The magnetic shear s_hat = ',f7.4)") shat
          write (report_unit, fmt="('The value of dp/drho found from the equilibrium file will be multiplied by',f10.4)") dp_mult
          write (report_unit, fmt="('to give beta gradient d beta / d rho = ',f8.4)") beta_prime_new

          if (abs(beta_prime_new - dbetadrho) > 1.e-2*abs(dbetadrho)) then
             write (report_unit, *)
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, fmt="('beta_prime_new is not consistent with beta and Lp.')")
             write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
             write (report_unit, fmt="('################# WARNING #######################')")
             write (report_unit, *)
          end if

       case default

          write (report_unit, *)
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, fmt="('You have selected a value for bishop that is not recommended.')")
          write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')")
          write (report_unit, fmt="('################# WARNING #######################')")
          write (report_unit, *)

       end select
    end if
  end subroutine check_theta_grid_eik