1821477209

Committed 16 May 2025 02:50PM UTC coverage: 8.139% (+0.2%) from 7.92%

Build Type

push

gitlab-ci

Committed by David Dickinson

Commit Message

Merged in bugfix/use_uv_in_coverage_test_to_install_packages (pull request #1142)

Run Details

3704 of 45511 relevant lines covered (8.14%)

122643.73 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/programs/eiktest.f90

!> FIXME : Add documentation
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
  use warning_helpers, only: not_exactly_equal, exactly_equal, is_not_zero
  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 (not_exactly_equal(rmin, -1.)) write(*,*) "Option deprecated : rmin set but not used"
  if (not_exactly_equal(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 (exactly_equal(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(is_not_zero(bb)) 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 = ',e11.4,' drhodpsi = ',e11.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

1	!> FIXME : Add documentation
2	program eiktest	×
3	use mp, only: init_mp, finish_mp, proc0	×
4	use constants, only: pi, run_name_size
5	!> Inputs
6	use geometry, only: surf, eqfile, equal_arc, use_large_aspect
7	use geometry, only: bishop, irho, force_sym, efit_eq, dfit_eq, writelots
8	use geometry, only: gen_eq, ppl_eq, local_eq, idfit_eq, chs_eq
9	use geometry, only: p_prime_input, invLp_input, beta_prime_input, big, gs2d_eq
10	use geometry, only: transp_eq, eqnormfile, alpha_input, dp_mult
11	!> Methods from geometry
12	use geometry, only: eikcoefs_output_type, eikcoefs, geom
13	use geometry, only: pbarfun, pbarofrho, rpofrho, diameter
14	use file_utils, only: init_file_utils, run_name
15	use gs2_io_grids, only: nc_write_grid_sizes, nc_write_grid_scalars, nc_write_grids
16	use gs2_io_grids, only: nc_grid_file_open, nc_grid_file_close
17	use warning_helpers, only: not_exactly_equal, exactly_equal, is_not_zero
18	implicit none
19
20	integer :: i, nth, ntgrid, ntheta, shotnum, geoType
21	real :: q, qq, qplus, qmnus
22	real :: bb, profile_fac
23	real :: rk, rkm, rkp, rkpri
24	real :: rt, rtm, rtp, rtpri
25	real :: rd, rdm, rdp, rdpri
26	real :: rhocm, rhocp
27	real :: bbar, rmin, rmax, pbar
28	real :: akappa, akappri, tri, tripri
29	real :: rhoc, rmaj, r_geo, qinp, s_hat_input, shift
30
31	real :: diffscheme, beta_p1, beta_p2, beta_prime_times, beta_prime_over
32	real :: tstar
33	integer :: nbeta, funit
34	logical :: no_end_point
35	logical :: fast, dipole, list, Xanthopoulos, in_nt
36	character(len=:), allocatable :: in_stem, out_stem	×
37
38	character(len=run_name_size) :: ncout
39	type(eikcoefs_output_type) :: eikcoefs_results	×
40	real, dimension(:), allocatable :: theta, bmag, aplot, aprime, bpol, cvdrift, grho	×
41	real, dimension(:), allocatable :: cvdrift0, gbdrift, gbdrift0, gds2, gds21, gds22	×
42	real, dimension(:), allocatable :: gradpar, zplot, zprime, rprime, rplot	×
43	real :: qsf, surfarea, kxfac, dvdrhon, shat, drhodpsi, dbetadrho, delrho
44	integer :: fort11_unit, isym, iflux, itor, nperiod
45	namelist/stuff/ntheta,nperiod,geoType,rmaj,akappri,akappa,shift,equal_arc, &
46	rhoc,rmin,rmax,itor,qinp,iflux,delrho,tri,bishop, &
47	irho,isym,tripri,efit_eq,dfit_eq,writelots,R_geo, &
48	gen_eq, ppl_eq, eqfile, local_eq, idfit_eq,&
49	chs_eq, force_sym, s_hat_input,p_prime_input,invLp_input,beta_prime_input, &
50	diffscheme,nbeta,beta_p1,beta_p2,alpha_input,big, &
51	beta_prime_times, beta_prime_over, fast, profile_fac, &
52	tstar, shotnum, gs2d_eq, transp_eq, Xanthopoulos, dipole, &
53	eqnormfile, no_end_point, in_nt
54
55	call init_mp	×
56
57	! set defaults
58	gen_eq = .false.	×
59	ppl_eq = .false.	×
60	chs_eq = .false.	×
61	idfit_eq = .false.	×
62	dfit_eq = .false.	×
63	efit_eq = .false.	×
64	gs2d_eq = .false.	×
65	transp_eq = .false.	×
66	Xanthopoulos = .false.	×
67	dipole = .false.	×
68	shotnum=1001109020	×
69	tstar = 1.0	×
70	ntheta=32; nperiod=1	×
71	equal_arc = .false. ; bishop = 0 ; fast = .true.	×
72	big = 8	×
73	eqfile='dskeq.cdf'	×
74	rhoc=0.5	×
75	rmin=-1.0 ; rmax=-1.0	×
76	beta_prime_times = -1.	×
77	profile_fac = 0.5	×
78	geoType = 0	×
79	rmaj = 3 ; akappa = 1 ; tri = 0 ; qinp = 2	×
80	R_geo = 3 ; akappri = 0 ; tripri=0 ; shat = 1	×
81	shift = 0	×
82	s_hat_input = shat	×
83	beta_prime_input = -1	×
84	p_prime_input = -2	×
85	invLp_input = 5	×
86	irho = 2 ; isym = -1 ; itor = -1 ; iflux=-1	×
87	delrho = 0.01	×
88	no_end_point = .false. ; in_nt = .false.	×
89
90	if (command_argument_count() == 0) then	×
91	allocate(character(len=3)::in_stem)	×
92	in_stem = 'eik'	×
93	allocate(character(len=0)::out_stem)	×
94	out_stem = ''	×
95	else
96	call init_file_utils(list)	×
97	allocate(character(len=len_trim(run_name))::in_stem)	×
98	allocate(character(len=len_trim(run_name)+1)::out_stem)	×
99	in_stem = trim(run_name)	×
100	out_stem = trim(in_stem)//'.'	×
101	endif
102
103	! read in data
104	open(newunit=funit,file=in_stem//'.in')	×
105	read(funit,nml=stuff)	×
106	close(funit)	×
107
108	if (local_eq .and. (geoType /= 0)) print*,'Warning : geoType /= 0 not supported'	×
109	if (Xanthopoulos) write(,) "Option deprecated : Xanthopoulos set but unused"	×
110	if (in_nt) write(,) "Option deprecated : in_nt set but unused"	×
111	if (.not. fast) write(,) "Option deprecated : fast set but doesn't do anything"	×
112	if (not_exactly_equal(rmin, -1.)) write(,) "Option deprecated : rmin set but not used"	×
113	if (not_exactly_equal(rmax, -1.)) write(,) "Option deprecated : rmax set but not used"	×
114	if (isym /= -1) write(,) "Option deprecated : isym set but not used (see force_sym)"	×
115	if (iflux /= -1) write(,) "Option deprecated : iflux set but not used"	×
116	if (itor /= -1) write(,) "Option deprecated : itor set but not used (see use_large_aspect)"	×
117
118	surf%geoType = geoType	×
119	surf%mMode = 2 ; surf%nMode = 3	×
120	surf%rmaj = rmaj ; surf%R_geo = R_geo ; surf%r = rhoc ; surf%dr = delrho	×
121	surf%sHorz = shift ; surf%delm = akappa	×
122	surf%deln = tri ; surf%delmp = akappri	×
123	surf%delnp = tripri ; surf%thm = 0	×
124	surf%thn = 0 ; surf%q = qinp ; surf%shat = s_hat_input	×
125
126	if(bishop >= 7) then	×
127	if (exactly_equal(beta_prime_times, -1.)) then	×
128	write(,) 'beta_prime_times should be set for bishop = 7 or 8'	×
129	beta_prime_times = 1.	×
130	endif
131	dp_mult = beta_prime_times	×
132	endif
133	fort11_unit = 12	×
134
135	! Note that if .not. local_eq then always choose use_large_aspect = F
136	if ((.not. local_eq) .and. use_large_aspect) then	×
137	use_large_aspect = .false.	×
138	write(,) 'Forcing use_large_aspect to false'	×
139	end if
140
141	! Catch if no q profile
142	if (local_eq .and. irho /= 2) then	×
143	irho = 2	×
144	write(,) 'Forcing irho=2'	×
145	end if
146
147	if(local_eq) then	×
148	if (gen_eq) write(,) 'Forcing gen_eq to be false'	×
149	if (efit_eq) write(,) 'Forcing efit_eq to be false'	×
150	if (dfit_eq) write(,) 'Forcing dfit_eq to be false'	×
151	if (gs2d_eq) write(,) 'Forcing gs2d_eq to be false'	×
152	if (idfit_eq) write(,) 'Forcing idfit_eq to be false'	×
153	if (ppl_eq) write(,) 'Forcing ppl_eq to be false'	×
154	if (chs_eq) write(,) 'Forcing chs_eq to be false'	×
155	if (transp_eq) write(,) 'Forcing transp_eq to be false'	×
156	if (any([gen_eq, ppl_eq, transp_eq, chs_eq, efit_eq, dfit_eq, gs2d_eq, idfit_eq])) &	×
157	write(,) 'because local_eq'	×
158	gen_eq = .false. ; ppl_eq = .false. ; transp_eq = .false. ; chs_eq = .false.	×
159	efit_eq = .false. ; dfit_eq = .false. ; gs2d_eq = .false. ; idfit_eq = .false.	×
160	endif
161
162	open(newunit=funit,file=out_stem//'eik.out',status='unknown')	×
163
164	call eikcoefs(ntheta, nperiod, eikcoefs_results)	×
165
166	! Unpack results for ease of use
167	allocate(theta, source = eikcoefs_results%theta)	×
168	allocate(bmag, source = eikcoefs_results%bmag)	×
169	allocate(gradpar, source = eikcoefs_results%gradpar)	×
170	allocate(grho, source = eikcoefs_results%grho)	×
171	allocate(cvdrift, source = eikcoefs_results%cvdrift)	×
172	allocate(gbdrift, source = eikcoefs_results%gbdrift)	×
173	allocate(cvdrift0, source = eikcoefs_results%cvdrift0)	×
174	allocate(gbdrift0, source = eikcoefs_results%gbdrift0)	×
175	allocate(gds2, source = eikcoefs_results%gds2)	×
176	allocate(gds21, source = eikcoefs_results%gds21)	×
177	allocate(gds22, source = eikcoefs_results%gds22)	×
178	allocate(rplot, source = eikcoefs_results%rplot)	×
179	allocate(zplot, source = eikcoefs_results%zplot)	×
180	allocate(aplot, source = eikcoefs_results%aplot)	×
181	allocate(rprime, source = eikcoefs_results%rprime)	×
182	allocate(zprime, source = eikcoefs_results%zprime)	×
183	allocate(aprime, source = eikcoefs_results%aprime)	×
184	allocate(bpol, source = eikcoefs_results%bpol)	×
185	kxfac = eikcoefs_results%kxfac	×
186	qsf = eikcoefs_results%qsf	×
187	surfarea = eikcoefs_results%surfarea	×
188	dvdrhon = eikcoefs_results%dvdrhon	×
189	shat = eikcoefs_results%shat	×
190	drhodpsi = eikcoefs_results%drhodpsi	×
191	dbetadrho = eikcoefs_results%dbetadrho	×
192
193	! note that ntheta may have changed
194	ntgrid = (size(theta)-1)/2 ; nth = ntheta / 2	×
195
196	if (writelots) then	×
197	bbar = 0.	×
198	do i=-ntgrid,ntgrid-1	×
199	bbar = bbar + bmag(i)*(theta(i+1)-theta(i))	×
200	end do
201	bbar = bbar/(2.*pi)	×
202
203	write(,) 'B_bar = ',bbar	×
204	end if
205
206	if (dipole) then	×
207	shat = max(1.e-7, shat)	×
208	end if
209
210	write(funit,*) ' ntgrid nperiod ntheta, drhodpsi, rmaj, shat, kxfac, q'	×
211	!CMR 14/01/05 change format to prevent line breaks with intel compiler and allow runngridgen to read the output files
212	write(funit,fmt='(3i8,1p5e12.4)') ntgrid, nperiod, ntheta, drhodpsi, rmaj, shat, kxfac, qsf	×
213	!CMR
214
215	write(funit,*) ' gbdrift gradpar grho tgrid'	×
216	do i= -ntgrid,ntgrid	×
217	write(funit,1000) gbdrift(i), gradpar(i),grho(i),theta(i)	×
218	enddo
219
220	write(funit,*) ' cvdrift gds2 bmag tgrid'	×
221	do i= -ntgrid,ntgrid	×
222	write(funit,1000) cvdrift(i), gds2(i),bmag(i),theta(i)	×
223	enddo
224
225	write(funit,*) ' gds21 gds22 tgrid'	×
226	do i= -ntgrid,ntgrid	×
227	write(funit,1000) gds21(i), gds22(i),theta(i)	×
228	enddo
229
230	write(funit,*) ' cvdrift0 gbdrift0 tgrid'	×
231	do i= -ntgrid,ntgrid	×
232	write(funit,1000) cvdrift0(i),gbdrift0(i),theta(i)	×
233	enddo
234
235	call geom%hahm_burrell(profile_fac)	×
236	dipole = dfit_eq .or. idfit_eq	×
237
238	if (.not. dipole) then	×
239	pbar = pbarfun(geom%rfun(rpofrho(rhoc),0.),0.)	×
240	bb = geom%betafun(pbar)	×
241	write(funit,) '-8pi/B_0*2 dp/drho= '	×
242	if(is_not_zero(bb)) then	×
243	write(funit,1000) -dbetadrho	×
244	write(funit,1000) -dbetadrho/bb, -0.5*dbetadrho/bb	×
245	write(funit,1000) bb, bb/2.	×
246	endif
247	end if
248	close(funit)	×
249
250	open(newunit=funit,file=out_stem//'eik2.out',status='unknown')	×
251	write(funit,*) 'dV/drhon= ',dvdrhon	×
252	pbar = pbarofrho(rhoc)	×
253	if (.not. dipole) then	×
254	write (funit,*)	×
255	write (funit,*) 'The following q and shat values'	×
256	write (funit,*) 'were calculated directly'	×
257	write (funit,*) 'from the equilibrium.'	×
258	write (funit,*) 'They should match the ones above' !EGH	×
259	q = geom%qfun(0.)	×
260	write(funit,*) 'q_0= ',q	×
261	q = geom%qfun(pbar)	×
262	write(funit,*) 'rho= ',rhoc,' q(rho)= ',q	×
263	q = geom%qfun(0.95)	×
264	write(funit,*) 'q_95= ',q	×
265
266	call geofax(rhoc, theta, rt, rk, rd, nth, ntgrid)	×
267	if (.not. local_eq) then	×
268	rhocp=rhoc+delrho	×
269	rhocm=rhoc-delrho	×
270	call geofax(rhocp, theta, rtp, rkp, rdp, nth, ntgrid)	×
271	call geofax(rhocm, theta, rtm, rkm, rdm, nth, ntgrid)	×
272	qplus = geom%qfun(pbarofrho(rhocp))	×
273	qmnus = geom%qfun(pbarofrho(rhocm))	×
274	rkpri=0.5*(rkp-rkm)/delrho	×
275	rtpri=0.5*(rtp-rtm)/delrho	×
276	rdpri=0.5*(rdp-rdm)/delrho	×
277	if (dipole) then	×
278	shat = 0.	×
279	else
280	shat = 0.5 * rhoc * (qplus - qmnus) / geom%qfun(pbar) / delrho	×
281	endif
282	else
283	rkpri=akappri	×
284	rtpri=tripri	×
285	rdpri=shift	×
286	endif
287
288	qq = geom%qfun(pbar)	×
289	write(funit,*) 'shat= ',shat	×
290	write(funit,*) 'qinp= ',qq	×
291	write(funit,*) 'shift= ',rdpri	×
292	write(funit,*) 'akappa= ',rk	×
293	write(funit,*) 'akappri= ',rkpri	×
294	write(funit,*) 'tri= ',rt	×
295	write(funit,*) 'tripri= ',rtpri	×
296	write(funit,) 'surface area= ',surfarea,' avgrmid*2'	×
297	write(funit, *) 'diameter/a = ', diameter(rpofrho(rhoc))	×
298	write(funit, *) 'diameter(rhoc=1.0)/a/2.0 = ', diameter(rpofrho(1.0))/2.0, ' ... should be 1'	×
299	write(funit, *) 'r/a = ', diameter(rpofrho(rhoc))/diameter(rpofrho(1.0))	×
300	end if
301
302	close(funit)	×
303
304	open (newunit=funit,file='eik5.out',status='unknown')	×
305	write (funit,'(''#theta1 gradpar2 bmag3 grho4 '', &
306	& ''gbdrift5 gbdrift06 cvdrift7 cvdrift08 '', &
307	& ''gds29 gds2110 gds2211'')')	×
308	do i = -ntgrid, ntgrid	×
309	write (funit,1000) theta(i),gradpar(i),bmag(i),grho(i), &	×
310	gbdrift(i),gbdrift0(i),cvdrift(i),cvdrift0(i), &	×
311	gds2(i),gds21(i),gds22(i)	×
312	end do
313	close (funit)	×
314
315	open (newunit=funit,file='eik6.out',status='unknown')	×
316	write (funit,fmt="('# shape: torus')")	×
317	write (funit,fmt="('# q = ',e11.4,' drhodpsi = ',e11.4)") qsf, drhodpsi	×
318	write (funit,fmt="('# theta1 R2 Z3 alpha4 ', &
319	& ' Rprime5 Zprime6 alpha_prime7 bpol8')")	×
320	do i = -ntgrid, ntgrid	×
321	write (funit,1000) theta(i),Rplot(i),Zplot(i),aplot(i), &	×
322	Rprime(i),Zprime(i),aprime(i),Bpol(i)	×
323	end do
324	close (funit)	×
325	1000 format(20(1x,1pg18.11))
326
327	if (proc0) then	×
328	ncout = trim(out_stem)//"out.nc"	×
329	call nc_grid_file_open(ncout, "w")	×
330	call nc_write_grid_sizes(ntheta, ntgrid, nperiod)	×
331	call nc_write_grid_scalars(shat, drhodpsi, kxfac, qsf, rmaj, geom%B_T, geom%aminor, surfarea/dvdrhon, surfarea)	×
332	call nc_write_grids(ntgrid, bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, &
333	gds2, gds21, gds22, grho, theta, &
334	Rplot=Rplot, Rprime=Rprime, Zplot=Zplot, Zprime=Zprime, aplot=aplot, aprime=aprime, &
335	Bpol=Bpol, no_end_point_in=no_end_point)	×
336	call nc_grid_file_close()	×
337	end if
338
339	call finish_mp	×
340
341	contains
342
343	!> FIXME : Add documentation
344	subroutine geofax(rho, theta, t, e, d, nth, ntgrid)	×
345	use geometry, only: rpofrho, geom
346	implicit none
347	integer, intent(in) :: nth, ntgrid
348	real, intent (in) :: rho
349	real, dimension(-ntgrid:ntgrid), intent(in) :: theta
350	real, intent (out) :: t, e, d
351	real, dimension(-nth:nth) :: rgrid	×
352	real :: a, b, aa, bb, rp, hmaxu, hmaxl, rmaxu, rmaxl, rmid, Rlo, Rhi, Rinnr, Routr
353	integer :: i
354	logical, save :: first=.true.
355	if (first) then	×
356	write(*,'("WARNING: GEOFAX USES UNINITIALISED VALUES -- RESULT MAY NOT BE SENSIBLE")')	×
357	first = .false.	×
358	end if
359	rp = rpofrho(rho)	×
360	do i = -nth, nth	×
361	rgrid(i) = geom%rfun(rp, theta(i))	×
362	end do
363
364	! do not treat up-down asymmetry completely for now
365	hmaxu = -1.e6 ; hmaxl = 1.e6	×
366	rmaxu = 0. ; rmaxl = 0.	×
367
368	do i = 0, nth	×
369	a = geom%Rpos(rgrid(i), theta(i)) ; b = geom%Zpos(rgrid(i), theta(i))	×
370
371	aa = geom%Rpos(rgrid(-i), theta(-i)) ; bb = geom%Zpos(rgrid(-i), theta(-i))	×
372
373	if (hmaxu < b ) then	×
374	hmaxu = b ; rmaxu = a	×
375	end if
376
377	if (hmaxl > bb ) then	×
378	hmaxl = bb ; rmaxl = aa	×
379	end if
380
381	Rlo = min(a, aa) ; Rhi = max(a, aa)	×
382	!<DD>WARNING : Rinnr has not been defined on first loop through
383	if (Rinnr > Rlo) Rinnr = Rlo	×
384	!<DD>WARNING : Routr has not been defined on first loop through
385	if (Routr < Rhi) Routr = Rhi	×
386	end do
387	a = (geom%Rpos(rgrid(0), theta(0)) - geom%Rpos(rgrid(nth), theta(nth)))	×
388	e = (hmaxu - hmaxl) / a	×
389	! Ro=(Routr+Rinnr)/2. should = rcenter(rp)
390	rmid = (Routr - Rinnr) / 2. ! midplane minor radius of this surface	×
391	! Using only the upper triangularity:
392	! wasn't properly normalized before 29 Feb 08, and don't use asin:
393	! t = asin(rcenter(rp)-rmax)
394	t = (geom%rcenter(rp) - rmaxu) / rmid	×
395	d = geom%rcenter(rp) - geom%rcenter(geom%psi_a)	×
396	end subroutine geofax	×
397	end program eiktest

gyrokinetics / gs2 / 1821477209

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous