eiktest Program

program eiktest
  use mp, only: init_mp, finish_mp, proc0
  use constants, only: pi, run_name_size
  !> Inputs
  use geometry, only: surf, eqfile, equal_arc, use_large_aspect
  use geometry, only: bishop, irho, force_sym, efit_eq, dfit_eq, writelots
  use geometry, only: gen_eq, ppl_eq, local_eq, idfit_eq, chs_eq
  use geometry, only: p_prime_input, invLp_input, beta_prime_input, big, gs2d_eq
  use geometry, only: transp_eq, eqnormfile, alpha_input, dp_mult
  !> Methods from geometry
  use geometry, only: eikcoefs_output_type, eikcoefs, geom
  use geometry, only: pbarfun, pbarofrho, rpofrho, diameter
  use file_utils, only: init_file_utils, run_name
  use gs2_io_grids, only: nc_write_grid_sizes, nc_write_grid_scalars, nc_write_grids
  use gs2_io_grids, only: nc_grid_file_open, nc_grid_file_close
  implicit none

  integer :: i, nth, ntgrid, ntheta, shotnum, geoType
  real :: q, qq, qplus, qmnus
  real :: bb, profile_fac
  real :: rk, rkm, rkp, rkpri
  real :: rt, rtm, rtp, rtpri
  real :: rd, rdm, rdp, rdpri
  real :: rhocm, rhocp
  real :: bbar, rmin, rmax, pbar
  real :: akappa, akappri, tri, tripri
  real :: rhoc, rmaj, r_geo, qinp, s_hat_input, shift

  real :: diffscheme, beta_p1, beta_p2, beta_prime_times, beta_prime_over
  real :: tstar
  integer :: nbeta, funit
  logical :: no_end_point
  logical :: fast, dipole, list, Xanthopoulos, in_nt
  character(len=:), allocatable :: in_stem, out_stem

  character(len=run_name_size) :: ncout
  type(eikcoefs_output_type) :: eikcoefs_results
  real, dimension(:), allocatable :: theta, bmag, aplot, aprime, bpol, cvdrift, grho
  real, dimension(:), allocatable :: cvdrift0, gbdrift, gbdrift0, gds2, gds21, gds22
  real, dimension(:), allocatable :: gradpar, zplot, zprime, rprime, rplot
  real :: qsf, surfarea, kxfac, dvdrhon, shat, drhodpsi, dbetadrho, delrho
  integer :: fort11_unit, isym, iflux, itor, nperiod
  namelist/stuff/ntheta,nperiod,geoType,rmaj,akappri,akappa,shift,equal_arc, &
       rhoc,rmin,rmax,itor,qinp,iflux,delrho,tri,bishop, &
       irho,isym,tripri,efit_eq,dfit_eq,writelots,R_geo, &
       gen_eq, ppl_eq, eqfile, local_eq, idfit_eq,&
       chs_eq, force_sym, s_hat_input,p_prime_input,invLp_input,beta_prime_input, &
       diffscheme,nbeta,beta_p1,beta_p2,alpha_input,big, &
       beta_prime_times, beta_prime_over, fast, profile_fac, &
       tstar, shotnum, gs2d_eq, transp_eq, Xanthopoulos, dipole, &
       eqnormfile, no_end_point, in_nt

  call init_mp

  ! set defaults
  gen_eq = .false.
  ppl_eq = .false.
  chs_eq = .false.
  idfit_eq = .false.
  dfit_eq = .false.
  efit_eq = .false.
  gs2d_eq = .false.
  transp_eq = .false.
  Xanthopoulos = .false.
  dipole = .false.
  shotnum=1001109020
  tstar = 1.0
  ntheta=32; nperiod=1
  equal_arc = .false. ; bishop = 0 ; fast = .true.
  big = 8
  eqfile='dskeq.cdf'
  rhoc=0.5
  rmin=-1.0 ; rmax=-1.0
  beta_prime_times = -1.
  profile_fac = 0.5
  geoType = 0
  rmaj = 3   ;     akappa = 1   ;     tri = 0   ;     qinp = 2
  R_geo = 3  ;     akappri = 0  ;     tripri=0  ;     shat = 1
  shift = 0
  s_hat_input = shat
  beta_prime_input = -1
  p_prime_input = -2
  invLp_input = 5
  irho = 2 ;      isym = -1 ; itor = -1 ; iflux=-1
  delrho = 0.01
  no_end_point = .false. ; in_nt = .false.

  if (command_argument_count() == 0) then
     allocate(character(len=3)::in_stem)
     in_stem = 'eik'
     allocate(character(len=0)::out_stem)
     out_stem = ''
  else
     call init_file_utils(list)
     allocate(character(len=len_trim(run_name))::in_stem)
     allocate(character(len=len_trim(run_name)+1)::out_stem)
     in_stem = trim(run_name)
     out_stem = trim(in_stem)//'.'
  endif

  !     read in data
  open(newunit=funit,file=in_stem//'.in')
  read(funit,nml=stuff)
  close(funit)

  if (local_eq .and. (geoType /= 0)) print*,'Warning : geoType /= 0 not supported'
  if (Xanthopoulos) write(*,*) "Option deprecated : Xanthopoulos set but unused"
  if (in_nt) write(*,*) "Option deprecated : in_nt set but unused"
  if (.not. fast) write(*,*) "Option deprecated : fast set but doesn't do anything"
  if (rmin /= -1) write(*,*) "Option deprecated : rmin set but not used"
  if (rmax /= -1) write(*,*) "Option deprecated : rmax set but not used"
  if (isym /= -1) write(*,*) "Option deprecated : isym set but not used (see force_sym)"
  if (iflux /= -1) write(*,*) "Option deprecated : iflux set but not used"
  if (itor /= -1) write(*,*) "Option deprecated : itor set but not used (see use_large_aspect)"

  surf%geoType = geoType
  surf%mMode = 2 ; surf%nMode = 3
  surf%rmaj = rmaj ; surf%R_geo = R_geo ; surf%r = rhoc ; surf%dr = delrho
  surf%sHorz = shift ; surf%delm = akappa
  surf%deln = tri ; surf%delmp = akappri
  surf%delnp = tripri ; surf%thm = 0
  surf%thn = 0 ; surf%q = qinp ; surf%shat = s_hat_input

  if(bishop >= 7) then
     if(beta_prime_times == -1) then
        write(*,*) 'beta_prime_times should be set for bishop = 7 or 8'
        beta_prime_times = 1.
     endif
     dp_mult = beta_prime_times
  endif
  fort11_unit = 12

  ! Note that if .not. local_eq then always choose use_large_aspect = F
  if ((.not. local_eq) .and. use_large_aspect) then
     use_large_aspect = .false.
     write(*,*) 'Forcing use_large_aspect to false'
  end if

  ! Catch if no q profile
  if (local_eq .and. irho /= 2) then
     irho = 2
     write(*,*) 'Forcing irho=2'
  end if

  if(local_eq) then
     if (gen_eq) write(*,*) 'Forcing gen_eq to be false'
     if (efit_eq) write(*,*) 'Forcing efit_eq to be false'
     if (dfit_eq) write(*,*) 'Forcing dfit_eq to be false'
     if (gs2d_eq) write(*,*) 'Forcing gs2d_eq to be false'
     if (idfit_eq) write(*,*) 'Forcing idfit_eq to be false'
     if (ppl_eq) write(*,*) 'Forcing ppl_eq to be false'
     if (chs_eq) write(*,*) 'Forcing chs_eq to be false'
     if (transp_eq) write(*,*) 'Forcing transp_eq to be false'
     if (any([gen_eq, ppl_eq, transp_eq, chs_eq, efit_eq, dfit_eq, gs2d_eq, idfit_eq])) &
          write(*,*) 'because local_eq'
     gen_eq = .false. ; ppl_eq = .false. ; transp_eq = .false. ; chs_eq = .false.
     efit_eq = .false. ; dfit_eq = .false. ; gs2d_eq = .false. ; idfit_eq = .false.
  endif

  open(newunit=funit,file=out_stem//'eik.out',status='unknown')

  call eikcoefs(ntheta, nperiod, eikcoefs_results)

  ! Unpack results for ease of use
  allocate(theta, source = eikcoefs_results%theta)
  allocate(bmag, source = eikcoefs_results%bmag)
  allocate(gradpar, source = eikcoefs_results%gradpar)
  allocate(grho, source = eikcoefs_results%grho)
  allocate(cvdrift, source = eikcoefs_results%cvdrift)
  allocate(gbdrift, source = eikcoefs_results%gbdrift)
  allocate(cvdrift0, source = eikcoefs_results%cvdrift0)
  allocate(gbdrift0, source = eikcoefs_results%gbdrift0)
  allocate(gds2, source = eikcoefs_results%gds2)
  allocate(gds21, source = eikcoefs_results%gds21)
  allocate(gds22, source = eikcoefs_results%gds22)
  allocate(rplot, source = eikcoefs_results%rplot)
  allocate(zplot, source = eikcoefs_results%zplot)
  allocate(aplot, source = eikcoefs_results%aplot)
  allocate(rprime, source = eikcoefs_results%rprime)
  allocate(zprime, source = eikcoefs_results%zprime)
  allocate(aprime, source = eikcoefs_results%aprime)
  allocate(bpol, source = eikcoefs_results%bpol)
  kxfac = eikcoefs_results%kxfac
  qsf = eikcoefs_results%qsf
  surfarea = eikcoefs_results%surfarea
  dvdrhon = eikcoefs_results%dvdrhon
  shat = eikcoefs_results%shat
  drhodpsi = eikcoefs_results%drhodpsi
  dbetadrho = eikcoefs_results%dbetadrho

  ! note that ntheta may have changed
  ntgrid = (size(theta)-1)/2 ; nth = ntheta / 2

  if (writelots) then
     bbar = 0.
     do i=-ntgrid,ntgrid-1
        bbar = bbar + bmag(i)*(theta(i+1)-theta(i))
     end do
     bbar = bbar/(2.*pi)

     write(*,*) 'B_bar = ',bbar
  end if

  if (dipole) then
     shat = max(1.e-7, shat)
  end if

  write(funit,*) ' ntgrid  nperiod  ntheta, drhodpsi, rmaj, shat, kxfac, q'
  !CMR 14/01/05 change format to prevent line breaks with intel compiler and allow runngridgen to read the output files
  write(funit,fmt='(3i8,1p5e12.4)') ntgrid, nperiod, ntheta, drhodpsi, rmaj, shat, kxfac, qsf
  !CMR

  write(funit,*) '    gbdrift     gradpar       grho    tgrid'
  do i= -ntgrid,ntgrid
     write(funit,1000) gbdrift(i), gradpar(i),grho(i),theta(i)
  enddo

  write(funit,*) '    cvdrift            gds2            bmag   tgrid'
  do i= -ntgrid,ntgrid
     write(funit,1000) cvdrift(i), gds2(i),bmag(i),theta(i)
  enddo

  write(funit,*) '    gds21             gds22  tgrid'
  do i= -ntgrid,ntgrid
     write(funit,1000) gds21(i), gds22(i),theta(i)
  enddo

  write(funit,*) '    cvdrift0          gbdrift0   tgrid'
  do i= -ntgrid,ntgrid
     write(funit,1000) cvdrift0(i),gbdrift0(i),theta(i)
  enddo

  call geom%hahm_burrell(profile_fac)
  dipole = dfit_eq .or. idfit_eq

  if (.not. dipole) then
     pbar = pbarfun(geom%rfun(rpofrho(rhoc),0.),0.)
     bb = geom%betafun(pbar)
     write(funit,*) '-8*pi/B_0**2 * dp/drho= '
     if(bb.ne.0.) then
        write(funit,1000) -dbetadrho
        write(funit,1000) -dbetadrho/bb, -0.5*dbetadrho/bb
        write(funit,1000) bb, bb/2.
     endif
  end if
  close(funit)

  open(newunit=funit,file=out_stem//'eik2.out',status='unknown')
  write(funit,*) 'dV/drhon= ',dvdrhon
  pbar = pbarofrho(rhoc)
  if (.not. dipole) then
     write (funit,*)
     write (funit,*) 'The following q and shat values'
     write (funit,*) 'were calculated directly'
     write (funit,*) 'from the equilibrium.'
     write (funit,*) 'They should match the ones above' !EGH
     q = geom%qfun(0.)
     write(funit,*) 'q_0= ',q
     q = geom%qfun(pbar)
     write(funit,*) 'rho= ',rhoc,' q(rho)= ',q
     q = geom%qfun(0.95)
     write(funit,*) 'q_95= ',q

     call geofax(rhoc, theta, rt, rk, rd, nth, ntgrid)
     if (.not. local_eq) then
        rhocp=rhoc+delrho
        rhocm=rhoc-delrho
        call geofax(rhocp, theta, rtp, rkp, rdp, nth, ntgrid)
        call geofax(rhocm, theta, rtm, rkm, rdm, nth, ntgrid)
        qplus = geom%qfun(pbarofrho(rhocp))
        qmnus = geom%qfun(pbarofrho(rhocm))
        rkpri=0.5*(rkp-rkm)/delrho
        rtpri=0.5*(rtp-rtm)/delrho
        rdpri=0.5*(rdp-rdm)/delrho
        if (dipole) then
           shat = 0.
        else
           shat = 0.5 * rhoc * (qplus - qmnus) / geom%qfun(pbar) / delrho
        endif
     else
        rkpri=akappri
        rtpri=tripri
        rdpri=shift
     endif

     qq = geom%qfun(pbar)
     write(funit,*) 'shat= ',shat
     write(funit,*) 'qinp= ',qq
     write(funit,*) 'shift= ',rdpri
     write(funit,*) 'akappa= ',rk
     write(funit,*) 'akappri= ',rkpri
     write(funit,*) 'tri= ',rt
     write(funit,*) 'tripri= ',rtpri
     write(funit,*) 'surface area= ',surfarea,' avgrmid**2'
     write(funit, *) 'diameter/a = ', diameter(rpofrho(rhoc))
     write(funit, *) 'diameter(rhoc=1.0)/a/2.0 = ', diameter(rpofrho(1.0))/2.0, ' ...  should be 1'
     write(funit, *) 'r/a = ', diameter(rpofrho(rhoc))/diameter(rpofrho(1.0))
  end if

  close(funit)

  open (newunit=funit,file='eik5.out',status='unknown')
  write (funit,'(''#theta1 gradpar2 bmag3 grho4 '', &
       &   ''gbdrift5 gbdrift06 cvdrift7 cvdrift08 '', &
       &    ''gds29 gds2110 gds2211'')')
  do i = -ntgrid, ntgrid
     write (funit,1000) theta(i),gradpar(i),bmag(i),grho(i), &
          gbdrift(i),gbdrift0(i),cvdrift(i),cvdrift0(i), &
          gds2(i),gds21(i),gds22(i)
  end do
  close (funit)

  open (newunit=funit,file='eik6.out',status='unknown')
  write (funit,fmt="('# shape: torus')")
  write (funit,fmt="('# q = ',e10.4,' drhodpsi = ',e10.4)") qsf, drhodpsi
  write (funit,fmt="('# theta1             R2                  Z3               alpha4      ', &
            &   '       Rprime5              Zprime6           alpha_prime7         bpol8')")
  do i = -ntgrid, ntgrid
     write (funit,1000) theta(i),Rplot(i),Zplot(i),aplot(i), &
          Rprime(i),Zprime(i),aprime(i),Bpol(i)
  end do
  close (funit)
1000 format(20(1x,1pg18.11))

  if (proc0) then
     ncout = trim(out_stem)//"out.nc"
     call nc_grid_file_open(ncout, "w")
     call nc_write_grid_sizes(ntheta, ntgrid, nperiod)
     call nc_write_grid_scalars(shat, drhodpsi, kxfac, qsf, rmaj, geom%B_T, geom%aminor, surfarea/dvdrhon, surfarea)
     call nc_write_grids(ntgrid, bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, &
          gds2, gds21, gds22, grho, theta, &
          Rplot=Rplot, Rprime=Rprime, Zplot=Zplot, Zprime=Zprime, aplot=aplot, aprime=aprime, &
          Bpol=Bpol, no_end_point_in=no_end_point)
     call nc_grid_file_close()
  end if

  call finish_mp

contains

  !> FIXME : Add documentation
  subroutine geofax(rho, theta, t, e, d, nth, ntgrid)
    use geometry, only: rpofrho, geom
    implicit none
    integer, intent(in) :: nth, ntgrid
    real, intent (in) :: rho
    real, dimension(-ntgrid:ntgrid), intent(in) :: theta
    real, intent (out) :: t, e, d
    real, dimension(-nth:nth) :: rgrid
    real :: a, b, aa, bb, rp, hmaxu, hmaxl, rmaxu, rmaxl, rmid, Rlo, Rhi, Rinnr, Routr
    integer :: i
    logical, save :: first=.true.
    if (first) then
       write(*,'("WARNING: GEOFAX USES UNINITIALISED VALUES -- RESULT MAY NOT BE SENSIBLE")')
       first = .false.
    end if
    rp = rpofrho(rho)
    do i = -nth, nth
       rgrid(i) = geom%rfun(rp, theta(i))
    end do

    ! do not treat up-down asymmetry completely for now
    hmaxu = -1.e6 ; hmaxl = 1.e6
    rmaxu = 0. ; rmaxl = 0.

    do i = 0, nth
       a = geom%Rpos(rgrid(i), theta(i)) ; b = geom%Zpos(rgrid(i), theta(i))

       aa = geom%Rpos(rgrid(-i), theta(-i)) ; bb = geom%Zpos(rgrid(-i), theta(-i))

       if (hmaxu < b ) then
          hmaxu = b ; rmaxu = a
       end if

       if (hmaxl > bb ) then
          hmaxl = bb ; rmaxl = aa
       end if

       Rlo = min(a, aa) ; Rhi = max(a, aa)
       !<DD>WARNING : Rinnr has not been defined on first loop through
       if (Rinnr > Rlo) Rinnr = Rlo
       !<DD>WARNING : Routr has not been defined on first loop through
       if (Routr < Rhi) Routr = Rhi
    end do
    a = (geom%Rpos(rgrid(0), theta(0)) - geom%Rpos(rgrid(nth), theta(nth)))
    e = (hmaxu - hmaxl) / a
    ! Ro=(Routr+Rinnr)/2. should = rcenter(rp)
    rmid = (Routr - Rinnr) / 2. ! midplane minor radius of this surface
    ! Using only the upper triangularity:
    ! wasn't properly normalized before 29 Feb 08, and don't use asin:
    !    t = asin(rcenter(rp)-rmax)
    t = (geom%rcenter(rp) - rmaxu) / rmid
    d = geom%rcenter(rp) - geom%rcenter(geom%psi_a)
  end subroutine geofax
end program eiktest