ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
mapping_mod.f90
Go to the documentation of this file.
1 module mapping_mod
2 
3  use itm_types
4  use mod_mesh
5 
6  implicit none
7 
8  type rz_coordinate
9  real(r8) :: R, Z
10  end type rz_coordinate
11 
12  real(r8), allocatable, private :: cchi(:, :)
13 
14  contains
15 
16  subroutine mapping(cx, cy, xaxis, yaxis, a, equilibrium_out)
17 
18 !-----------------------------------------------------------------------
19 ! subroutine to calculate the metric coefficients needed for castor
20 !-----------------------------------------------------------------------
21 
22  use itm_types
23  use mod_parameter
24  use mod_dat
25  use mod_dete
26  use mod_equilibrium
27  use mod_gauss
28  use mod_map
29  use mod_nodes
30  use mod_output
31  use mod_profiles
32  use mod_suydam
33  use mod_helena_io
34  use mod_interpolation
35  use plot_data
36  use mod_ne_te_profile
37  use math_functions
38  use phys_functions
39  use helena_spline
40 
41  use euitm_schemas
42 
43  implicit none
44 
45  integer(itm_i4), parameter :: out_dat = 16, out_vec = 21
46 
47  real(r8), intent(in) :: cx, cy, xaxis, yaxis, a
48 
49  type(type_equilibrium) :: equilibrium_out
50  type(rz_coordinate) :: top_point, bottom_point
51 
52  character(len = 132) :: grid_type(4)
53 
54  type(spline_coefficients) :: coeff
55  real(r8), dimension(nr) :: pskn1, work1
56  real(r8), dimension(nr) :: p0tmp, ftmp, dptmp, dftmp
57  real(r8), dimension(nr) :: ppelite, felite, ffpelite, pppelite, fppelite
58 ! type(spline_coefficients) :: dd_spline, f_spline
59  real(r8), dimension(nr) :: nne, dnne, tte, dtte
60  real(r8), dimension(2 * nr - 1) :: cuplot, xplot, pplot, &
61  psiplot, qplot, dqplot
62  real(r8), dimension((nr - 1) * nchi) :: xchi, ychi
63  real(r8) :: abltg(3)
64 
65  integer(itm_i4) :: i, j, k
66  integer(itm_i4) :: n1, n2, n3, n4, no
67  integer(itm_i4) :: nom, nop
68  integer(itm_i4) :: ierr, jerr
69  integer(itm_i4) :: jbase
70  integer(itm_i4) :: ifail
71  integer(itm_i4) :: nedgelab
72  integer(itm_i4) :: ileft, iright
73 
74  real(r8) :: maxerr, errj, serr, cerr
75  real(r8) :: factas
76  real(r8) :: pax, paxis, bm, b0, r0
77  real(r8) :: rbscale, pscale0, rbscale0
78  real(r8) :: rminor, rmag, raxis
79  real(r8) :: sumq, sumqr, zsumq, zsumqr
80  real(r8) :: zpsi, dzpsi, ddzpsi, zchi
81  real(r8) :: psir, psirr
82  real(r8) :: s, s2, s3
83  real(r8) :: x, xr, xs, xrs, xrr, xss
84  real(r8) :: y, yr, ys, yrs, yrr, yss
85  real(r8) :: dchi, chir, chis, chirs, chirr, chiss
86  real(r8) :: jacobian, jacobian2
87  real(r8) :: er, bigr
88  real(r8) :: dum, dum1, dum2, dum3, dum4, dummy
89  real(r8) :: psi_dummy
90 ! real :: xdum, ydum
91  real(r8) :: qq, dqq, rb, drb
92  real(r8) :: psix, psiy, chix, chiy
93  real(r8) :: grps2, grpgrc, grchi2
94  real(r8) :: a0, a1, a2, a3
95  real(r8) :: omega, rold, zold, bpold, r, z, leng, bp
96  real(r8) :: dl
97  real(r8) :: qcalc, qnorm
98  real(r8) :: dndpsico, dtdpsico
99  real(r8) :: defte, tefte
100  real(r8) :: conste, consde
101  real(r8) :: ps
102  real(r8) :: psr, pss, psrs, psrr, psss
103  real(r8) :: grad_psi2
104  real(r8) :: xleft, xright
105  real(r8) :: f_dia
106 
107  real(r8), dimension(nr, nchi) :: xout, yout
108  real(r8), dimension(nr, nchi) :: cr, cz
109  real(r8), dimension(nr, nchi) :: br, bz, bphi
110  real(r8), dimension(nchi) :: vx, vy
111  real(r8), dimension(nr) :: psi_out, p_out, rbphi_out, dp_out
112  real(r8), dimension(nr, 1) :: psi_out_2d
113  real(r8) :: cpsurf_out
114 
115  real(r8) :: r_geo_flux_surface, a_flux_surface
116  real(r8) :: r_min_flux_surface, r_max_flux_surface
117  real(r8) :: delta_upper, delta_lower
118 
119 !-- variables for extrapolation onto magnetic axis
120  real(r8), dimension(nr) :: c1
121 
122  real(r8) :: dgamma_dpsi, dp_dpsi
123 
124  logical chin
125 
126  real(r8) :: psi_bar
127  real(r8), dimension(nr) :: drbphi_out, rho, r_in
128  real(r8), dimension(nr) :: r_geo, rho_vol, dq_dpsi, dv_dpsi, i_tor
129  real(r8), dimension(nr) :: stability_alpha, global_shear
130  integer(itm_i4), parameter :: n_string = 39, len_string = 16
131  character, dimension(n_string) :: string*(len_string) = &
132  (/' psi_bar ', ' radius s ', ' psi ', &
133  ' psi_tilde ', ' r_inboard ', ' R_inboard ', &
134  ' r_outboard ', ' R_outboard ', ' Rgeo ', &
135  ' rho_vol ', ' rho_vol_bar ', ' rho_vol_tilde ', &
136  ' phi ', ' rho_tor ', &
137  ' area ', ' volume ', ' dV/dpsi ', &
138  ' p ', ' dp/dpsi ', ' dp/dpsi_tilde ', &
139  ' F ', ' dF/dpsi ', ' dF/dpsi_tilde ', &
140  ' FdF/dpsi_tilde', ' q ', ' dq/dpsi ', &
141  ' <j_tor> ', ' I_tor ', ' alpha ', &
142  ' global shear ', ' gm1 ', ' gm2 ', &
143  ' gm3 ', ' gm4 ', ' gm5 ', &
144  ' gm6 ', ' gm7 ', ' gm8 ', &
145  ' gm9 '/)
146  character :: header_string*(len_string * n_string)
147 
148 !-- define grid type
149  grid_type(1) = '2'
150  grid_type(2) = 'inverse'
151  grid_type(3) = '1'
152  grid_type(4) = 'straight field line'
153 
154 ! --- initialize coordinate system
155  xout = 0._r8
156  yout = 0._r8
157 
158  g11_hel = 0._r8
159  g12_hel = 0._r8
160  g33_hel = 0._r8
161  cr = 0._r8
162  cz = 0._r8
163  br = 0._r8
164  bz = 0._r8
165  bphi = 0._r8
166 
167 !---------------------------------------------------------------------
168  maxerr = -1.e20_r8
169  factas = 2._r8
170  if (ias == 1) factas =1._r8
171 !--------------------------------------------- norm factors ----------
172 ! profiles yields the pressure normalized to eps*Bvac^2/(mu_0 * alfa^2)
173 ! and R Bphi normalized to Rvac*Bvac
174  call profiles(p0, rbphi, dp0, drbphi, a)
175 
176  do i = 1, nr
177  pskn1(i) = psikn(nr + 1 - i)
178  end do
179 
180  if (standard_output) then
181  write(out_he, *)
182  write(out_he, *) '*************************************'
183  write(out_he, *) '* pressure profile before norm. : *'
184  write(out_he, *) '*************************************'
185  write(out_he, 91) (p0(i) * eps / alfa**2, i = 1, nr)
186  write(out_he, *) '*************************************'
187  write(out_he, *)
188  write(out_he, *)
189  write(out_he, *) '*************************************'
190  write(out_he, *) '* denormalized pressure profile : *'
191  write(out_he, *) '*************************************'
192  write(out_he, 91) (p0(i) * eps * bvac**2 / (alfa**2 * 4 * pi &
193  * 1e-07_r8), i = 1, nr)
194  write(out_he, *) '*************************************'
195  write(out_he, *)
196  end if
197 
198  pax = p0(1) * eps / (alfa * alfa)
199 
200  raxis = 1._r8 + eps * xaxis
201  bm = rbphi(1) / raxis
202 
203  if (verbosity > 2) then
204  write(iu6, *) '********magnetic/vacuum axis****************************'
205  write(iu6, *) 'vacuum magnetic axis rvac =', rvac,' m'
206  write(iu6, *) 'plasma magnetic axis rmag =', raxis * rvac,' m'
207  write(iu6, *) 'shift rmag/rvac =', raxis
208  write(iu6, *) 'vacuum magnetic field on axis bvac =', bvac,'T'
209  write(iu6, *) 'plasma magnetic field on axis bmag =', bm * bvac,'T'
210  write(iu6, *) 'bmag/bvac =', bm
211  write(iu6, *) '********************************************************'
212  end if
213 
214  b0 = 1._r8 / bm
215  r0 = 1._r8 / raxis
216  radius = eps * r0
217  pscale = b0**2 * eps / alfa**2
218  rbscale = b0 * r0
219  do i = 1, nr
220  p0(i) = p0(i) * pscale
221  dp0(i) = dp0(i) * pscale
222  rbphi(i) = rbphi(i) * rbscale
223  drbphi(i) = drbphi(i) * rbscale
224  end do
225  cpsurf = radius**2 * b0 / alfa
226  cpsurf_out = twopi * eps**2 * rvac**2 * bvac / alfa
227 
228  if (verbosity > 2) then
229  write(iu6, *) '++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
230  write(iu6, *) 'consistency check for total poloidal ', &
231  'at the plasma boundary (cpsurf_out):'
232  write(iu6, *) 'poloidal flux normalized to 2*Pi*bmag*rmag**2 ', &
233  'psi_1_hat = ', cpsurf
234  write(iu6, *) 'unnormalized poloidal flux psi_1 = ', cpsurf_out
235  write(iu6, *) 'unnormalized poloidal flux on input psi_1 = ', cpsurfin
236  write(iu6,*) '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
237  end if
238 
239  raxis = 1._r8
240  paxis = p0(1)
241 
242  if (standard_output) then
243  write(out_he, 2) radius
244  write(out_he, 3) b0
245  write(out_he, 4) cpsurf
246  end if
247 
248 !-- data for vector plot
249  if ((trim(output) == 'equilibrium' .or. trim(output) == 'full') &
250  .and. vec_file) then
251  open (unit = out_vec, file = trim(adjustl(path)) // file_vec, &
252  status = 'replace', form = 'formatted', action = 'write', &
253  iostat = i_error)
254  write(out_vec, *) nr, nchi, eps
255  end if
256 
257 !-- q profile
258 ! dq/ds profile
259 ! chi values at theta nodes
260  allocate(cchi(4, nr * np))
261  cchi = 0._r8
262  cs = 0._r8
263 
264  do i = 1, nr - 1
265  sumq = 0._r8
266  sumqr = 0._r8
267  zpsi = psikn(i)
268 ! call radial_mesh(radpsi(i), zpsi, dzpsi, ddzpsi)
269  psir = -dpsikn(i) / (2._r8 * (nr - 1))
270 ! minus sign may account for the sign change in psi when normalising ???
271  psirr = ddpsikn(i) / (2._r8 * (nr - 1))**2
272  cs(nr - i + 1) = sqrt(zpsi)
273  do j = 1, np - 1
274  n1 = (i - 1) * np + j
275  n2 = n1 + 1
276  n3 = n2 + np
277  n4 = n1 + np
278  do k = 1, 4
279  s = xgauss(k)
280  call interpolation(2, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
281  -1._r8, s, x, xr, xs, xrs, xrr, xss)
282  call interpolation(2, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
283  -1._r8, s, y, yr, ys, yrs, yrr, yss)
284  call interpolation(2, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) &
285  + 1), psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1), -1._r8, s, &
286  ps, psr, pss, psrs, psrr, psss)
287  jacobian = xr * ys - xs * yr
288  er = xrr * ys + xr * yrs - xrs * yr - xs * yrr
289  bigr = 1._r8 + eps * x
290 
291  dl = sqrt(xs**2 + ys**2)
292  grad_psi2 = psr**2 * (xs**2 + ys**2) / jacobian**2
293  sumq = sumq + wgauss(k) / (bigr * sqrt(grad_psi2)) * dl
294 ! sumq = sumq - wgauss(k) * jacobian / (bigr * abs(psir))
295  sumqr = sumqr + psirr * jacobian / ((psir**2) * bigr) * wgauss(k)
296  sumqr = sumqr - er / (bigr * psir) * wgauss(k)
297  sumqr = sumqr + jacobian * eps * xr / ((bigr**2) * psir) * wgauss(k)
298  end do
299  cchi(1, (i - 1) * np + j + 1) = sumq
300  cchi(2, (i - 1) * np + j + 1) = sumqr
301  call interpolation(2, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
302  -1._r8, 1._r8, x, xr, xs, xrs, xrr, xss)
303  call interpolation(2, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
304  -1._r8, 1._r8, y, yr, ys, yrs, yrr, yss)
305  jacobian = xr * ys - xs * yr
306  er = xrr * ys + xr * yrs - xrs * yr - xs * yrr
307  bigr = 1._r8 + eps * x
308  zsumq = -jacobian / (bigr * abs(psir))
309  zsumqr = psirr * jacobian / (psir**2 * bigr)
310  zsumqr = zsumqr - er / (bigr * psir)
311  zsumqr = zsumqr + jacobian * eps * xr / (bigr**2 * psir)
312  cchi(3, (i - 1) * np + j + 1) = zsumq
313  cchi(4, (i - 1) * np + j + 1) = zsumqr
314  end do
315  cchi(1, (i - 1) * np + 1) = 0._r8
316  cchi(2, (i - 1) * np + 1) = 0._r8
317  n1 = (i - 1) * np + 1
318  n2 = n1 + 1
319  n3 = n2 + np
320  n4 = n1 + np
321  call interpolation(2, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
322  -1._r8, -1._r8, x, xr, xs, xrs, xrr, xss)
323  call interpolation(2, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
324  -1._r8, -1._r8, y, yr, ys, yrs, yrr, yss)
325  jacobian = xr * ys - xs * yr
326  er = xrr * ys + xr * yrs - xrs * yr - xs * yrr
327  bigr = 1._r8 + eps * x
328  zsumq = -jacobian / (bigr * abs(psir))
329  zsumqr = psirr * jacobian / (psir**2 * bigr)
330  zsumqr = zsumqr - er / (bigr * psir)
331  zsumqr = zsumqr + jacobian * eps * xr / (bigr**2 * psir)
332  cchi(3, (i - 1) * np + 1) = zsumq
333  cchi(4, (i - 1) * np + 1) = zsumqr
334  qs(nr - i + 1) = 0.5_r8 * factas * sumq * rbphi(nr - i + 1) / rbscale
335  dqs(nr - i + 1) = sumqr * rbphi(nr - i + 1) + sumq &
336  * drbphi(nr - i + 1) / (2._r8 * (nr - 1))
337  dqs(nr - i + 1)= 0.5_r8 * factas * dqs(nr - i + 1) / rbscale
338  end do
339  if (eps <= 0.01_r8 .and. p0(1) <= 0._r8) then
340  s2 = sqrt(psikn(nr - 1))
341  s3 = sqrt(psikn(nr - 2))
342  qs(1) = qs(2) - (qs(3) - qs(2)) / (s3 - s2) * s2
343  else
344 ! qs(1) = rbphi(1) * pi / (2._r8 * sqrt(cx * cy) * (1. + eps * xaxis) &
345 ! * rbscale)
346  qs(1) = pi / (2._r8 * sqrt(cx * cy) * (1. + eps * xaxis))
347  end if
348  dqs(1) = 0._r8
349  do i = 1, nr - 1
350  do j = 1, np
351  dum = cchi(1, i * np)
352  no = (i - 1) * np + j
353  cchi(1, no) = (1._r8 + ias) * pi * cchi(1, no) / dum
354  dum2 = cchi(2, i * np)
355  cchi(2, no) = (1._r8 + ias) * pi * cchi(2, no) / dum
356  cchi(3, no) = (1._r8 + ias) * pi * cchi(3, no) / dum
357  cchi(4, no) = (1._r8 + ias) * pi * cchi(4, no) / dum
358 !-- final values of chi
359  qq = qs(nr - i + 1)
360  dqq = dqs(nr - i + 1)
361  rb = rbphi(nr - i + 1)
362  drb = drbphi(nr - i + 1) / (2 * (nr - 1))
363  cchi(2, no)= (dqq / qq - drb / rb) * cchi(1, no) - cchi(2, no)
364  cchi(4, no)= (dqq / qq - drb / rb) * cchi(3, no) - cchi(4, no)
365  end do
366  end do
367  do i = 1, nr
368  qs(i) = qs(i) * alfa / pi
369  qplot(nr + 1 - i) = qs(i)
370  dqs(i) = dqs(i) * alfa / pi
371  dqplot(nr + 1 - i) = dqs(i)
372  end do
373  qaxis = qs(1)
374 
375  if (standard_output) then
376  write(out_he, *)
377  write(out_he, 31) qs(1)
378  write(out_he, 32) qs(nr)
379  end if
380 
381 !-- determine positions of equidistant chi's
382 ! and calculate matrix elements
383  if (ias == 0) then
384  do j = 1, nchi
385  chikn(j) = pi * (j - 1) / dble(nchi - 1)
386  chi(j) = pi * (j - 1) / dble(nchi - 1)
387  end do
388  else
389  do j = 1, nchi
390  chikn(j) = 2._r8 * pi * (j - 1) / dble(nchi)
391  chi(j) = 2._r8 * pi * (j - 1) / dble(nchi)
392  end do
393  end if
394  do i = 1, nr - 1
395  call radial_mesh(radpsi(i), zpsi, dzpsi, ddzpsi)
396  psir = -dpsikn(i) / (2._r8 * dble(nr - 1))
397 !-- first point is known
398  no = (i - 1) * nchi + 1
399  k = 1
400  s = -1._r8
401  n1 = (i - 1) * np + k
402  n2 = n1 + 1
403  n3 = n2 + np
404  n4 = n1 + np
405  call interpolation(2, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
406  -1._r8, s, xchi(no), xr, xs, xrs, xrr, xss)
407  call interpolation(2, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
408  -1._r8, s, ychi(no), yr, ys, yrs, yrr, yss)
409  call interpolation(2, cchi(1, n1), cchi(1, n2), cchi(1, n3), &
410  cchi(1, n4), -1._r8, s, dchi, chir, chis, chirs, chirr, chiss)
411 !-- vacuum data
412  if (i == 1) then
413  vx(1) = xchi(no)
414  vy(1) = ychi(no)
415  endif
416 !-----------------------------------------------------------------------
417  jacobian = xr * ys - xs * yr
418  jacobian2 = jacobian**2
419 !-- data for vector file
420  psix = psir * ys / jacobian
421  psiy = -psir * xs / jacobian
422  chix = (chir * ys - chis * yr) / jacobian
423  chiy = (-chir * xs + chis * xr)/ jacobian
424  if (vec_file) write(out_vec, 61) sqrt(zpsi), dchi, xchi(no), &
425  ychi(no), psix, psiy, chix, chiy
426 !----------------------------------------------------------------------
427  grps2 = psir**2 * (xs**2 + ys**2) / jacobian2
428 !-- multiplying by cpsurf changes normalization from psi_hat (i.e.,
429 ! psi_out / cpsurf_out) to psi_bar
430 ! (i.e., psi_out / (twopi * bmag * rmag**2)) which means the normalized
431 ! poloidal flux per radian
432  g11_hel(nr - i + 1, 1) = grps2 * (cpsurf / radius)**2
433  grpgrc = (chir * psir * (xs**2 + ys**2) - chis * psir * (xr * xs &
434  + yr * ys)) / jacobian2
435  g12_hel(nr - i + 1, 1) = grpgrc * cpsurf / radius**2
436  g33_hel(nr - i + 1, 1) = ((1._r8 + eps * xaxis) / (1._r8 + eps &
437  * xchi(no)))**2
438 
439 !-- now: g11_hel normalized to (twopi * bmag * rmag**2 / rmag)**2
440 ! g12_hel normalized to twopi * bmag * rmag**2 / rmag**2
441 ! g33_hel normalized to 1 / rmag**2
442  xout(nr - i + 1, 1) = xchi(no)
443  yout(nr - i + 1, 1) = ychi(no)
444  cr(nr - i + 1, 1) = xchi(no)
445  cz(nr - i + 1, 1) = ychi(no)
446 !-- calculate magnetic field component
447 ! br, bz normalized to rvac**2 / (bmag * rmag**2)
448 ! bphi normalized to rvac / (bmag * rmag)
449  br(nr - i + 1, 1) = -psiy / (eps * (1._r8 + eps * xchi(no))) &
450  * cpsurf
451  bz(nr - i + 1, 1) = psix / (eps * (1._r8 + eps * xchi(no))) &
452  * cpsurf
453  bphi(nr - i + 1, 1) = rbphi(nr - i + 1) / (1._r8 + eps * xchi(no))
454 !-- check jacobian
455  grchi2 = chix**2 + chiy**2
456  dum1 = rbphi(nr - i + 1)**2 / (qs(nr - i + 1)**2 &
457  / g33_hel(nr - i + 1, 1))
458  dum2 = g11_hel(nr - i + 1, 1) * grchi2 / radius**2
459  dum3 = g12_hel(nr - i + 1, 1)
460  dum4 = dum2 - dum3 * dum3
461  errj = abs(dum4 - dum1)
462  if (errj > maxerr) then
463  maxerr = errj
464  ierr = i
465  jerr = j
466  serr = sqrt(zpsi)
467  cerr = dchi
468  end if
469 !---------------------------------------------------------------------
470  jbase = 2
471  do k = 1, np - 1
472  chin = .false.
473  do j = jbase, nchi
474  zchi = chikn(j)
475  no = (i - 1) * nchi + j
476  if ((cchi(1, (i - 1) * np + k) <= zchi .and. cchi(1, &
477  (i - 1) * np + k + 1) >= zchi) .or. (j == nchi &
478  .and. k == np - 1 .and. ias == 0)) then
479  chin = .true.
480  nom = (i - 1) * np + k
481  nop = nom + 1
482  a3 = (cchi(1, nom) + cchi(3, nom) - cchi(1, nop) &
483  + cchi(3, nop)) / 4._r8
484  a2 = (-cchi(3, nom) + cchi(3, nop)) / 4._r8
485  a1 = (-3 * cchi(1, nom) - cchi(3, nom) + 3 * cchi(1, nop) &
486  - cchi(3, nop)) / 4._r8
487  a0 = (2 * cchi(1, nom) + cchi(3, nom) + 2 * cchi(1, nop) &
488  - cchi(3, nop)) / 4._r8 - zchi
489  call solve_cubic(a0, a1, a2, a3, s, s2, s3, ifail)
490  if (ias == 0 .and. j == nchi) then
491  s = 1._r8
492  ifail = 0
493  end if
494  if (ifail == 0) then
495  n1 = (i - 1) * np + k
496  n2 = n1 + 1
497  n3 = n2 + np
498  n4 = n1 + np
499  call interpolation(2, xx(1, n1), xx(1, n2), xx(1, n3), &
500  xx(1, n4), -1._r8, s, xchi(no), xr, xs, xrs, xrr, xss)
501  call interpolation(2, yy(1, n1), yy(1, n2), yy(1, n3), &
502  yy(1, n4), -1._r8, s, ychi(no), yr, ys, yrs, yrr, yss)
503  call interpolation(2, cchi(1, n1), cchi(1, n2), cchi(1, n3), &
504  cchi(1, n4), -1._r8, s, dchi, chir, chis, chirs, chirr, &
505  chiss)
506 !-- vacuum data
507  if (i == 1) then
508  vx(j) = xchi(no)
509  vy(j) = ychi(no)
510  end if
511 !-----------------------------------------------------------------------
512  jacobian = xr * ys - xs * yr
513  jacobian2 = jacobian**2
514 !-- data for vector file
515  psix = psir * ys / jacobian
516  psiy = -psir * xs / jacobian
517  chix = (chir * ys - chis * yr) / jacobian
518  chiy = (-chir * xs + chis * xr) / jacobian
519  if (vec_file) write(out_vec, 61) sqrt(zpsi), dchi, xchi(no), &
520  ychi(no), psix, psiy, chix, chiy
521 !----------------------------------------------------------------------
522  grps2 = psir**2 * (xs**2 + ys**2) / jacobian2
523  g11_hel(nr - i + 1, j) = grps2 * (cpsurf / radius)**2
524  grpgrc = (chir * psir * (xs**2 + ys**2) - chis * psir &
525  * (xr * xs + yr * ys)) / jacobian2
526  g12_hel(nr - i + 1, j) = grpgrc * cpsurf / radius**2
527  g33_hel(nr - i + 1, j) = ((1._r8 + eps * xaxis) &
528  / (1._r8 + eps * xchi(no)))**2
529  xout(nr - i + 1, j) = xchi(no)
530  yout(nr - i + 1, j) = ychi(no)
531  cr(nr - i + 1, j) = xchi(no)
532  cz(nr - i + 1, j) = ychi(no)
533 !-- calculate magnetic field component
534 ! br, bz normalized to rvac**2 / (bmag * rmag**2)
535 ! bphi normalized to rvac / (bmag * rmag)
536  br(nr - i + 1, j) = -psiy / (eps * (1._r8 + eps * xchi(no))) &
537  * cpsurf
538  bz(nr - i + 1, j) = psix / (eps * (1._r8 + eps * xchi(no))) &
539  * cpsurf
540  bphi(nr - i + 1, j) = rbphi(nr - i + 1) / (1._r8 + eps &
541  * xchi(no))
542 !-- check jacobian
543  grchi2 = chix**2 + chiy**2
544  dum1 = rbphi(nr - i + 1)**2 / (qs(nr - i + 1)**2 &
545  / g33_hel(nr - i + 1, j))
546  dum2 = g11_hel(nr - i + 1, j) * grchi2 / radius**2
547  dum3 = g12_hel(nr - i + 1, j)
548  dum4 = dum2 - dum3 * dum3
549  errj = abs(dum4 - dum1)
550  if (errj > maxerr) then
551  maxerr = errj
552  ierr = i
553  jerr = j
554  serr = sqrt(zpsi)
555  cerr = dchi
556  end if
557 !---------------------------------------------------------------------
558  else
559  if (standard_output) &
560  write(out_he, *) 'error in solve_cubic i,j,k : ', i, j, k, &
561  s, s2, s3
562  if (verbosity > 3) then
563  write(iu6, *) a0, a1, a2, a3, zchi
564  write(iu6, *) cchi(1, (i - 1) * np + k), cchi(1, (i - 1) &
565  * np + k + 1)
566  end if
567  end if
568  else if (chin) then
569  jbase = j
570  goto 70
571  end if
572  end do
573  70 end do
574  end do
575 
576  if (standard_output) then
577  write(out_he, *)
578  write(out_he, *) '***************************************************'
579  write(out_he, 62) maxerr, serr, cerr, ierr, jerr
580  write(out_he, *) '***************************************************'
581  end if
582 
583  if (verbosity > 2) then
584  write(iu6, *) 'a = ', a
585  write(iu6, *) 'b = ', b
586  write(iu6, *) 'alfa = ', alfa
587  end if
588 
589 !-- write normalized profiles to vector file
590  if (vec_file) then
591  write(out_vec, 11) (p0(i), rbphi(i), qs(i), i = 1, nr)
592  write(out_vec, 11) cpsurf
593  close(out_vec)
594  end if
595 
596 !-- derivatives for splines (quantity in SI units, derivative wrt s)
597  bmag = bm * bvac
598 
599 !-- calculate fields (in SI units) for output
600  rminor = eps * rvac
601  rmag = rminor / radius
602  do i = 1, nr
603  psi_out(i) = cpsurf_out * psikn(nr + 1 - i)
604  p_out(i) = p0(i) * bmag**2 / mu0
605  rbphi_out(i) = rbphi(i) * rmag * bmag
606  end do
607  psi_out(1) = 0._r8
608 
609 !-- add missing data on magnetic axis
610  cr(1, : ) = xaxis
611  cz(1, : ) = yaxis
612  g11_hel(1, : ) = 0._r8
613  g33_hel(1, : ) = 1._r8
614 !-- g12_hel has to be extrapolated
615  call allocate_spline_coefficients(coeff, nr - 1)
616  do j = 1, nchi
617  c1(1 : nr - 1) = g12_hel(2 : nr, j)
618  call spline(nr - 1, cs(2 : nr), c1, 0._r8, 0._r8, 3, coeff)
619  g12_hel(1, j) = spwert(nr - 1, 0._r8, coeff, cs(2 : nr), abltg, 0)
620  end do
621  call deallocate_spline_coefficients(coeff)
622 !-- no poloidal field on magnetic axis
623  br(1, : ) = 0._r8
624  bz(1, : ) = 0._r8
625  bphi(1, : ) = rvac / rmag
626 !-- g11_hel, g12_hel, and g33_hel should be in SI units with psi the
627 ! total poloidal flux, i.e., in Wb
628  g11_hel = g11_hel * (twopi * rmag * bmag)**2 ! [Wb^2/m^2]
629  g12_hel = g12_hel * twopi * bmag ! [Wb/m^2]
630  g33_hel = g33_hel / rmag**2 ! [m^-2]
631 !-- denormalize cr and cz
632  cr = rvac + rminor * cr ! [m]
633  cz = zgeo + rminor * cz ! [m]
634 !-- denormalize magnetic field
635  br = br * bmag * rmag**2 / rvac**2 ! [T]
636  bz = bz * bmag * rmag**2 / rvac**2 ! [T]
637  bphi = bphi * bmag * rmag / rvac ! [T]
638 
639 !-- caculate dp_hat/dpsi_hat
640  do i = 1, nr
641  if (hbt) then
642  dp_out(i) = -0.5_r8 * a * b * dp_dpsi(psikn(nr - i + 1))
643  else
644  dp_out(i) = -eps * a * b * dp_dpsi(psikn(nr - i + 1))
645  end if
646  end do
647 !-- denormalize dp_out
648  dp_out = dp_out * eps * cpsurf_out / (mu0 * rminor**4 * twopi**2)
649 
650  vx = rvac + rminor * vx
651  vy = zgeo + rminor * vy
652  xout(1, :) = rvac + rminor * xaxis
653  yout(1, :) = zgeo + rminor * yaxis
654  xout(2 : nr, :) = rvac + rminor * xout(2 : nr, :)
655  yout(2 : nr, :) = zgeo + rminor * yout(2 : nr, :)
656 
657 !-- calculate phi and rho_tor
658  call allocate_spline_coefficients(coeff, nr)
659  call spline(nr, psi_out, qs, 0._r8, 0._r8, 3, coeff)
660  equilibrium_out%profiles_1d%phi = integrate(coeff, psi_out, nr, .false.)
661  call deallocate_spline_coefficients(coeff)
662 
663 !--------- write equilibrium into equilibrium_out ----------------
664 
665  call fill_equilibrium(a, xaxis, yaxis, current_averaging, equilibrium_out)
666 
667  equilibrium_out%profiles_1d%rho_tor &
668  = sqrt(equilibrium_out%profiles_1d%phi / pi &
669  / equilibrium_out%global_param%toroid_field%b0)
670 
671 !-- calculate V' = dV/dPsi, A' = dA/dPsi, and dpsidrho_tor = dPsi/drho_tor
672  call allocate_spline_coefficients(coeff, nr)
673  call spline(nr, equilibrium_out%profiles_1d%psi, &
674  equilibrium_out%profiles_1d%volume, 0._r8, 0._r8, 3, coeff)
675  do i = 1, nr
676  dummy = spwert(nr, equilibrium_out%profiles_1d%psi(i), coeff, &
677  equilibrium_out%profiles_1d%psi, abltg, 1)
678  equilibrium_out%profiles_1d%vprime(i) = abltg(1)
679  end do
680  call spline(nr, equilibrium_out%profiles_1d%psi, &
681  equilibrium_out%profiles_1d%area, 0._r8, 0._r8, 3, coeff)
682  do i = 1, nr
683  dummy = spwert(nr, equilibrium_out%profiles_1d%psi(i), coeff, &
684  equilibrium_out%profiles_1d%psi, abltg, 1)
685  equilibrium_out%profiles_1d%aprime(i) = abltg(1)
686  end do
687  call spline(nr, equilibrium_out%profiles_1d%rho_tor, &
688  equilibrium_out%profiles_1d%psi, 0._r8, 0._r8, 3, coeff)
689  do i = 1, nr
690  dummy = spwert(nr, equilibrium_out%profiles_1d%rho_tor(i), coeff, &
691  equilibrium_out%profiles_1d%rho_tor, abltg, 1)
692  equilibrium_out%profiles_1d%dpsidrho_tor(i) = abltg(1)
693  end do
694  call deallocate_spline_coefficients(coeff)
695 
696 !-- rescale psi on axis to psi_axis
697  equilibrium_out%profiles_1d%psi = psi_axis &
698  + equilibrium_out%profiles_1d%psi - equilibrium_out%profiles_1d%psi(1)
699 
700 !-- scalar quantities
701  equilibrium_out%global_param%psi_bound = &
702  equilibrium_out%profiles_1d%psi(nr)
703  equilibrium_out%global_param%psi_ax = &
704  equilibrium_out%profiles_1d%psi(1)
705 
706 !-- 1d profiles
707  equilibrium_out%profiles_1d%q = qs
708 
709 !-- safety factors
710  equilibrium_out%global_param%mag_axis%q = equilibrium_out%profiles_1d%q(1)
711  equilibrium_out%global_param%q_min = minval(abs(equilibrium_out%profiles_1d%q)) &
712  * sign(1.0_r8, equilibrium_out%profiles_1d%q(1))
713 !-- interpolate q vs. psi to determine q_95
714  call allocate_spline_coefficients(coeff, nr)
715 
716  call spline(nr, (equilibrium_out%profiles_1d%psi - equilibrium_out%profiles_1d%psi(1)) &
717  / (equilibrium_out%profiles_1d%psi(nr) - equilibrium_out%profiles_1d%psi(1)), &
718  equilibrium_out%profiles_1d%q, 0._r8, 0._r8, 2, coeff)
719 
720  equilibrium_out%global_param%q_95 = spwert(nr, 0.95_r8, coeff, &
721  (equilibrium_out%profiles_1d%psi - equilibrium_out%profiles_1d%psi(1)) &
722  / (equilibrium_out%profiles_1d%psi(nr) - equilibrium_out%profiles_1d%psi(1)), abltg, 0)
723 
724  call deallocate_spline_coefficients(coeff)
725 
726 
727 !-- 2d profiles
728  equilibrium_out%profiles_2d(1)%grid%dim1 = equilibrium_out%profiles_1d%psi
729  do j = 1, nchi
730  equilibrium_out%profiles_2d(1)%grid%dim2(j) = (j - 1) * twopi &
731  / dble(nchi)
732  equilibrium_out%profiles_2d(1)%theta( : , j) = equilibrium_out%profiles_2d(1)%grid%dim2(j)
733  equilibrium_out%profiles_2d(1)%psi( : , j) = equilibrium_out%profiles_1d%psi
734  end do
735  equilibrium_out%profiles_2d(1)%br = br
736  equilibrium_out%profiles_2d(1)%bz = bz
737  equilibrium_out%profiles_2d(1)%bphi = bphi
738 
739 !-- boundary
740  equilibrium_out%eqgeometry%boundary(1)%r = vx
741  equilibrium_out%eqgeometry%boundary(1)%z = vy
742  equilibrium_out%eqgeometry%boundarytype = 1 ! separatrix
743 
744 !-- grid
745  equilibrium_out%coord_sys%grid_type = grid_type
746  equilibrium_out%coord_sys%grid%dim1 = equilibrium_out%profiles_1d%psi
747  do j = 1, nchi
748  equilibrium_out%coord_sys%grid%dim2(j) = (j - 1) * twopi &
749  / dble(nchi)
750  end do
751  equilibrium_out%coord_sys%position%r = cr
752  equilibrium_out%coord_sys%position%z = cz
753  equilibrium_out%profiles_2d(1)%r = cr
754  equilibrium_out%profiles_2d(1)%z = cz
755  equilibrium_out%profiles_2d(1)%grid_type = grid_type
756  equilibrium_out%coord_sys%jacobian = cr**2 / twopi
757  do i = 1, nr
758  equilibrium_out%coord_sys%jacobian(i, : ) = &
759  equilibrium_out%profiles_1d%q(i) &
760  / equilibrium_out%profiles_1d%F_dia(i) &
761  * equilibrium_out%coord_sys%jacobian(i, : )
762  end do
763  equilibrium_out%coord_sys%g_11 = g11_hel
764  equilibrium_out%coord_sys%g_12 = g12_hel
765  equilibrium_out%coord_sys%g_13 = 0._r8
766 !-- g_22 is infinity on axis, use itm_r8_invalid from itm_types
767  equilibrium_out%coord_sys%g_22(1, : ) = itm_r8_invalid
768  do i = 2, nr
769  equilibrium_out%coord_sys%g_22(i, : ) = (1._r8 &
770  / equilibrium_out%coord_sys%jacobian(i, : )**2 &
771  / g33_hel(i, : ) + g12_hel(i, : )**2) / g11_hel(i, : )
772  end do
773  equilibrium_out%coord_sys%g_23 = 0._r8
774  equilibrium_out%coord_sys%g_33 = g33_hel
775 
776 !-- calculate 1D metric coefficients gm1 through gm9
777 
778  do i = 1, nr
779  f_dia = equilibrium_out%profiles_1d%F_dia(i) &
780  * equilibrium_out%eqgeometry%a_minor &
781  / (equilibrium_out%profiles_1d%psi(nr) - psi_axis)
782  equilibrium_out%profiles_1d%gm1(i) = flux_surface_average(i, xaxis, &
783  1._r8, f_dia, 'theta ', one_over_r2) &
784  / equilibrium_out%eqgeometry%geom_axis%r**2
785  equilibrium_out%profiles_1d%gm4(i) = flux_surface_average(i, xaxis, &
786  1._r8, f_dia, 'theta ', one_over_b2) &
787  * (equilibrium_out%eqgeometry%a_minor &
788  * equilibrium_out%eqgeometry%geom_axis%r &
789  / (equilibrium_out%profiles_1d%psi(nr) - psi_axis))**2
790  equilibrium_out%profiles_1d%gm5(i) = flux_surface_average(i, xaxis, &
791  1._r8, f_dia, 'theta ', b2) / (equilibrium_out%eqgeometry%a_minor &
792  * equilibrium_out%eqgeometry%geom_axis%r &
793  / (equilibrium_out%profiles_1d%psi(nr) - psi_axis))**2
794  equilibrium_out%profiles_1d%gm8(i) = flux_surface_average(i, xaxis, &
795  1._r8, f_dia, 'theta ', r_major) &
796  * equilibrium_out%eqgeometry%geom_axis%r
797  equilibrium_out%profiles_1d%gm9(i) = flux_surface_average(i, xaxis, &
798  1._r8, f_dia, 'theta ', one_over_r) &
799  / equilibrium_out%eqgeometry%geom_axis%r
800  if (i > 1) then
801  equilibrium_out%profiles_1d%gm2(i) = flux_surface_average(i, xaxis, &
802  1._r8, f_dia, 'theta ', gradpsi2_over_r2) &
803  / (equilibrium_out%eqgeometry%a_minor &
804  * equilibrium_out%eqgeometry%geom_axis%r &
805  / (equilibrium_out%profiles_1d%psi(nr) - psi_axis))**2 &
806  * (equilibrium_out%profiles_1d%q(i) &
807  / (equilibrium_out%profiles_1d%rho_tor(i) &
808  * equilibrium_out%global_param%toroid_field%b0))**2
809  equilibrium_out%profiles_1d%gm3(i) = flux_surface_average(i, xaxis, &
810  1._r8, f_dia, 'theta ', gradpsi2) &
811  / (equilibrium_out%eqgeometry%a_minor &
812  / (equilibrium_out%profiles_1d%psi(nr) - psi_axis))**2 &
813  * (equilibrium_out%profiles_1d%q(i) &
814  / (equilibrium_out%profiles_1d%rho_tor(i) &
815  * equilibrium_out%global_param%toroid_field%b0))**2
816  equilibrium_out%profiles_1d%gm6(i) = flux_surface_average(i, xaxis, &
817  1._r8, f_dia, 'theta ', gradrho2_over_b2) &
818  * equilibrium_out%eqgeometry%geom_axis%r**2 &
819  * (equilibrium_out%profiles_1d%q(i) &
820  / (equilibrium_out%profiles_1d%rho_tor(i) &
821  * equilibrium_out%global_param%toroid_field%b0))**2
822  equilibrium_out%profiles_1d%gm7(i) = flux_surface_average(i, xaxis, &
823  1._r8, f_dia, 'theta ', gradpsi) &
824  * (equilibrium_out%profiles_1d%psi(nr) &
825  - psi_axis) / equilibrium_out%eqgeometry%a_minor &
826  * equilibrium_out%profiles_1d%q(i) &
827  / (equilibrium_out%profiles_1d%rho_tor(i) &
828  * equilibrium_out%global_param%toroid_field%b0)
829  end if
830 ! warning: added by DPC but I'm not a HELENA expert
831  equilibrium_out%profiles_1d%b_av(i) = flux_surface_average(i, xaxis, &
832  1._r8, f_dia, 'theta ', b_fun) / abs(equilibrium_out%eqgeometry%a_minor &
833  * equilibrium_out%eqgeometry%geom_axis%r &
834  / (equilibrium_out%profiles_1d%psi(nr) - psi_axis))
835  equilibrium_out%profiles_1d%b_min(i) = minval(sqrt(br(i,:)*br(i,:)+bz(i,:)*bz(i,:)+bphi(i,:)*bphi(i,:)))
836  equilibrium_out%profiles_1d%b_max(i) = maxval(sqrt(br(i,:)*br(i,:)+bz(i,:)*bz(i,:)+bphi(i,:)*bphi(i,:)))
837  end do
838 !-- extrapolate gm2, gm3, gm6, and gm7 onto axis
839  call allocate_spline_coefficients(coeff, nr - 1)
840  call spline(nr - 1, psi_out(2 : nr), &
841  equilibrium_out%profiles_1d%gm2(2 : nr), 0._r8, 0._r8, 2, coeff)
842  equilibrium_out%profiles_1d%gm2(1) = spwert(nr - 1, 0._r8, coeff, &
843  psi_out(2 : nr), abltg, 0)
844  call spline(nr - 1, psi_out(2 : nr), &
845  equilibrium_out%profiles_1d%gm3(2 : nr), 0._r8, 0._r8, 2, coeff)
846  equilibrium_out%profiles_1d%gm3(1) = spwert(nr - 1, 0._r8, coeff, &
847  psi_out(2 : nr), abltg, 0)
848  call spline(nr - 1, psi_out(2 : nr), &
849  equilibrium_out%profiles_1d%gm6(2 : nr), 0._r8, 0._r8, 2, coeff)
850  equilibrium_out%profiles_1d%gm6(1) = spwert(nr - 1, 0._r8, coeff, &
851  psi_out(2 : nr), abltg, 0)
852  call spline(nr - 1, psi_out(2 : nr), &
853  equilibrium_out%profiles_1d%gm7(2 : nr), 0._r8, 0._r8, 2, coeff)
854  equilibrium_out%profiles_1d%gm7(1) = spwert(nr - 1, 0._r8, coeff, &
855  psi_out(2 : nr), abltg, 0)
856  call deallocate_spline_coefficients(coeff)
857 
858 !-- calculation of r_inboard
859  if (ias == 0) then
860  r_in = xout( : , nchi)
861  else
862  call r_inboard
863  end if
864 
865 !-- calculation of elongation and triangularities
866  equilibrium_out%profiles_1d%elongation(1) = 1._r8 ! magnetic axis
867  equilibrium_out%profiles_1d%tria_upper(1) = 0._r8 ! magnetic axis
868  equilibrium_out%profiles_1d%tria_lower(1) = 0._r8 ! magnetic axis
869  do i = 2, nr
870  r_min_flux_surface = r_min(rminor, xaxis, i)
871  r_max_flux_surface = r_max(rminor, xaxis, i)
872  a_flux_surface = (r_max_flux_surface - r_min_flux_surface) / 2._r8
873  r_geo_flux_surface = (r_max_flux_surface + r_min_flux_surface) / 2._r8
874  equilibrium_out%profiles_1d%elongation(i) = equilibrium_out%profiles_1d%area(i) &
875  / (pi * a_flux_surface**2)
876  top_point = rz_top(rminor, xaxis, yaxis, i)
877  bottom_point = rz_bottom(rminor, xaxis, yaxis, i)
878  delta_upper = r_geo_flux_surface - top_point%R
879  delta_lower = r_geo_flux_surface - bottom_point%R
880  equilibrium_out%profiles_1d%tria_upper(i) = delta_upper / a_flux_surface
881  equilibrium_out%profiles_1d%tria_lower(i) = delta_lower / a_flux_surface
882  end do
883  equilibrium_out%eqgeometry%elongation = equilibrium_out%profiles_1d%elongation(nr)
884  equilibrium_out%eqgeometry%tria_upper = equilibrium_out%profiles_1d%tria_upper(nr)
885  equilibrium_out%eqgeometry%tria_lower = equilibrium_out%profiles_1d%tria_lower(nr)
886 
887  equilibrium_out%profiles_1d%r_inboard = r_in
888  equilibrium_out%profiles_1d%r_outboard = xout( : , 1)
889 
890 !TODO: This is legacy code. Should be removed after testing is completed.
891  if (trim(output) == 'equilibrium' .or. trim(output) == 'full') &
892  call write_old_equilibrium(path, equilibrium_out)
893 
894  if (xmgrace_output) then
895  call plot_data_1d('line', equilibrium_out%profiles_1d%psi,&
896  equilibrium_out%profiles_1d%F_dia, nr, 'xmgr_final_F_psi')
897  call plot_data_1d('line', psi_out, p_out, nr, 'xmgr_final_p_psipol')
898  call plot_data_1d('line', psi_out, dp_out, nr, &
899  'xmgr_final_dp_dpsi_psipol')
900  work1 = dp0 * bmag**2 / mu0
901  call plot_data_1d('line', psi_out, work1, nr, &
902  'xmgr_final_dp_drho_psipol')
903  work1 = drbphi * rmag * bmag
904  call plot_data_1d('line', psi_out, work1, nr, &
905  'xmgr_final_drbphi_drho_psipol')
906  call plot_data_1d('line', psi_out, qs, nr, 'xmgr_q_psipol')
907  call plot_data_1d('line', vx, vy, nchi, 'xmgr_boundary_x_y_real')
908  call plot_data_2d('2d', xout, yout, 'xmgr_final_flux_surfaces_real')
909  call plot_data_2d('2d_t', xout, yout, 'xmgr_final_radial_real')
910  call plot_data_2d('grid', xout, yout, 'xmgr_final_grid_x_y_real')
911  call plot_data_1d('line', psi_out, p0, nr, 'xmgr_p_bar_psipol')
912  call plot_data_1d('line', psi_out, rbphi, nr, 'xmgr_rbphi_bar_psipol')
913  call plot_data_2d('3d', cr, cz, 'idl_g11_x_y', g11_hel)
914  call plot_data_2d('3d', cr, cz, 'idl_g12_x_y', g12_hel)
915  call plot_data_2d('3d', cr, cz, 'idl_g33_x_y', g33_hel)
916  call plot_data_2d('3d', cr, cz, 'idl_br_x_y', br)
917  call plot_data_2d('3d', cr, cz, 'idl_bz_x_y', bz)
918  call plot_data_2d('3d', cr, cz, 'idl_bphi_x_y', bphi)
919  psi_out_2d(:, 1) = psi_out(:)
920  call plot_data_2d('line', psi_out_2d, g11_hel, 'xmgr_g11_spol')
921  call plot_data_2d('line', psi_out_2d, g12_hel, 'xmgr_g12_spol')
922  call plot_data_2d('line', psi_out_2d, g33_hel, 'xmgr_g33_spol')
923  end if
924 
925  if (profiles_output) then
926 !-- write functions and their derivatives
927  open (unit = out_dat, file = trim(adjustl(path)) // file_dat, &
928  status = 'replace', form = 'formatted', action = 'write', &
929  iostat = i_error)
930  do i = 1, n_string
931  header_string((i - 1) * len_string + 1 : i * len_string) = &
932  string(i)(1 : len_string)
933  end do
934  write(out_dat, '(416a)') header_string
935 
936  rho(1) = sqrt(psikn(nr))
937 !-- calculation of drbphi
938  do i = 2, nr
939  rho(i) = sqrt(psikn(nr + 1 - i))
940  drbphi_out(i) = drbphi(i) * rmag * bmag / (2._r8 * rho(i)) &
941  / cpsurf_out
942  end do
943 !-- extrapolation to axis
944  call allocate_spline_coefficients(coeff, nr - 1)
945  call spline(nr - 1, psi_out(2 : nr), drbphi_out(2 : nr), 0._r8, &
946  0._r8, 2, coeff)
947  drbphi_out(1) = spwert(nr - 1, 0._r8, coeff, psi_out(2 : nr), abltg, 0)
948  call deallocate_spline_coefficients(coeff)
949 !-- calculation of R_geo and rho_vol
950  r_geo = 0.5_r8 * (r_in + xout( : , 1))
951  rho_vol = sqrt(equilibrium_out%profiles_1d%volume / (2._r8 * pi**2 &
952  * r_geo))
953 !TODO: I_tor from current_calculations could be used but is not allocated
954 !-- calculate I_tor
955  i_tor(1) = equilibrium_out%profiles_1d%jphi(1) &
956  * equilibrium_out%profiles_1d%area(1)
957  do i = 2, nr
958  i_tor(i) = i_tor(i - 1) + (equilibrium_out%profiles_1d%area(i) &
959  - equilibrium_out%profiles_1d%area(i - 1)) &
960  * equilibrium_out%profiles_1d%jphi(i)
961  end do
962 !-- calculate dq/dpsi and dV/dpsi
963  call allocate_spline_coefficients(coeff, nr)
964  call spline(nr, psi_out, qs, 0._r8, 0._r8, 2, coeff)
965  do i = 1, nr
966  psi_dummy = spwert(nr, psi_out(i), coeff, psi_out, abltg, 1)
967  dq_dpsi(i) = twopi * abltg(1)
968  end do
969  call spline(nr, psi_out, equilibrium_out%profiles_1d%volume, &
970  0._r8, 0._r8, 2, coeff)
971  do i = 1, nr
972  psi_dummy = spwert(nr, psi_out(i), coeff, psi_out, abltg, 1)
973  dv_dpsi(i) = twopi * abltg(1)
974  end do
975  call deallocate_spline_coefficients(coeff)
976 !-- calculate stability alpha and global shear (following Miller et al.,
977 ! Phys. Plasmas 5(4), 1998)
978  stability_alpha = -dv_dpsi / pi * rho_vol * mu0 * dp_out
979  global_shear = 2._r8 * equilibrium_out%profiles_1d%volume &
980  * dq_dpsi / dv_dpsi
981 
982  do i = 1, nr
983  psi_bar = psikn(nr + 1 - i)
984  write(out_dat, 16) psi_bar, rho(i), psi_out(i) / twopi, psi_out(i),&
985  rmag - r_in(i), r_in(i), xout(i, 1) - rmag, xout(i, 1), r_geo(i), &
986  rho_vol(i), equilibrium_out%profiles_1d%rho_vol(i), rho_vol(i) &
987  * sqrt(r_geo(i) / r_geo(nr)), equilibrium_out%profiles_1d%phi(i), &
988  equilibrium_out%profiles_1d%rho_tor(i), &
989  equilibrium_out%profiles_1d%area(i), &
990  equilibrium_out%profiles_1d%volume(i), dv_dpsi(i), &
991  p_out(i), twopi * dp_out(i), dp_out(i), rbphi_out(i), twopi &
992  * drbphi_out(i), drbphi_out(i), rbphi_out(i) * drbphi_out(i), &
993  qs(i), dq_dpsi(i), equilibrium_out%profiles_1d%jphi(i), &
994  i_tor(i), stability_alpha(i), global_shear(i), &
995  equilibrium_out%profiles_1d%gm1(i), &
996  equilibrium_out%profiles_1d%gm2(i), &
997  equilibrium_out%profiles_1d%gm3(i), &
998  equilibrium_out%profiles_1d%gm4(i), &
999  equilibrium_out%profiles_1d%gm5(i), &
1000  equilibrium_out%profiles_1d%gm6(i), &
1001  equilibrium_out%profiles_1d%gm7(i), &
1002  equilibrium_out%profiles_1d%gm8(i), &
1003  equilibrium_out%profiles_1d%gm9(i)
1004  end do
1005 
1006  close (out_dat)
1007 
1008  end if
1009 
1010  if (elite_output) then
1011 
1012 !-- normalize g11_hel, g12_hel, and g33_hel
1013  g11_hel = g11_hel / (twopi * rmag * bmag )**2 ! [Wb^2/m^2]
1014  g12_hel = g12_hel / (twopi * bmag) ! [Wb/m^2]
1015  g33_hel = g33_hel * rmag**2 ! [m^-2]
1016 
1017 !-- write elite input
1018  call profiles(p0tmp, ftmp, dptmp, dftmp, a)
1019 
1020 !-- ELITE uses poloidal flux per radian instead of total poloidal flux
1021  cpsurf_out = cpsurf_out / twopi
1022  psi_out = psi_out / twopi
1023 
1024  if (verbosity > 3) &
1025  write(iu6, *) 'elite: cpsurf = ', cpsurf, ' , cpsurf_out = ', &
1026  cpsurf_out
1027  pscale0 = pscale
1028  pscale = bvac**2 * eps / alfa**2 / cpsurf_out
1029  rbscale0 = rbscale
1030  rbscale = bvac * rvac
1031 
1032  do i = 1, nr
1033  pskn1(i) = psikn(nr + 1 - i)
1034  end do
1035 
1036  do i = 2, nr
1037  ppelite(i) = dptmp(i) * pscale / 2._r8 / sqrt(psikn(nr + 1 - i))
1038  felite(i) = ftmp(i) * rbscale
1039  ffpelite(i) = dftmp(i) * ftmp(i) * rbscale**2 / cpsurf_out / 2._r8 &
1040  / sqrt(psikn(nr + 1 - i))
1041  end do
1042 
1043  felite(1) = ftmp(1) * rbscale
1044  ppelite(1) = ppelite(2) - (ppelite(2) - ppelite(3)) &
1045  / (psi_out(3) - psi_out(2)) * (psi_out(2) - psi_out(1))
1046  ffpelite(1) = ffpelite(2) - (ffpelite(2) - ffpelite(3)) &
1047  / (psi_out(3) - psi_out(2)) * (psi_out(2) - psi_out(1))
1048 
1049 ! xdum = 0._r8
1050 ! ydum = (ppelite(nr) - ppelite(nr - 1)) / (psi_out(nr) - psi_out(nr - 1))
1051 ! call allocate_spline_coefficients(dd_spline, nr)
1052 ! call spline(nr, pskn1, ppelite, xdum, ydum, 1, dd_spline)
1053 ! xdum = (ffpelite(2) - ffpelite(1)) / (psi_out(2) - psi_out(1))
1054 ! ydum = (ffpelite(nr) - ffpelite(nr - 1)) / (psi_out(nr) - psi_out(nr &
1055 ! - 1))
1056 ! call allocate_spline_coefficients(f_spline, nr)
1057 ! call spline(nr, pskn1, ffpelite, xdum, ydum, 1, f_spline)
1058 
1059 ! Samuli's final ffprime and f^2 profiles
1060 
1061  if (xmgrace_output) then
1062  call plot_data_1d('line', psi_out, ppelite, nr, &
1063  'xmgr_elite_pprime_realpsipol')
1064  call plot_data_1d('line', psi_out, felite, nr, &
1065  'xmgr_elite_f_realpsipol')
1066  call plot_data_1d('line', psi_out, ffpelite, nr, &
1067  'xmgr_elite_ffprime_realpsipol')
1068  end if
1069 
1070  do i = 2, nr - 1
1071  pppelite(i) = ((ppelite(i) - ppelite(i - 1)) / (psi_out(i) &
1072  - psi_out(i - 1)) * (psi_out(i + 1) - psi_out(i)) &
1073  + (ppelite(i + 1) - ppelite(i)) / (psi_out(i + 1) &
1074  - psi_out(i)) * (psi_out(i) - psi_out(i - 1))) &
1075  / (psi_out(i + 1) - psi_out(i - 1))
1076  fppelite(i) = ((ffpelite(i) - ffpelite(i - 1)) / (psi_out(i) &
1077  - psi_out(i - 1)) * (psi_out(i + 1) - psi_out(i)) &
1078  + (ffpelite(i + 1) - ffpelite(i)) / (psi_out(i + 1) &
1079  - psi_out(i)) * (psi_out(i) - psi_out(i - 1))) &
1080  / (psi_out(i + 1) - psi_out(i - 1))
1081 ! xdum = spwert(nr, psi_out(i), dd_spline, pskn1, abltg, 1)
1082 ! pppelite(i) = abltg(1)
1083 ! xdum = spwert(nr, psi_out(i), f_spline, pskn1, abltg, 1)
1084 ! fppelite(i) = abltg(1)
1085 ! if (verbosity > 3) &
1086 ! write(iu6, *) psi_out(i), fppelite(i), (ffpelite(i) &
1087 ! - ffpelite(i - 1)) / (psi_out(i) - psi_out(i - 1))
1088  end do
1089 
1090 ! call deallocate_spline_coefficients(dd_spline)
1091 ! call deallocate_spline_coefficients(f_spline)
1092 
1093 ! pppelite(nr) = 2 * pppelite(nr - 1) - pppelite(nr - 2)
1094 ! fppelite(nr) = 2 * fppelite(nr - 1) - fppelite(nr - 2)
1095 ! pppelite(nr) = (ppelite(nr) - ppelite(nr - 1)) / (psi_out(nr) &
1096 ! - psi_out(nr - 1))
1097 ! fppelite(nr) = (ffpelite(nr) - ffpelite(nr - 1)) / (psi_out(nr) &
1098 ! - psi_out(nr - 1))
1099 
1100  pppelite(1) = pppelite(2) - (pppelite(2) - pppelite(3)) &
1101  / (psi_out(3) - psi_out(2)) * (psi_out(2) - psi_out(1))
1102  fppelite(1) = fppelite(2) - (fppelite(2) - fppelite(3)) &
1103  / (psi_out(3) - psi_out(2)) * (psi_out(2) - psi_out(1))
1104  fppelite(nr) = fppelite(nr - 1) - (fppelite(nr - 1) &
1105  - fppelite(nr - 2)) / (psi_out(nr - 2) - psi_out(nr - 1)) &
1106  * (psi_out(nr - 1) - psi_out(nr))
1107  pppelite(nr) = pppelite(nr - 1) - (pppelite(nr - 1) &
1108  - pppelite(nr - 2)) / (psi_out(nr - 2) - psi_out(nr - 1)) &
1109  * (psi_out(nr - 1) - psi_out(nr))
1110 ! pppelite(1) = 2 * pppelite(2) - pppelite(3)
1111 ! fppelite(1) = 2 * fppelite(2) - fppelite(3)
1112 
1113  if (verbosity > 2) then
1114  write(iu6, *) 'writing elite input'
1115  write(iu6, *) 'alfa =', alfa
1116  write(iu6, *) 'b0, minor radius, rvac, axis =', b0, eps * rvac, &
1117  rvac, raxis / (1 + eps * xaxis)
1118  write(iu6, *) 'rmag =', rmag
1119  write(iu6, *) 'cpsurf_out, psi1 = ', cpsurf_out, (rvac * eps)**2 &
1120  * bvac / alfa
1121  end if
1122 
1123 ! qs(1) = rbphi(1) * pi / (2. * sqrt(cx * cy) * (1. + eps * xaxis) &
1124 ! * rbscale)
1125 
1126 ! calculate q using bp
1127 
1128  omega = 0.0_r8
1129  rold = vx(1)
1130  zold = vy(1)
1131  bpold = sqrt(g11_hel(nr, 1)) * cpsurf_out / cpsurf / (vx(1) * rmag)
1132 
1133  leng = 0
1134  do i = j, nchi
1135  r = (1._r8 + xx(1, j) * eps) * rvac
1136  z = zgeo + (yy(1, j) * eps) * rvac
1137 ! r = vx(j)
1138 ! z = vy(j)
1139  dl = sqrt((r - rold)**2 + (z - zold)**2)
1140  bp = sqrt(g11_hel(nr, j)) * cpsurf_out / cpsurf &
1141  / (vx(j) * rmag)
1142  omega = omega + 0.5_r8 * dl * (1._r8 / (r**2 * bp) + 1._r8 &
1143  / (rold**2 * bpold))
1144 ! if (verbosity > 3) write(iu6, *) omega, r, bp
1145  rold = r
1146  zold = z
1147  bpold = bp
1148  end do
1149  qcalc = factas * omega * felite(nr) / pi / 2._r8
1150  qnorm = qcalc / qs(nr)
1151  if (verbosity > 2) &
1152  write(iu6, *) 'q_calc = ', omega * felite(nr) / pi / 2._r8, &
1153  'q_helena = ', qs(nr)
1154 
1155  open (unit = out_elite, file = trim(adjustl(path)) // file_elite, &
1156  status = 'replace', form = 'formatted', action = 'write', &
1157  iostat = i_error)
1158  open (unit = out_elite_b, file = trim(adjustl(path)) // file_elite_b, &
1159  status = 'replace', form = 'formatted', action = 'write', &
1160  iostat = i_error)
1161  open (unit = out_elite_b9, file = trim(adjustl(path)) // file_elite_b9, &
1162  status = 'replace', form = 'formatted', action = 'write', &
1163  iostat = i_error)
1164 
1165  write(out_elite, 1000) 'helena generated input'
1166  if (ias == 1) then
1167  write(out_elite, 1001) nr, nchi + 1
1168  else
1169  write(out_elite, 1001) nr, 2 * nchi - 1
1170  end if
1171  write(out_elite, *) 'psi:'
1172  write(out_elite, 1002) (psi_out(i), i = 1, nr)
1173  write(out_elite, *) 'dp/dpsi:'
1174  write(out_elite, 1002) (ppelite(i), i = 1, nr)
1175  write(out_elite, *) 'd2p/dpsi:'
1176  write(out_elite, 1002) (pppelite(i), i = 1, nr)
1177  write(out_elite, *) 'fpol:'
1178  write(out_elite, 1002) (felite(i), i = 1, nr)
1179  write(out_elite, *) 'ffp:'
1180  write(out_elite, 1002) (ffpelite(i), i = 1, nr)
1181  write(out_elite, *) 'dffp:'
1182  write(out_elite, 1002) (fppelite(i), i = 1, nr)
1183  write(out_elite, *) 'q:'
1184  write(out_elite, 1002) (qs(i), i = 1, nr)
1185  write(out_elite, *) 'r:'
1186  do j = 0, nchi - 1
1187  write(out_elite, 1002) (xout(i, j + 1), i = 1, nr)
1188  end do
1189  if (ias == 1) then
1190  write(out_elite, 1002) (xout(i, 1), i = 1, nr)
1191  else
1192  do j = nchi - 2, 0, -1
1193  write(out_elite, 1002) (xout(i, j + 1), i = 1, nr)
1194  end do
1195  end if
1196  write(out_elite, *) 'z:'
1197  if (verbosity > 3) &
1198  write(iu6, *) 'xaxis, yaxis=', xaxis, yaxis
1199  do j = 0, nchi - 1
1200  write(out_elite, 1002) (yout(i, j + 1), i = 1, nr)
1201  end do
1202  if (ias == 1) then
1203  write(out_elite, 1002) (yout(i, 1), i = 1, nr)
1204  else
1205  do j = nchi - 2, 0, -1
1206  write(out_elite, 1002) (-yout(i, j + 1), i = 1, nr)
1207  end do
1208  end if
1209  jacobian = xr * ys - xs * yr
1210  write(out_elite, *) 'bp:' ! Bpol
1211  do j = 0, nchi - 1
1212  g11_hel(1, j + 1) = 0._r8
1213  write(out_elite, 1002) &
1214  (abs(sqrt(g11_hel(i, j + 1)) * cpsurf_out / cpsurf &
1215  / (cr(i, j + 1) * rmag)), i = 1, nr)
1216  end do
1217  if (ias == 1) then
1218  write(out_elite, 1002) &
1219  (abs(sqrt(g11_hel(i, 1)) * cpsurf_out / cpsurf &
1220  / (cr(i, 1) * rmag)), i = 1, nr)
1221  else
1222  do j = nchi - 2, 0, -1
1223  write(out_elite, 1002) &
1224  (abs(sqrt(g11_hel(i, j + 1)) * cpsurf_out / cpsurf &
1225  / (cr(i, j + 1) * rmag)), i = 1, nr)
1226  end do
1227  end if
1228 
1229 ! write dummy ne, te and ti
1230 
1231  if (te%shape == 0) then
1232  write(out_elite, *) 'ne:'
1233  write(out_elite, 1002) (3.0d19, i = 1, nr)
1234  write(out_elite, *)'dne/dpsi:'
1235  write(out_elite, 1002) (0.0, i = 1, nr)
1236  write(out_elite, *) 'te:'
1237  write(out_elite, 1002) (1.0, i = 1, nr)
1238  write(out_elite, *) 'dte/dpsi:'
1239  write(out_elite, 1002) (0.0, i = 1, nr)
1240  write(out_elite, *) 'ti:'
1241  write(out_elite, 1002) (1.0, i = 1, nr)
1242  write(out_elite, *) 'dti/dpsi:'
1243  write(out_elite, 1002) (0.0, i = 1, nr)
1244  else
1245  call ne_te_profile(ne, ne%psi_0, defte, dndpsico, consde)
1246  call ne_te_profile(te, te%psi_0, tefte, dtdpsico, conste)
1247 
1248  do i = 1, nr
1249  ps = pskn1(i)
1250  if (ps >= ne%psi_0) then
1251  call ne_te_profile(ne, ps, nne(i), dnne(i))
1252  else
1253  call ne_te_profile(ne, ps, nne(i), dnne(i), consde, dndpsico)
1254  end if
1255  dnne(i) = -dnne(i)
1256 ! te given in ev (profiles defined in kev)
1257  if (ps >= te%psi_0) then
1258  call ne_te_profile(te, ps, tte(i), dtte(i))
1259  else
1260  call ne_te_profile(te, ps, tte(i), dtte(i), conste, dtdpsico)
1261  end if
1262  dtte(i) = -dtte(i) * 1000.0_r8
1263  dtte(i) = dtte(i) / cpsurf_out
1264  dnne(i) = dnne(i) / cpsurf_out * 1.e19_r8
1265  nne(i) = nne(i) * 1.e19_r8
1266  end do
1267 
1268  write(out_elite, *) 'ne:'
1269  write(out_elite, 1002) (nne(i), i = 1, nr)
1270  write(out_elite, *) 'dne/dpsi:'
1271  write(out_elite, 1002) (dnne(i), i = 1, nr)
1272  write(out_elite, *) 'te:'
1273  write(out_elite, 1002) (tte(i), i = 1, nr)
1274  write(out_elite, *) 'dte/dpsi:'
1275  write(out_elite, 1002) (dtte(i), i = 1, nr)
1276 ! assume ti = te
1277  write(out_elite, *) 'ti:'
1278  write(out_elite, 1002) (tte(i), i = 1, nr)
1279  write(out_elite, *) 'dti/dpsi:'
1280  write(out_elite, 1002) (dtte(i), i = 1, nr)
1281 
1282  end if
1283 
1284  write(out_elite_b, 1003) &
1285  (xout(nr, j), yout(nr, j), j = 1, nchi)
1286 
1287  nedgelab = int(edgelabe * (dble(nr - 1)) / 100._r8)
1288  if (verbosity > 3) &
1289  write(iu6, *) 'writing surface: ', nedgelab
1290  if (nedgelab >= nr) then
1291  nedgelab = nr - 1
1292  end if
1293  write(out_elite_b9, 1003) &
1294  (xout(nedgelab, j), yout(nedgelab, j), j = 1, nchi)
1295  pscale = pscale0
1296  rbscale = rbscale0
1297 
1298  close(out_elite)
1299  close(out_elite_b)
1300  close(out_elite_b9)
1301 
1302  end if
1303 
1304  if (diagnostics_on) then
1305 !---------------------------------------- plot profiles --------------
1306  do i = 1, nr
1307  if (ias == 1) then
1308  ileft = (i - 1) * np + (np + 1) / 2
1309  else
1310  ileft = (i - 1) * np + np
1311  end if
1312  iright = (i - 1) * np + 1
1313  xleft = xx(1, ileft)
1314  xright = xx(1, iright)
1315  zpsi = psikn(i)
1316  if (hbt) then
1317  cuplot(i) = dgamma_dpsi(zpsi) + b * xleft * (1._r8 + eps &
1318  * xleft / 2._r8) * dp_dpsi(zpsi)
1319  cuplot(2 * nr - i) = dgamma_dpsi(zpsi) + b * xright * (1._r8 &
1320  + eps * xright / 2._r8) * dp_dpsi(zpsi)
1321  else
1322  cuplot(i) = dgamma_dpsi(zpsi) + b * (1._r8 + eps * xleft)**2 &
1323  * dp_dpsi(zpsi)
1324  cuplot(2 * nr - i) = dgamma_dpsi(zpsi) + b * (1._r8 + eps &
1325  * xright)**2 * dp_dpsi(zpsi)
1326  end if
1327  cuplot(i) = a * cuplot(i) / (1._r8 + eps * xleft)
1328  if (i /= nr) then
1329  cuplot(2 * nr - i) = a * cuplot(2 * nr - i) / (1._r8 + eps &
1330  * xright)
1331  end if
1332  end do
1333  do i = 1, nr
1334  if (ias == 1) then
1335  ileft = (i - 1) * np + (np + 1) / 2
1336  else
1337  ileft = (i - 1) * np + np
1338  end if
1339  iright= (i - 1) * np + 1
1340  xplot(i) = xx(1, ileft)
1341  xplot(2 * nr - i) = xx(1, iright)
1342  pplot(i) = p0(nr - i + 1) / p0(1)
1343  pplot(2 * nr - i) = p0(nr - i + 1) / p0(1)
1344  psiplot(i) = psikn(i)
1345  psiplot(2 * nr - i) = psikn(i)
1346  qplot(2 * nr - i) = qplot(i)
1347  dqplot(2 * nr - i) = dqplot(i)
1348  end do
1349  if (standard_output) then
1350  write(out_he, *)
1351  write(out_he, *) '***************************************************'
1352  write(out_he, *) '* i, x, psi, p, q *'
1353  write(out_he, *) '***************************************************'
1354  write(out_he, 101) (i, xplot(i), psiplot(i), pax * pplot(i), &
1355  qplot(i), i = 1, 2 * nr - 1)
1356  write(out_he, *) '***************************************************'
1357  write(out_he, *)
1358  write(out_he, *) '***********************************************'
1359  write(out_he, *) ' toroidal current density on the midplane'
1360  write(out_he, *) '***********************************************'
1361  write(out_he, *) ' x jz '
1362  write(out_he, 77) (xplot(i), cuplot(i), i = 1, 2 * nr - 1)
1363  write(out_he, *)
1364  end if
1365 
1366  if (additional_output) then
1367  write(25, *) ' i, x, psi, p, q '
1368  write(25, 106) (xplot(i), psiplot(i), paxis * pplot(i), qplot(i), &
1369  i = 1, 2 * nr - 1)
1370  end if
1371 
1372  end if
1373 
1374  2 format(' radius : ',e12.4)
1375  3 format(' b0 : ',e12.4)
1376  4 format(' cpsurf : ',e12.4)
1377  11 format(3e16.8)
1378  16 format(39es16.8)
1379  31 format(' q on axis = ', f7.4)
1380  32 format(' q at boundary = ', f7.4)
1381  61 format(8e16.8)
1382  62 format(' max. error in jacobian after mapping : ', &
1383  es10.3, 2f10.6, 2i4)
1384  77 format(2e12.4)
1385  91 format(4e16.8)
1386 101 format(i4,4e12.4)
1387 106 format(4e12.4)
1388 1000 format(a30)
1389 1001 format(2i5)
1390 1002 format(5g20.12)
1391 1003 format(2g20.12)
1392 
1393  return
1394 
1395  contains
1396 
1397  subroutine r_inboard
1398 
1399 !-----------------------------------------------------------------------
1400 ! subroutine to calculate values of R in the midplane inboard
1401 !-----------------------------------------------------------------------
1402 
1403  real(r8), dimension(nchi) :: x_val, y_val
1404  integer(itm_i4) :: i_start, i_end, i_diff, i_1, i_2, ind(1), j
1405  integer(itm_i4) :: n_spline
1406  real(r8) :: zaxis
1407 
1408  zaxis = yout(1, 1)
1409  i_start = nchi / 4
1410  i_end = (3 * nchi) / 4
1411  i_diff = i_end - i_start + 1
1412  n_spline = nchi / 5
1413 
1414  call allocate_spline_coefficients(coeff, n_spline + 1)
1415 
1416  r_in(1) = xout(1, 1)
1417 
1418  do j = 2, nr
1419  y_val(1 : i_diff) = abs(yout(j, i_start : i_end) - zaxis )
1420  ind = minloc(y_val(1 : i_diff)) + i_start - 1
1421  i_1 = ind(1) - n_spline / 2
1422  i_2 = ind(1) + n_spline / 2
1423  x_val(1 : n_spline + 1) = xout(j, i_2 : i_1 : -1)
1424  y_val(1 : n_spline + 1) = yout(j, i_2 : i_1 : -1)
1425  call spline(n_spline + 1, y_val, x_val, 0._r8, 0._r8, 2, coeff)
1426  r_in(j) = spwert(n_spline + 1, zaxis, coeff, y_val, abltg, 0)
1427  end do
1428  call deallocate_spline_coefficients(coeff)
1429  return
1430  end subroutine r_inboard
1431 
1432  end subroutine mapping
1433 
1434 
1435 !-----------------------------------------------------------------------
1436  function r_min(rminor, xaxis, i) result(f_R_min)
1437 !-----------------------------------------------------------------------
1438 ! function to determine the smallest value of the major radius
1439 ! on a flux surface
1440 ! dependencies: xx
1441 !-----------------------------------------------------------------------
1442 
1443  use math_functions, only : root
1444  use mod_dat, only : rvac
1445  use mod_helena_io
1446 
1447  implicit none
1448 
1449  integer(itm_i4), intent(in) :: i
1450  real(r8), intent(in) :: rminor, xaxis
1451 
1452  real(r8) :: f_r_min
1453 
1454  integer(itm_i4) :: i_tilde
1455  integer(itm_i4) :: n1, n2
1456  integer(itm_i4) :: j, j_min
1457  integer(itm_i4) :: j_next, j_previous
1458  integer(itm_i4) :: ifail
1459  real(r8) :: x_min, s_min
1460  real(r8) :: c1, c2, c3, c4
1461  real(r8) :: aa, bb, cc, det ! annotation as in solve_cubic
1462  real(r8), parameter :: tol = 1.e-8_r8
1463 
1464 
1465  if (ias == 0) then ! up-down-symmetric case
1466 
1467  i_tilde = nr - i + 1
1468  j_min = i_tilde * np ! last point in half circle = inboard side
1469  x_min = xx(1, j_min)
1470 
1471  else
1472 
1473  if (i == 1) then ! magnetic axis
1474  x_min = xaxis
1475  else
1476 
1477 !-- radial index runs from boundary to axis
1478  i_tilde = nr - i + 1
1479 
1480 !-- search for minimum in xx(1, (i_tilde - 1) * np + 1 : i_tilde * np)
1481  j_min = (i_tilde - 1) * np + 1
1482  x_min = xx(1, j_min)
1483 
1484  do j = (i_tilde - 1) * np + 2, i_tilde * np
1485  if (xx(1, j) < x_min) then
1486  x_min = xx(1, j)
1487  j_min = j
1488  end if
1489  end do
1490 
1491  j_next = mod(j_min, np) + 1 + (i_tilde - 1) * np
1492  j_previous = modulo(j_min - 2, np) + 1 + (i_tilde - 1) * np
1493 
1494 !-- search for extremum in interval [j_previous, j_min]
1495 
1496 !-- determine coefficients of quadratic equation to solve (xx'=0)
1497  n1 = j_previous
1498  n2 = n1 + 1
1499  c1 = xx(1, n1)
1500  c2 = xx(3, n1)
1501  c3 = xx(1, n2)
1502  c4 = xx(3, n2)
1503  aa = 3._r8 * (c1 + c2 - c3 + c4)
1504  bb = 2._r8 * (c4 - c2)
1505  cc = 3._r8 * (c3 - c1) - c2 - c4
1506  det = bb**2 - 4._r8 * aa * cc
1507  ifail = 0
1508  if (det > 0._r8) then ! has real root
1509  s_min = root(aa, bb, cc, det, 1._r8)
1510  if (abs(s_min) > 1._r8 + tol) then
1511  s_min = root(aa, bb, cc, det, -1._r8)
1512  if (abs(s_min) > 1._r8 + tol) then
1513  ifail = 1 ! no root inside interval s=[-1,1]
1514  end if
1515  end if
1516  else
1517  ifail = 1 ! no real root
1518  end if
1519 
1520  if (ifail == 1) then
1521 !-- search for extremum in interval [j_min, j_next]
1522 
1523  ifail = 0
1524 !-- determine coefficients of quadratic equation to solve (xx'=0)
1525  n1 = j_min
1526  n2 = n1 + 1
1527  c1 = xx(1, n1)
1528  c2 = xx(3, n1)
1529  c3 = xx(1, n2)
1530  c4 = xx(3, n2)
1531  aa = 3._r8 * (c1 + c2 - c3 + c4)
1532  bb = 2._r8 * (c4 - c2)
1533  cc = 3._r8 * (c3 - c1) - c2 - c4
1534  det = bb**2 - 4._r8 * aa * cc
1535  if (det > 0._r8) then ! has real root
1536  s_min = root(aa, bb, cc, det, 1._r8)
1537  if (abs(s_min) > 1._r8 + tol) then
1538  s_min = root(aa, bb, cc, det, -1._r8)
1539  if (abs(s_min) > 1._r8 + tol) then
1540  ifail = 1 ! no root inside interval s=[-1,1]
1541  end if
1542  end if
1543  else
1544  ifail = 1
1545  end if
1546  end if
1547 
1548  if (ifail == 1) then
1549  write(iu6, *) 'Warning: No root found for R_min'
1550  write(iu6, *) 'root replaced by node value'
1551  s_min = -1._r8
1552  end if
1553 
1554 !-- evaluate Hermite polynomial at s_min
1555  x_min = 0.25_r8 * (c1 * (s_min - 1)**2 * (s_min + 2) &
1556  + c2 * (s_min - 1)**2 * (s_min + 1) - c3 * (s_min + 1)**2 &
1557  * (s_min - 2) + c4 * (s_min + 1)**2 * (s_min - 1))
1558  end if
1559 
1560  end if
1561 
1562  f_r_min = rvac + rminor * x_min
1563 
1564  end function r_min
1565 
1566 !-----------------------------------------------------------------------
1567  function r_max(rminor, xaxis, i) result(f_R_max)
1568 !-----------------------------------------------------------------------
1569 ! function to determine the largest value of the major radius
1570 ! on a flux surface
1571 ! dependencies: xx
1572 !-----------------------------------------------------------------------
1573 
1574  use math_functions, only : root
1575  use mod_dat, only : rvac
1576  use mod_helena_io
1577 
1578  implicit none
1579 
1580  integer(itm_i4), intent(in) :: i
1581  real(r8), intent(in) :: rminor, xaxis
1582 
1583  real(r8) :: f_r_max
1584 
1585  integer(itm_i4) :: i_tilde
1586  integer(itm_i4) :: n1, n2
1587  integer(itm_i4) :: j, j_max
1588  integer(itm_i4) :: j_next, j_previous
1589  integer(itm_i4) :: ifail
1590  real(r8) :: x_max, s_max
1591  real(r8) :: c1, c2, c3, c4
1592  real(r8) :: aa, bb, cc, det ! annotation as in solve_cubic
1593  real(r8), parameter :: tol = 1.e-8_r8
1594 
1595  if (ias == 0) then ! up-down-symmetric case
1596 
1597  i_tilde = nr - i + 1
1598  j_max = (i_tilde - 1) * np + 1 ! first point in half circle = outboard side
1599  x_max = xx(1, j_max)
1600 
1601  else
1602 
1603  if (i == 1) then ! magnetic axis
1604  x_max = xaxis
1605  else
1606 
1607 !-- radial index runs from boundary to axis
1608  i_tilde = nr - i + 1
1609 
1610 !-- search for maximum in xx(1, (i_tilde - 1) * np + 1 : i_tilde * np)
1611  j_max = (i_tilde - 1) * np + 1
1612  x_max = xx(1, j_max)
1613 
1614  do j = (i_tilde - 1) * np + 2, i_tilde * np
1615  if (xx(1, j) > x_max) then
1616  x_max = xx(1, j)
1617  j_max = j
1618  end if
1619  end do
1620 
1621  j_next = mod(j_max, np) + 1 + (i_tilde - 1) * np
1622  j_previous = modulo(j_max - 2, np) + 1 + (i_tilde - 1) * np
1623 
1624 !-- search for extremum in interval [j_previous, j_max]
1625 
1626 !-- determine coefficients of quadratic equation to solve (xx'=0)
1627  n1 = j_previous
1628  n2 = n1 + 1
1629  c1 = xx(1, n1)
1630  c2 = xx(3, n1)
1631  c3 = xx(1, n2)
1632  c4 = xx(3, n2)
1633  aa = 3._r8 * (c1 + c2 - c3 + c4)
1634  bb = 2._r8 * (c4 - c2)
1635  cc = 3._r8 * (c3 - c1) - c2 - c4
1636  det = bb**2 - 4._r8 * aa * cc
1637  ifail = 0
1638  if (det > 0._r8) then ! has real root
1639  s_max = root(aa, bb, cc, det, 1._r8)
1640  if (abs(s_max) > 1._r8 + tol) then
1641  s_max = root(aa, bb, cc, det, -1._r8)
1642  if (abs(s_max) > 1._r8 + tol) then
1643  ifail = 1 ! no root inside interval s=[-1,1]
1644  end if
1645  end if
1646  else
1647  ifail = 1 ! no real root
1648  end if
1649 
1650  if (ifail == 1) then
1651 !-- search for extremum in interval [j_max, j_next]
1652 
1653  ifail = 0
1654 !-- determine coefficients of quadratic equation to solve (xx'=0)
1655  n1 = j_max
1656  n2 = n1 + 1
1657  c1 = xx(1, n1)
1658  c2 = xx(3, n1)
1659  c3 = xx(1, n2)
1660  c4 = xx(3, n2)
1661  aa = 3._r8 * (c1 + c2 - c3 + c4)
1662  bb = 2._r8 * (c4 - c2)
1663  cc = 3._r8 * (c3 - c1) - c2 - c4
1664  det = bb**2 - 4._r8 * aa * cc
1665  if (det > 0._r8) then ! has real root
1666  s_max = root(aa, bb, cc, det, 1._r8)
1667  if (abs(s_max) > 1._r8 + tol) then
1668  s_max = root(aa, bb, cc, det, -1._r8)
1669  if (abs(s_max) > 1._r8 + tol) then
1670  ifail = 1 ! no root inside interval s=[-1,1]
1671  end if
1672  end if
1673  else
1674  ifail = 1
1675  end if
1676  end if
1677 
1678  if (ifail == 1) then
1679  write(iu6, *) 'Warning: No root found for R_max'
1680  write(iu6, *) 'root replaced by node value'
1681  s_max = -1._r8
1682  end if
1683 
1684 !-- evaluate Hermite polynomial at s_max
1685  x_max = 0.25_r8 * (c1 * (s_max - 1)**2 * (s_max + 2) &
1686  + c2 * (s_max - 1)**2 * (s_max + 1) - c3 * (s_max + 1)**2 &
1687  * (s_max - 2) + c4 * (s_max + 1)**2 * (s_max - 1))
1688  end if
1689 
1690  end if
1691 
1692  f_r_max = rvac + rminor * x_max
1693 
1694  end function r_max
1695 
1696 !-----------------------------------------------------------------------
1697  function rz_top(rminor, xaxis, yaxis, i) result(f_RZ_top)
1698 !-----------------------------------------------------------------------
1699 ! function to determine the R,Z coordinates of the uppermost point
1700 ! of a flux surface
1701 ! dependencies: xx, yy
1702 !-----------------------------------------------------------------------
1703 
1704  use math_functions, only : root
1705  use mod_dat, only : rvac, zgeo
1706  use mod_helena_io
1707 
1708  implicit none
1709 
1710  integer(itm_i4), intent(in) :: i
1711  real(r8), intent(in) :: rminor, xaxis, yaxis
1712 
1713  type(rz_coordinate) :: f_rz_top
1714 
1715  integer(itm_i4) :: i_tilde
1716  integer(itm_i4) :: n1, n2
1717  integer(itm_i4) :: j, j_max
1718  integer(itm_i4) :: j_next, j_previous
1719  integer(itm_i4) :: ifail
1720  real(r8) :: x_max, y_max, s_max
1721  real(r8) :: cx1, cx2, cx3, cx4
1722  real(r8) :: cy1, cy2, cy3, cy4
1723  real(r8) :: aa, bb, cc, det ! annotation as in solve_cubic
1724  real(r8), parameter :: tol = 1.e-8_r8
1725 
1726  if (i == 1) then ! magnetic axis
1727  x_max = xaxis
1728  y_max = yaxis
1729  else
1730 
1731 !-- radial index runs from boundary to axis
1732  i_tilde = nr - i + 1
1733 
1734 !-- search for maximum in yy(1, (i_tilde - 1) * np + 1 : i_tilde * np)
1735  j_max = (i_tilde - 1) * np + 1
1736  y_max = yy(1, j_max)
1737 
1738  do j = (i_tilde - 1) * np + 2, i_tilde * np
1739  if (yy(1, j) > y_max) then
1740  y_max = yy(1, j)
1741  j_max = j
1742  end if
1743  end do
1744 
1745  j_next = mod(j_max, np) + 1 + (i_tilde - 1) * np
1746  j_previous = modulo(j_max - 2, np) + 1 + (i_tilde - 1) * np
1747 
1748 !-- search for extremum in interval [j_previous, j_max]
1749 
1750 !-- determine coefficients of quadratic equation to solve (yy'=0)
1751  n1 = j_previous
1752  n2 = n1 + 1
1753  cx1 = xx(1, n1)
1754  cx2 = xx(3, n1)
1755  cx3 = xx(1, n2)
1756  cx4 = xx(3, n2)
1757  cy1 = yy(1, n1)
1758  cy2 = yy(3, n1)
1759  cy3 = yy(1, n2)
1760  cy4 = yy(3, n2)
1761  aa = 3._r8 * (cy1 + cy2 - cy3 + cy4)
1762  bb = 2._r8 * (cy4 - cy2)
1763  cc = 3._r8 * (cy3 - cy1) - cy2 - cy4
1764  det = bb**2 - 4._r8 * aa * cc
1765  ifail = 0
1766  if (det > 0._r8) then ! has real root
1767  s_max = root(aa, bb, cc, det, 1._r8)
1768  if (abs(s_max) > 1._r8 + tol) then
1769  s_max = root(aa, bb, cc, det, -1._r8)
1770  if (abs(s_max) > 1._r8 + tol) then
1771  ifail = 1 ! no root inside interval s=[-1,1]
1772  end if
1773  end if
1774  else
1775  ifail = 1 ! no real root
1776  end if
1777 
1778  if (ifail == 1) then
1779 !-- search for extremum in interval [j_max, j_next]
1780 
1781  ifail = 0
1782 !-- determine coefficients of quadratic equation to solve (yy'=0)
1783  n1 = j_max
1784  n2 = n1 + 1
1785  cx1 = xx(1, n1)
1786  cx2 = xx(3, n1)
1787  cx3 = xx(1, n2)
1788  cx4 = xx(3, n2)
1789  cy1 = yy(1, n1)
1790  cy2 = yy(3, n1)
1791  cy3 = yy(1, n2)
1792  cy4 = yy(3, n2)
1793  aa = 3._r8 * (cy1 + cy2 - cy3 + cy4)
1794  bb = 2._r8 * (cy4 - cy2)
1795  cc = 3._r8 * (cy3 - cy1) - cy2 - cy4
1796  det = bb**2 - 4._r8 * aa * cc
1797  if (det > 0._r8) then ! has real root
1798  s_max = root(aa, bb, cc, det, 1._r8)
1799  if (abs(s_max) > 1._r8 + tol) then
1800  s_max = root(aa, bb, cc, det, -1._r8)
1801  if (abs(s_max) > 1._r8 + tol) then
1802  ifail = 1 ! no root inside interval s=[-1,1]
1803  end if
1804  end if
1805  else
1806  ifail = 1
1807  end if
1808  end if
1809 
1810  if (ifail == 1) then
1811  write(iu6, *) 'Warning: No root found for RZ_top'
1812  write(iu6, *) 'root replaced by node value'
1813  s_max = -1._r8
1814  end if
1815 
1816 !-- evaluate Hermite polynomial at s_max for xx and yy
1817  x_max = 0.25_r8 * (cx1 * (s_max - 1)**2 * (s_max + 2) &
1818  + cx2 * (s_max - 1)**2 * (s_max + 1) - cx3 * (s_max + 1)**2 &
1819  * (s_max - 2) + cx4 * (s_max + 1)**2 * (s_max - 1))
1820  y_max = 0.25_r8 * (cy1 * (s_max - 1)**2 * (s_max + 2) &
1821  + cy2 * (s_max - 1)**2 * (s_max + 1) - cy3 * (s_max + 1)**2 &
1822  * (s_max - 2) + cy4 * (s_max + 1)**2 * (s_max - 1))
1823  end if
1824 
1825  f_rz_top%R = rvac + rminor * x_max
1826  f_rz_top%Z = zgeo + rminor * y_max
1827 
1828  end function rz_top
1829 
1830 !-----------------------------------------------------------------------
1831  function rz_bottom(rminor, xaxis, yaxis, i) result(f_RZ_bottom)
1832 !-----------------------------------------------------------------------
1833 ! function to determine the R,Z coordinates of the lowermost point
1834 ! of a flux surface
1835 ! dependencies: xx, yy
1836 !-----------------------------------------------------------------------
1837 
1838  use math_functions, only : root
1839  use mod_dat, only : rvac, zgeo
1840  use mod_helena_io
1841 
1842  implicit none
1843 
1844  integer(itm_i4), intent(in) :: i
1845  real(r8), intent(in) :: rminor, xaxis, yaxis
1846 
1847  type(rz_coordinate) :: f_rz_bottom
1848 
1849  integer(itm_i4) :: i_tilde
1850  integer(itm_i4) :: n1, n2
1851  integer(itm_i4) :: j, j_min
1852  integer(itm_i4) :: j_next, j_previous
1853  integer(itm_i4) :: ifail
1854  real(r8) :: x_min, y_min, s_min
1855  real(r8) :: cx1, cx2, cx3, cx4
1856  real(r8) :: cy1, cy2, cy3, cy4
1857  real(r8) :: aa, bb, cc, det ! annotation as in solve_cubic
1858  real(r8), parameter :: tol = 1.e-8_r8
1859 
1860 
1861  if (ias == 0) then ! up-down-symmetric case (top -> bottom)
1862 
1863  f_rz_bottom = rz_top(rminor, xaxis, yaxis, i)
1864  f_rz_bottom%Z = -f_rz_bottom%Z
1865 
1866  else
1867 
1868  if (i == 1) then ! magnetic axis
1869  x_min = xaxis
1870  y_min = yaxis
1871  else
1872 
1873 !-- radial index runs from boundary to axis
1874  i_tilde = nr - i + 1
1875 
1876 !-- search for minimum in yy(1, (i_tilde - 1) * np + 1 : i_tilde * np)
1877  j_min = (i_tilde - 1) * np + 1
1878  y_min = yy(1, j_min)
1879 
1880  do j = (i_tilde - 1) * np + 2, i_tilde * np
1881  if (yy(1, j) < y_min) then
1882  y_min = yy(1, j)
1883  j_min = j
1884  end if
1885  end do
1886 
1887  j_next = mod(j_min, np) + 1 + (i_tilde - 1) * np
1888  j_previous = modulo(j_min - 2, np) + 1 + (i_tilde - 1) * np
1889 
1890 !-- search for extremum in interval [j_previous, j_min]
1891 
1892 !-- determine coefficients of quadratic equation to solve (yy'=0)
1893  n1 = j_previous
1894  n2 = n1 + 1
1895  cx1 = xx(1, n1)
1896  cx2 = xx(3, n1)
1897  cx3 = xx(1, n2)
1898  cx4 = xx(3, n2)
1899  cy1 = yy(1, n1)
1900  cy2 = yy(3, n1)
1901  cy3 = yy(1, n2)
1902  cy4 = yy(3, n2)
1903  aa = 3._r8 * (cy1 + cy2 - cy3 + cy4)
1904  bb = 2._r8 * (cy4 - cy2)
1905  cc = 3._r8 * (cy3 - cy1) - cy2 - cy4
1906  det = bb**2 - 4._r8 * aa * cc
1907  ifail = 0
1908  if (det > 0._r8) then ! has real root
1909  s_min = root(aa, bb, cc, det, 1._r8)
1910  if (abs(s_min) > 1._r8 + tol) then
1911  s_min = root(aa, bb, cc, det, -1._r8)
1912  if (abs(s_min) > 1._r8 + tol) then
1913  ifail = 1 ! no root inside interval s=[-1,1]
1914  end if
1915  end if
1916  else
1917  ifail = 1 ! no real root
1918  end if
1919 
1920  if (ifail == 1) then
1921 !-- search for extremum in interval [j_min, j_next]
1922 
1923  ifail = 0
1924 !-- determine coefficients of quadratic equation to solve (yy'=0)
1925  n1 = j_min
1926  n2 = n1 + 1
1927  cx1 = xx(1, n1)
1928  cx2 = xx(3, n1)
1929  cx3 = xx(1, n2)
1930  cx4 = xx(3, n2)
1931  cy1 = yy(1, n1)
1932  cy2 = yy(3, n1)
1933  cy3 = yy(1, n2)
1934  cy4 = yy(3, n2)
1935  aa = 3._r8 * (cy1 + cy2 - cy3 + cy4)
1936  bb = 2._r8 * (cy4 - cy2)
1937  cc = 3._r8 * (cy3 - cy1) - cy2 - cy4
1938  det = bb**2 - 4._r8 * aa * cc
1939  if (det > 0._r8) then ! has real root
1940  s_min = root(aa, bb, cc, det, 1._r8)
1941  if (abs(s_min) > 1._r8 + tol) then
1942  s_min = root(aa, bb, cc, det, -1._r8)
1943  if (abs(s_min) > 1._r8 + tol) then
1944  ifail = 1 ! no root inside interval s=[-1,1]
1945  end if
1946  end if
1947  else
1948  ifail = 1
1949  end if
1950  end if
1951 
1952  if (ifail == 1) then
1953  write(iu6, *) 'Warning: No root found for RZ_bottom'
1954  write(iu6, *) 'root replaced by node value'
1955  s_min = -1._r8
1956  end if
1957 
1958 !-- evaluate Hermite polynomial at s_min for xx and yy
1959  x_min = 0.25_r8 * (cx1 * (s_min - 1)**2 * (s_min + 2) &
1960  + cx2 * (s_min - 1)**2 * (s_min + 1) - cx3 * (s_min + 1)**2 &
1961  * (s_min - 2) + cx4 * (s_min + 1)**2 * (s_min - 1))
1962  y_min = 0.25_r8 * (cy1 * (s_min - 1)**2 * (s_min + 2) &
1963  + cy2 * (s_min - 1)**2 * (s_min + 1) - cy3 * (s_min + 1)**2 &
1964  * (s_min - 2) + cy4 * (s_min + 1)**2 * (s_min - 1))
1965  end if
1966 
1967  f_rz_bottom%R = rvac + rminor * x_min
1968  f_rz_bottom%Z = zgeo + rminor * y_min
1969 
1970  end if
1971 
1972  end function rz_bottom
1973 
1974 
1975 !---------------------------------------------------------------------------
1976  subroutine current_at_nodes(a, equilibrium)
1977 !---------------------------------------------------------------------------
1978 ! This subroutine calculates the toroidal and parallel current densities
1979 ! j_phi and j_par at the nodes of the (psi, theta) grid.
1980 ! j_phi = j_tor
1981 ! j_par = j . B / B_0
1982 ! It requires a prior call to mapping because various fields in the
1983 ! equilibrium CPO and cchi have to be filled.
1984 !---------------------------------------------------------------------------
1985  use mod_dat
1986  use mod_helena_io, only : iu6, out_he
1987  use mod_map
1988  use mod_nodes, only : radpsi
1989  use mod_output
1990  use mod_parameter
1991  use phys_functions, only : rh_side, b2
1992 
1993  use euitm_schemas
1994 
1995  implicit none
1996 
1997  real(r8), intent(in) :: a
1998  type(type_equilibrium) :: equilibrium
1999 
2000  real(r8) :: j0, f_dia
2001  real(r8) :: b_field
2002  real(r8) :: a0, a1, a2, a3
2003  real(r8) :: s, s2, s3
2004  real(r8) :: zpsi, dzpsi, ddzpsi, zchi
2005  integer(itm_i4) :: nelm
2006  integer(itm_i4) :: n1, n2, n3, n4
2007  integer(itm_i4) :: no, nom, nop
2008  integer(itm_i4) :: jbase
2009  integer(itm_i4) :: ifail
2010  integer(itm_i4) :: i, j, k
2011  logical chin
2012 
2013 !-- check allocation of required quantities
2014  if (.not. associated(equilibrium%profiles_1d%psi) &
2015  .or. .not. associated(equilibrium%profiles_1d%F_dia)) then
2016  write(iu6, *) 'ERROR: psi or F_dia not defined'
2017  stop
2018  end if
2019 
2020 !-- normalization constants
2021  j0 = bvac / (mu0 * eps * rvac * alfa)
2022 
2023 !-- determine positions of equidistant chi's
2024 ! and calculate matrix elements
2025  if (ias == 0) then
2026  do j = 1, nchi
2027  chikn(j) = pi * (j - 1) / dble(nchi - 1)
2028  chi(j) = pi * (j - 1) / dble(nchi - 1)
2029  end do
2030  else
2031  do j = 1, nchi
2032  chikn(j) = 2._r8 * pi * (j - 1) / dble(nchi)
2033  chi(j) = 2._r8 * pi * (j - 1) / dble(nchi)
2034  end do
2035  end if
2036 
2037  do i = 1, nr - 1
2038  call radial_mesh(radpsi(i), zpsi, dzpsi, ddzpsi)
2039 
2040 !-- first point is known
2041  no = (i - 1) * nchi + 1
2042  k = 1
2043  s = -1._r8
2044 
2045  n1 = (i - 1) * np + k
2046  n2 = n1 + 1
2047  n3 = n2 + np
2048  n4 = n1 + np
2049 
2050  equilibrium%profiles_2d(1)%jphi(nr - i + 1, 1) = - j0 * eps &
2051  * rh_side(a, xx(1, n1), 0._r8, 0._r8, 0._r8, 0._r8, &
2052  psi(4 * (n1 - 1) + 1), 0._r8, 0._r8)
2053  f_dia = equilibrium%profiles_1d%F_dia(nr - i + 1) &
2054  * equilibrium%eqgeometry%a_minor / equilibrium%profiles_1d%psi(nr)
2055  b_field = sqrt(b2(a, xx(1, n1), 0._r8, 0._r8, 0._r8, 0._r8, &
2056  psi(4 * (n1 - 1) + 1), 0._r8, f_dia) &
2057  / (equilibrium%eqgeometry%a_minor * equilibrium%eqgeometry%geom_axis%r &
2058  / equilibrium%profiles_1d%psi(nr))**2)
2059  equilibrium%profiles_2d(1)%jpar(nr - i + 1, 1) = &
2060  (equilibrium%profiles_1d%F_dia(nr - i + 1) &
2061  / equilibrium%coord_sys%position%r(nr - i + 1, 1) &
2062  * equilibrium%profiles_2d(1)%jphi(nr - i + 1, 1) &
2063  + equilibrium%profiles_1d%ffprime(nr - i + 1) &
2064  / equilibrium%profiles_1d%F_dia(nr - i + 1) / (twopi**2 * mu0 &
2065  * equilibrium%coord_sys%position%r(nr - i + 1, 1)**2) &
2066  * equilibrium%coord_sys%g_11(nr - i + 1, 1)) / b_field
2067 
2068  jbase = 2
2069 
2070  do k = 1, np - 1
2071 
2072  chin = .false.
2073 
2074  do j = jbase, nchi
2075 
2076  zchi = chikn(j)
2077  no = (i - 1) * nchi + j
2078 
2079  if ((cchi(1, (i - 1) * np + k) <= zchi .and. cchi(1, &
2080  (i - 1) * np + k + 1) >= zchi) .or. (j == nchi &
2081  .and. k == np - 1 .and. ias == 0)) then
2082 
2083  chin = .true.
2084 
2085  nom = (i - 1) * np + k
2086  nop = nom + 1
2087 
2088  a3 = (cchi(1, nom) + cchi(3, nom) - cchi(1, nop) &
2089  + cchi(3, nop)) / 4._r8
2090  a2 = (-cchi(3, nom) + cchi(3, nop)) / 4._r8
2091  a1 = (-3 * cchi(1, nom) - cchi(3, nom) + 3 * cchi(1, nop) &
2092  - cchi(3, nop)) / 4._r8
2093  a0 = (2 * cchi(1, nom) + cchi(3, nom) + 2 * cchi(1, nop) &
2094  - cchi(3, nop)) / 4._r8 - zchi
2095 
2096  call solve_cubic(a0, a1, a2, a3, s, s2, s3, ifail)
2097 
2098  if (ias == 0 .and. j == nchi) then
2099  s = 1._r8
2100  ifail = 0
2101  end if
2102 
2103  if (ifail == 0) then
2104  n1 = (i - 1) * np + k
2105  n2 = n1 + 1
2106  n3 = n2 + np
2107  n4 = n1 + np
2108 
2109  equilibrium%profiles_2d(1)%jphi(nr - i + 1, j) = - j0 * eps &
2110  * rh_side(a, xx(1, n1), 0._r8, 0._r8, 0._r8, 0._r8, &
2111  psi(4 * (n1 - 1) + 1), 0._r8, 0._r8)
2112  f_dia = equilibrium%profiles_1d%F_dia(nr - i + 1) &
2113  * equilibrium%eqgeometry%a_minor &
2114  / equilibrium%profiles_1d%psi(nr)
2115  b_field = sqrt(b2(a, xx(1, n1), 0._r8, 0._r8, 0._r8, 0._r8, &
2116  psi(4 * (n1 - 1) + 1), 0._r8, f_dia) &
2117  / (equilibrium%eqgeometry%a_minor &
2118  * equilibrium%eqgeometry%geom_axis%r &
2119  / equilibrium%profiles_1d%psi(nr))**2)
2120  equilibrium%profiles_2d(1)%jpar(nr - i + 1, j) = &
2121  (equilibrium%profiles_1d%F_dia(nr - i + 1) &
2122  / equilibrium%coord_sys%position%r(nr - i + 1, j) &
2123  * equilibrium%profiles_2d(1)%jphi(nr - i + 1, j) &
2124  + equilibrium%profiles_1d%ffprime(nr - i + 1) * twopi &
2125  / equilibrium%profiles_1d%F_dia(nr - i + 1) / (twopi**2 * mu0 &
2126  * equilibrium%coord_sys%position%r(nr - i + 1, j)**2) &
2127  * equilibrium%coord_sys%g_11(nr - i + 1, j)) / b_field
2128 
2129  else
2130  if (standard_output) &
2131  write(out_he, *) 'error in solve_cubic i,j,k : ', i, j, k, &
2132  s, s2, s3
2133  if (verbosity > 3) then
2134  write(iu6, *) a0, a1, a2, a3, zchi
2135  write(iu6, *) cchi(1, (i - 1) * np + k), cchi(1, (i - 1) &
2136  * np + k + 1)
2137  end if
2138  end if
2139 
2140  else if (chin) then
2141  jbase = j
2142  goto 70
2143  end if
2144  end do
2145  70 end do
2146  end do
2147 
2148 !-- cchi not needed anymore
2149  deallocate(cchi)
2150 
2151  end subroutine current_at_nodes
2152 
2153 end module mapping_mod
write_old_equilibrium
subroutine write_old_equilibrium(path, equilibrium_out)
Definition: write_old_equilibrium.f90:1
mapping_mod::mapping
subroutine mapping(cx, cy, xaxis, yaxis, a, equilibrium_out)
Definition: mapping_mod.f90:16
mod_dat
Definition: mod_dat.f90:1
mapping_mod::r_min
real(r8) function r_min(rminor, xaxis, i)
Definition: mapping_mod.f90:1436
radial_mesh
subroutine radial_mesh(rpsi, zpsi, dzpsi, ddzpsi)
Definition: radial_mesh.f90:1
plot_data::plot_data_2d
subroutine plot_data_2d(type_plot, x, y, printname, z)
Definition: plot_data.f90:109
mod_equilibrium
Definition: mod_equilibrium.f90:1
profiles
subroutine profiles(p0, rbphi, dp0, drbphi, a)
Definition: profiles.f90:1
plot_data
Definition: plot_data.f90:1
mod_mesh
Definition: mod_mesh.f90:1
phys_functions::b_fun
real(r8) function b_fun(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:265
mapping_mod::rz_top
type(rz_coordinate) function rz_top(rminor, xaxis, yaxis, i)
Definition: mapping_mod.f90:1697
mapping_mod::rz_bottom
type(rz_coordinate) function rz_bottom(rminor, xaxis, yaxis, i)
Definition: mapping_mod.f90:1831
math_functions::integrate
real(r8) function, dimension(n), public integrate(f_spline, x, n, inv)
Definition: math_functions.f90:17
mod_interpolation::interpolation
subroutine interpolation(type_interp, xn1, xn2, xn3, xn4, r, s, x, xr, xs, xrs, xrr, xss, yn1, yn2, yn3, yn4, pn1, pn2, pn3, pn4, yr, ys, ps)
Definition: mod_interpolation.f90:10
output
subroutine output(NGRID, betpol, zli3)
Definition: Eq2_m.f:1
mod_profiles::allocate_spline_coefficients
subroutine allocate_spline_coefficients(spline, n)
Definition: mod_profiles.f90:27
mod_output
Definition: mod_output.f90:1
mod_interpolation
Definition: mod_interpolation.f90:1
phys_functions::gradpsi2
real(r8) function gradpsi2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:160
mod_equilibrium::fill_equilibrium
subroutine, public fill_equilibrium(a, xaxis, yaxis, current_averaging, equilibrium)
Definition: mod_equilibrium.f90:256
phys_functions::one_over_r2
real(r8) function one_over_r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:140
phys_functions::gradpsi
real(r8) function gradpsi(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:196
spline
subroutine spline(N, X, Y, ALFA, BETA, TYP, A, B, C, D)
Definition: solution6.f90:655
spwert
REAL *8 function spwert(N, XWERT, A, B, C, D, X, ABLTG)
Definition: solution6.f90:155
phys_functions::rh_side
real(r8) function rh_side(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:30
phys_functions::one_over_b2
real(r8) function one_over_b2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:313
mod_profiles
Definition: mod_profiles.f90:1
dgamma_dpsi
real(r8) function dgamma_dpsi(flux)
Definition: dgamma_dpsi.f90:1
r_inboard
subroutine r_inboard
Definition: mapping_mod.f90:1397
math_functions::flux_surface_average
real(r8) function, public flux_surface_average(i, xaxis, a, F_dia, type, func)
Definition: math_functions.f90:436
phys_functions::b2
real(r8) function b2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:289
plot_data::plot_data_1d
subroutine plot_data_1d(type_plot, x, y, np, printname, z)
Definition: plot_data.f90:11
mod_nodes
Definition: mod_nodes.f90:1
math_functions
Definition: math_functions.f90:1
mod_parameter
Definition: mod_parameter.f90:1
phys_functions::gradrho2_over_b2
real(r8) function gradrho2_over_b2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:336
solve_cubic
subroutine solve_cubic(c0, c1, c2, c3, x1, x2, x3, ifail)
Definition: solve_cubic.f90:1
mod_profiles::deallocate_spline_coefficients
subroutine deallocate_spline_coefficients(spline)
Definition: mod_profiles.f90:43
mapping_mod::rz_coordinate
Definition: mapping_mod.f90:8
math_functions::root
real(r8) function, public root(a, b, c, d, sgn)
Definition: math_functions.f90:203
phys_functions::one_over_r
real(r8) function one_over_r(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:61
mapping_mod
Definition: mapping_mod.f90:1
phys_functions::r_major
real(r8) function r_major(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:102
mod_helena_io
Definition: mod_helena_io.f90:1
dp_dpsi
real(r8) function dp_dpsi(flux)
Definition: dp_dpsi.f90:1
mapping_mod::r_max
real(r8) function r_max(rminor, xaxis, i)
Definition: mapping_mod.f90:1567
mapping_mod::current_at_nodes
subroutine current_at_nodes(a, equilibrium)
Definition: mapping_mod.f90:1976
mod_suydam
Definition: mod_suydam.f90:1
mod_gauss
Definition: mod_gauss.f90:1
mod_ne_te_profile
Definition: mod_ne_te_profile.f90:1
helena_spline::spline_coefficients
Definition: helena_spline.f90:7
phys_functions
Definition: phys_functions.f90:1
mod_ne_te_profile::ne_te_profile
subroutine ne_te_profile(dt, psi, f, df, const, df_in)
Definition: mod_ne_te_profile.f90:12
mod_dete
Definition: mod_dete.f90:1
phys_functions::gradpsi2_over_r2
real(r8) function gradpsi2_over_r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:218
mod_map
Definition: mod_map.f90:1
pplot
subroutine pplot(MX, MY, X, Y, NPTS, INC)
Definition: ppplib.f:599
helena_spline
Definition: helena_spline.f90:1