V3FIT
profile_functions.f
1 ! File profile_functions.f contains
2 ! FUNCTION pcurr
3 ! FUNCTION piota
4 ! FUNCTION pmass
5 ! (JDH 2010-03-30)
6 !******************************************************************************
7 
8  FUNCTION pcurr (xx)
9 ! Function to compute the current profile
10 
11 ! Variables declared in vmec_input:
12 ! ac array (0:20) of coefficients
13 ! ac_aux_s Auxiliary array s, s-values used for splines
14 ! ac_aux_f Auxiliary array f, function-values used for splines
15 ! bloat used to expand the current profile
16 ! pcurr_type character, specifies the parameterization of the profile
17 ! | - X for parametrization of I-prime(s), _ for I(s)
18 ! sum_cossq_s X Sum of cos**2-waves - I-prime
19 ! sum_cossq_sqrts X Sum of cos**2-waves with respect to sqrt(s) - I-prime
20 ! sum_cossq_s_free X Sum of cos**2-waves up to 7, free position and width - I-prime
21 ! gauss_trunc X Truncated Gaussian - I-prime
22 ! two_power X Two powers - ac(0) * (1 - s ** ac(1)) ** ac(2)
23 ! two_power_gs X Two powers with gaussian peaks -
24 ! ac(0) * ((1 - s ** ac(1)) ** ac(2))*(1 + Sum[ac(i)*Exp(-(s - ac(i+1))/ac(i+2)) ** 2])
25 ! sum_atan _ sum of arctangents
26 ! power_series_I _ Power series for I(s) (NOT default)
27 ! Akima_spline_Ip X Akima spline for I-prime(s)
28 ! Akima_spline_I _ Akima spline for I(s)
29 ! cubic_spline_Ip X cubic spline for I-prime(s)
30 ! cubic_spline_I _ cubic spline for I(s)
31 ! pedestal _ Pedestal profile
32 ! rational _ Rational function (ratio of polynomials)
33 ! line_segment_Ip X Line segments for I-prime(s)
34 ! line_segment_I _ Line segments for I(s)
35 ! power_series X Power Series for I-prime(s) (Default)
36 
37 ! Local Variables
38 ! i integer counter
39 ! ioff offset for ac array. Should be zero
40 ! iflag error flag for spline call
41 ! xx real argument
42 ! x constrained to be between 0 and 1
43 ! xp variable for Gauss_legendre quadrature
44 ! pcurr_type_lc character, pcurr_type -> lower case
45 ! gli index for Gauss-Legendre quadrature loop
46 ! gln order of Gauss-Legendre quadrature
47 ! glx array of abscissa values for Gauss-Legendre quadrature
48 ! glw array of wieghts for Gauss-Legendre quadrature
49 
50 ! Note that the profile that is parameterized is often I-prime, whereas
51 ! I(x) (= Integral_from_0_to_x I-prime(s) ds) is the function that pcurr
52 ! returns. For the default case of a power series, the integral can be
53 ! computed analytically. For other cases, a numerical quadrature is done,
54 ! using a 10-point Gauss-Legendre quadrature.
55  USE stel_kinds
56  USE stel_constants, ONLY: zero, one, pi
57  USE vmec_input, ONLY: ac, bloat, pcurr_type, ac_aux_s, ac_aux_f
58  USE line_segment
59  USE functions
60  IMPLICIT NONE
61 ! ac assumed to be dimensioned (0:n), with n >= 20
62 !-----------------------------------------------
63  INTEGER :: i, ioff, iflag
64  REAL(rprec) :: xx, pcurr, x, xp, temp_num, temp_denom
65  CHARACTER(len=20) :: pcurr_type_lc
66 
67  INTEGER, PARAMETER :: gln = 10
68  INTEGER :: gli
69  REAL(rprec), DIMENSION(gln), PARAMETER :: glx = (/ &
70  & 0.01304673574141414, 0.06746831665550774, 0.1602952158504878, &
71  & 0.2833023029353764, 0.4255628305091844, 0.5744371694908156, &
72  & 0.7166976970646236, 0.8397047841495122, 0.9325316833444923, &
73  & 0.9869532642585859 /)
74  REAL(rprec), DIMENSION(gln), PARAMETER :: glw = (/ &
75  & 0.03333567215434407, 0.0747256745752903, 0.1095431812579910, &
76  & 0.1346333596549982, 0.1477621123573764, 0.1477621123573764, &
77  & 0.1346333596549982, 0.1095431812579910, 0.0747256745752903, &
78  & 0.03333567215434407 /)
79  REAL(rprec) :: g1,g2,g3,g4,a8,a12
80  REAL(rprec), dimension(21) :: xi, bsta, bend, wd
81  REAL(rprec) :: sqx,delx,delxsq,pisq
82  INTEGER :: ni
83  INTEGER :: ncssq
84 
85 !-----------------------------------------------
86 !
87 ! NOTE: AC COEFFICIENTS OBTAINED IN THREED1 FILE
88 ! BY MATCHING TO <JTOR> * dV/dPHI ~ SUM[x^(i+1) * ac(i)/(i+1)], i=0...UBOUND(ac)
89 !
90 ! Start of executable code
91 
92  x = min(abs(xx * bloat), one)
93  ioff = lbound(ac,1) ! Expected to be zero.
94 
95  pcurr = 0
96 
97 ! Convert to Lower Case, to avoid typo issues
98  pcurr_type_lc = pcurr_type
99  CALL tolower(pcurr_type_lc)
100  SELECT CASE(trim(pcurr_type_lc))
101 
102  CASE ('sum_cossq_s')
103 ! 20180218, Joachim Geiger
104 ! The idea was to use the combination of cos**2-waves located at
105 ! different radial locations to build up a current density in a
106 ! hopefully quite flexible way. Since the current is needed, the
107 ! integration is done analytically.
108 ! Sum of ncssq cos**2 terms put at different radial locations
109 ! and with a windowing around the maxima as current density:
110 ! ac(0) holds ncssq
111 ! ac(1)*(H(x)-H(x-dx))*(cos(pi*(x-xi(1))/(2*dx))**2
112 ! +ac(2)*(H(x-(xi(2)-dx))-H(x-(xi(2)+dx)))*cos(pi*(x-xi(2))/(2*dx))**2
113 ! + ...
114 ! +ac(ncssq)*(H(x-(xi(ncssq)-dx))-H(x-(xi(ncssq))))*cos(pi*(x-xi(ncssq))/(2*dx))**2
115 ! windowing is done with the Heavyside-functions producing the
116 ! Boxcar Function: H(x-a)-H(x-b) is, for a<b, 1 for a<x<b and zero otherwise.
117 ! use of variables:
118 ! half of interval width : dx=1/(ncssq-1.0) .
119 ! location of ith: xi=(i-1)*dx = (i-1)/(ncssq-1)
120 ! example: ncssq=5, dx=1/4.0=0.25, locations:0.0, 0.25, 0.5, ...
121 ! interval for first wave: 0.00 to 0.25, maximum at 0.00
122 ! interval for second wave: 0.00 to 0.50, maximum at 0.25
123 ! interval for third wave: 0.25 to 0.75, maximum at 0.50
124 ! interval for fourth wave: 0.50 to 1.00, maximum at 0.75
125 ! interval for fifth wave: 0.75 to 1.00, maximum at 1.00
126 ! Analytic integration results in a sum of sine-functions with
127 ! linear increases:
128 ! int(x(i)-dy to x) cos(pi*(t-xi(i))/(2*dx))**2 dx =
129 ! =0.5*(x-(x(i)-dx)+(dx/pi)*sin(pi(x-x(i))/dx))
130 ! For the very first part integration start at mid-interval,
131 ! i.e. x=0, x(i)=0:
132 ! int(0 to x) cos(pi*(t)/(2*dx))**2 dx =
133 ! =0.5*(x+(dx/pi)*sin(pi(x)/dx))
134 ! When adding up fully integrated intervals the first needs
135 ! also to be accounted for differently.
136 !
137  ncssq=int(ac(0))
138  if (ncssq .lt. 1 .or. ncssq .gt. 20) then
139  write(6,*) "Number of coeff.s for sum_cos_sq :",
140  & "1<= sum_cos_sq <= 21!"
141  stop 'Check input!'
142  endif
143  delx=1.0_dp/(ncssq-1.0_dp)
144  bsta = 2.0_dp
145  bend = 2.0_dp
146  do i=1,ncssq
147  xi(i)=(i-1)*delx ! location of maximum curr. dens.
148  bsta(i)=max(0.0_dp,xi(i)-delx) ! interval start of curr.dens.part
149  bend(i)=min(xi(i)+delx,1.0_dp) ! interval end of curr.dens.part
150  enddo
151  ni=1 ! first interval initially.
152  do i=1,ncssq ! check the interval of x in ( bsta(i) : bend(i) ]
153  if( x .le. bend(i) .and. x .gt. bsta(i))then
154  ni=i
155  endif
156  enddo
157  pcurr=0.0
158  do i=1,ni
159  if( x .gt. bsta(i) .and. x .le. bend(i) ) then
160  if(i .eq. 1) then ! the first has no linear term.
161  pcurr=pcurr + 0.5*ac(i)*(x+ &
162  & delx*sin(pi*(x-xi(i))/delx)/pi)
163  else
164  pcurr=pcurr + 0.5*ac(i)*(x-xi(i)+delx+ &
165  & delx*sin(pi*(x-xi(i))/delx)/pi)
166  endif
167  else
168  if(i .eq. 1) then !the first is only a half interval
169  pcurr=pcurr+0.5*ac(i)*delx
170  else
171  pcurr=pcurr+ac(i)*delx
172  endif
173  endif
174  enddo
175 
176  CASE ('sum_cossq_sqrts')
177 ! 20180219, Joachim Geiger
178 ! This is the implementation for using sqrt(s) as coordinate.
179 ! The same idea as for sum_cos_sq.
180 ! The analytical integration give different integrals.
181 ! Sum of ncssq cos**2 terms put at different radial locations
182 ! and with a windowing around the maxima as current density:
183 ! ac(0) holds ncssq
184 ! ac(1)*(H(x)-H(x-dx))*(cos(pi*(x-xi(1))/(2*dx))**2
185 ! +ac(2)*(H(x-(xi(2)-dx))-H(x-(xi(2)+dx)))*cos(pi*(x-xi(2))/(2*dx))**2
186 ! + ...
187 ! +ac(ncssq)*(H(x-(xi(ncssq)-dx))-H(x-(xi(ncssq))))*cos(pi*(x-xi(ncssq))/(2*dx))**2
188 ! windowing is done with the Heavyside-functions producing the
189 ! Boxcar Function: H(x-a)-H(x-b) is, for a<b, 1 for a<x<b and zero otherwise.
190 ! use of variables:
191 ! half of interval width : dx=1/(ncssq-1.0) .
192 ! location of ith: xi=(i-1)*dx = (i-1)/(ncssq-1)
193 ! example: ncssq=5, dx=1/4.0=0.25, locations:0.0, 0.25, 0.5, ...
194 ! Note the locations are with respect to sqrt(s), the s-value is given in brackets
195 ! interval for first wave: 0.00 to 0.25, maximum at 0.00 (0.0)
196 ! interval for second wave: 0.00 to 0.50, maximum at 0.25 (0.0625)
197 ! interval for third wave: 0.25 to 0.75, maximum at 0.50 (0.25)
198 ! interval for fourth wave: 0.50 to 1.00, maximum at 0.75 (0.5625)
199 ! interval for fifth wave: 0.75 to 1.00, maximum at 1.00 (1.0)
200 ! Analytic integration now wrt sqrt(s)=x results in a more complicated sum of
201 ! sine-functions and cosine-functions with and a linear increase:
202 ! int(x(i)-dy to x) cos(pi*(t-xi(i))/(2*dx))**2 x dx =
203 ! = 0.25*(x-(x(i))+dx**2/(2*pi**2)+
204 ! + (dx/(2*pi))*x*sin(pi(x-x(i))/dx)) +
205 ! + (dx**2/(2*pi**2))*x*cos(pi(x-x(i))/dx))
206 ! For the very first part integration start at mid-interval,
207 ! i.e. x=0, x(i)=0:
208 ! int(0 to x) cos(pi*(t)/(2*dx))**2 x dx =
209 ! =0.25*x+(dx/(2*pi)) * x * sin(pi(x)/dx)
210 ! +(dx**2/(2*pi**2))*(cos(pi(x)/dx)-1)
211 ! When adding up fully integrated intervals the first needs
212 ! also to be accounted for differently.
213 !
214  ncssq=int(ac(0))
215  if (ncssq .lt. 1 .or. ncssq .gt. 20) then
216  write(6,*) "Number of coeff.s for sum_cos_sq :",
217  & "1<= sum_cos_sq <= 21!"
218  stop 'Check input!'
219  endif
220  sqx=sqrt(x)
221  delx=1.0_dp/(ncssq-1.0_dp)
222  delxsq=delx*delx
223  pisq=pi*pi
224  bsta = 2.0_dp
225  bend = 2.0_dp
226  do i=1,ncssq
227  xi(i)=(i-1)*delx ! location of maximum curr. dens.
228  bsta(i)=max(0.0_dp,xi(i)-delx) ! interval start of curr.dens.part
229  bend(i)=min(xi(i)+delx,1.0_dp) ! interval end of curr.dens.part
230  enddo
231  ni=1 ! first interval initially.
232  do i=1,ncssq ! check the interval of x in ( bsta(i) : bend(i) ]
233  if( sqx .le. bend(i) .and. sqx .gt. bsta(i))then
234  ni=i
235  endif
236  enddo
237  pcurr=0.0
238  do i=1,ni
239  if( sqx .gt. bsta(i) .and. sqx .le. bend(i) ) then
240  if(i .eq. 1) then ! the first has no linear term.
241  pcurr=pcurr + ac(i)*(0.25_dp*x+ &
242  & (delx/(2*pi)*sqx*sin(pi*sqx/delx))+ &
243  & delxsq/(2*pisq)*(cos(pi*sqx/delx)-1.0_dp))
244  else
245  pcurr=pcurr + ac(i)*(0.25_dp*x+ &
246  & delx/(2*pi)*sqx*sin(pi*(sqx-xi(i))/delx)+ &
247  & delxsq/(2*pisq)*cos(pi*(sqx-xi(i))/delx)- &
248  & 0.25_dp*bsta(i)**2 + delxsq/(2*pisq))
249  endif
250  else
251  if(i .eq. 1) then !the first is only a half interval
252 ! if(x .gt. 0.0_dp )
253 ! & pcurr=pcurr+ac(i)*(0.25_dp-1/pisq)*delxsq
254  pcurr=pcurr+ac(i)*(0.25_dp-1/pisq)*delxsq
255  else
256  pcurr=pcurr+ac(i)*delx*xi(i)
257  endif
258  endif
259  enddo
260 
261  CASE ('sum_cossq_s_free')
262 ! 20180628, Joachim Geiger
263 ! The idea was to use the combination of cos**2-waves located at
264 ! different radial locations to build up a current density in a
265 ! hopefully quite flexible way. Since the current is needed, the
266 ! integration is done analytically.
267 ! This implementation allows a free positioning of up to seven
268 ! waves with individual periodicity length. This means broad and
269 ! small ones can be combined. The waves have not to fit fully or
270 ! with one half into the s-range from 0 to 1. This is the reason
271 ! that no version with sqrt(s) is used.
272 ! Every three coefficients hold the amplitude, the max-location
273 ! and the width, i.e. min. to min. distance:
274 ! AC(0+3*i) = amplitude
275 ! AC(1+3*i) = location
276 ! AC(2+3*i) = width
277 ! with i=0 to 6, the index thus runs from 0 to 20!
278 ! Sum of ncssq=7 cos**2 terms put at different radial locations
279 ! and with a windowing around the maxima as current density:
280 ! ac(0) holds ncssq
281 ! ac(0)*(H(x-(ac(1)-ac(2)/2.))-H(x-(ac(1)+ac(2)/2.)))*cos(pi*(x-ac(1))/(ac(2)))**2
282 ! +ac(3)*(H(x-(ac(4)-ac(5)/2.))-H(x-(ac(4)+ac(5)/2.)))*cos(pi*(x-ac(4))/(ac(5)))**2
283 ! + ...
284 ! +ac(18)*(H(x-(ac(19)-ac(20)/2.))-H(x-(ac(19)+ac(20)/2.)))*cos(pi*(x-ac(19))/(ac(20)))**2
285 ! windowing is done with the Heavyside-functions producing the
286 ! Boxcar Function: H(x-a)-H(x-b) is, for a<b, 1 for a<x<b and zero otherwise.
287 !
288  ncssq=7
289  bsta = 2.0_dp
290  bend = 2.0_dp
291  xi = 0.0_dp ! initialize to a finite value
292  wd = 1.0_dp ! initialize to a finite value
293  do i=1,ncssq
294  xi(i)=ac(1+3*(i-1)) ! location of maximum curr. dens.
295  wd(i)=ac(2+3*(i-1)) ! half width of interval
296 ! bsta(i)=max(0.0_dp,xi(i)-wd(i)) ! interval start of curr.dens.part
297 ! bend(i)=min(xi(i)+wd(i),1.0_dp) ! interval end of curr.dens.part
298  bsta(i)=max(0.0_dp,xi(i)-wd(i)) ! interval start of curr.dens.part
299  bend(i)=min(xi(i)+wd(i),1.0_dp) ! interval end of curr.dens.part
300  enddo
301 ! write(6,*) xi(1:ncssq)
302 ! write(6,*) wd(1:ncssq)
303 ! write(6,*) bsta(1:ncssq)
304 ! write(6,*) bend(1:ncssq)
305 
306  pcurr=0.0
307  do i=1,ncssq
308  if(ac(3*(i-1)) .ne. 0.0_dp)then ! do only if the amplitude contributes
309  if(x .gt. bsta(i) .and. x .le. bend(i))then
310  pcurr = pcurr + ac(3*(i-1))*0.5*( (x-bsta(i)) &
311  & + wd(i)/pi*( sin(pi*( x -xi(i))/wd(i) ) - &
312  & sin(pi*(bsta(i)-xi(i))/wd(i) ) ) )
313  elseif(x .gt. bend(i))then
314  pcurr = pcurr + ac(3*(i-1))*0.5*((bend(i)-bsta(i)) &
315  & + wd(i)/pi*( sin(pi*(bend(i)-xi(i))/wd(i) ) - &
316  & sin(pi*(bsta(i)-xi(i))/wd(i) ) ) )
317  endif
318  endif
319  enddo
320 ! write(6,*)x,pcurr
321 ! stop "Test"
322 
323  CASE ('gauss_trunc')
324 ! Truncated Gaussian
325 ! I-prime(s) = ac(0) * (exp(-(s/ac(1))**2) - exp(-(1/ac(1))**2)
326 ! Gauss-Legendre Quadrature to get I(s)
327  DO gli = 1,gln
328  xp = x * glx(gli)
329  pcurr = pcurr + glw(gli) * ac(0) * (exp(-(xp / ac(1)) ** 2) &
330  & - exp(-(1 / ac(1)) ** 2))
331  END DO
332  pcurr = pcurr * x ! correct for x interval
333 
334  CASE ('two_power')
335 ! Two power profile
336 ! I-prime(s) = [1 - s**ac(0)]**ac(1) !! Old as of 2010-05-26
337 ! I-prime(s) = ac(0) * [1 - s**ac(1)]**ac(2) !! Current as of 2010-05-26
338 ! Gauss-Legendre Quadrature to get I(s)
339  DO gli = 1,gln
340  xp = x * glx(gli)
341  pcurr = pcurr + glw(gli) * two_power(xp, ac)
342  END DO
343  pcurr = pcurr * x ! correct for x interval
344 
345  CASE ('two_power_gs')
346 ! Two power with Gaussian peak profile
347 ! I-prime(s) = ac(0) * [1 - s**ac(1)]**ac(2) * (1 + Sum[ac(i)*Exp(-(s - ac(i+1))/ac(i+2)) ** 2])
348 ! Gauss-Legendre Quadrature to get I(s)
349  DO gli = 1,gln
350  xp = x * glx(gli)
351  pcurr = pcurr + glw(gli) * two_power_gs(xp, ac)
352  END DO
353  pcurr = pcurr * x ! correct for x interval
354 
355  CASE('sum_atan')
356 ! I(s) - not I-prime(s)
357 ! pcurr = ac(0) + &
358 ! & ac(1) * (2/pi) * atan(ac(2)*x**ac(3)/(1-x)**ac(4)) + &
359 ! & ac(5) * (2/pi) * atan(ac(6)*x**ac(7)/(1-x)**ac(8)) + &
360 ! & ac(9) * (2/pi) * atan(ac(10)*x**ac(11)/(1-x)**ac(12)) + &
361 ! & ac(13) * (2/pi) * atan(ac(14)*x**ac(15)/(1-x)**ac(16)) + &
362 ! & ac(17) * (2/pi) * atan(ac(18)*x**ac(19)/(1-x)**ac(20))
363  IF (x .ge. one) THEN
364  pcurr = ac(0) + ac(1) + ac(5) + ac(9) + ac(13) + ac(17)
365  ELSE
366  pcurr = ac(0) + &
367  & (2/pi) * (ac(1) * atan(ac(2)*x**ac(3)/(1-x)**ac(4)) + &
368  & ac(5) * atan(ac(6)*x**ac(7)/(1-x)**ac(8)) + &
369  & ac(9) * atan(ac(10)*x**ac(11)/(1-x)**ac(12)) + &
370  & ac(13) * atan(ac(14)*x**ac(15)/(1-x)**ac(16)) + &
371  & ac(17) * atan(ac(18)*x**ac(19)/(1-x)**ac(20)))
372  ENDIF
373 
374  CASE('power_series_I','power_series_i') ! former ipcurr=3
375 ! I(s) - not I-prime(s)
376 ! I(s) = Sum(i,0,-)[ac(i) * s ** (i + 1)]
377 ! Note that coefficient interpretation here is different from other
378 ! polynomial parameterizations, where the terms are a_(i) * s ** i
379 ! Note that I(0) = 0.
380  DO i = ubound(ac,1), ioff, -1
381  pcurr = (pcurr + ac(i))*x
382  END DO
383 
384  CASE('Akima_spline_I','akima_spline_i') ! former ipcurr=11
385 ! I(s) with Akima splines
386 ! Akima-spline for current profile
387 ! i = minloc(scugrid(2:ndatafmax),dim=1) ! this is where zeros are again
388  i = minloc(ac_aux_s(2:),dim=1) ! this is where zeros are again
389  IF (i < 4) stop 'pcurr: check s-grid for curr-grid values!'
390 ! call spline_akima(x,pcurr,scugrid,cugrid,i,iflag)
391  CALL spline_akima(x,pcurr,ac_aux_s,ac_aux_f,i,iflag)
392  IF(iflag < 0) &
393  & stop 'pcurr: outside value from spline_akima requested'
394 
395  CASE('Akima_spline_Ip','akima_spline_ip') ! former ipcurr=12
396 ! I(s) from I-prime(s)-values
397 ! uses an integrated Akima-spline to deliver the values of the toroidal current
398 ! i = minloc(scdgrid(2:ndatafmax),dim=1) ! this is where zeros are again
399  i = minloc(ac_aux_s(2:),dim=1) ! this is where zeros are again
400  IF(i < 4)stop 'pcurr: check s-grid for curr.dens.-grid values!'
401 ! call spline_akima_int(x,pcurr,scdgrid,cdgrid,i,iflag)
402  CALL spline_akima_int(x,pcurr,ac_aux_s,ac_aux_f,i,iflag)
403  IF (iflag < 0) THEN
404  WRITE(6,*)"x = ", x
405  stop 'pcurr: outside value from spline_akima_int requested'
406  ENDIF
407 ! Note that spline_akima_int integrates with ac_aux_s(1) as the lower limit.
408 ! Assumption here is that lower limit is s = 0. Check this assumption.
409  IF (abs(ac_aux_s(1)) .gt. 0.0001) THEN
410  WRITE(*,*) ' Error in profile_functions, Akima_spline_Ip'
411  WRITE(*,*) .gt.' ABS(ac_aux_s(1)) 0.0001'
412  WRITE(*,*) ' Fix This. Stopping'
413  stop
414  ENDIF
415 
416  CASE('cubic_spline_I','cubic_spline_i') ! former ipcurr=13
417 ! I(s) with Cubic splines
418 ! Cubic-spline for current profile
419 ! i = minloc(scugrid(2:ndatafmax),dim=1) ! this is where zeros are again
420  i = minloc(ac_aux_s(2:),dim=1) ! this is where zeros are again
421  IF(i < 4) stop 'pcurr: check s-grid for curr-grid values!'
422 ! call spline_cubic(x,pcurr,scugrid,cugrid,i,iflag)
423  CALL spline_cubic(x,pcurr,ac_aux_s,ac_aux_f,i,iflag)
424  IF (iflag < 0) THEN
425  IF(iflag == -1) THEN
426  stop'pcurr: outside value from spline_cubic requested'
427  ELSEIF(iflag == -2) THEN
428  stop 'pcurr: identical scugrid-values in spline_cubic'
429  ELSE
430  stop 'pcurr: unknown error from spline_cubic'
431  ENDIF
432  ENDIF
433 
434  CASE('cubic_spline_Ip','cubic_spline_ip') ! former ipcurr=14
435 ! I(s) from I-prime(s)-values
436 ! uses an integrated Cubic-spline to deliver the values of the toroidal current
437 ! i = minloc(scdgrid(2:ndatafmax),dim=1) ! this is where zeros are again
438  i = minloc(ac_aux_s(2:),dim=1) ! this is where zeros are again
439  IF(i < 4) stop 'pcurr: check s-grid for curr.dens.-grid values!'
440 ! call spline_cubic_int(x,pcurr,scdgrid,cdgrid,i,iflag)
441  CALL spline_cubic_int(x,pcurr,ac_aux_s,ac_aux_f,i,iflag)
442  IF (iflag < 0)THEN
443  IF (iflag == -1) THEN
444  stop'pcurr: outside value from spline_cubic_int requested'
445  ELSEIF (iflag == -2) THEN
446  stop 'pcurr: identical scdgrid-values in spline_cubic_int'
447  ELSE
448  stop 'pcurr: unknown error from spline_cubic_int'
449  ENDIF
450  ENDIF
451 ! Note that spline_cubic_int integrates with ac_aux_s(1) as the lower limit.
452 ! Assumption here is that lower limit is s = 0. Check this assumption.
453  IF (abs(ac_aux_s(1)) .gt. 0.0001) THEN
454  WRITE(*,*) ' Error in profile_functions, cubic_spline_Ip'
455  WRITE(*,*) .gt.' ABS(ac_aux_s(1)) 0.0001'
456  WRITE(*,*) ' Fix This. Stopping'
457  stop
458  ENDIF
459 
460  CASE ('pedestal')
461 ! Parameterization for I(s)
462 ! Added by SPH 2010-05-26
463  DO i = 7, ioff, -1
464  pcurr = x*pcurr + ac(i)/(i-ioff+1)
465  END DO
466  pcurr = x*pcurr
467 
468  i=8
469  IF (ac(i+3) .le. 0._dp ) THEN
470  ac(i:i+4) = 0
471  ac(i+3) = 1.0e30_dp
472  ELSE
473  ac(i+4) = 1.0_dp/
474  1 (tanh(2*ac(i+2)/ac(i+3))-tanh(2*(ac(i+2)-1)/ac(i+3)))
475  END IF
476 
477  a8 =max(ac(i+8),0.01_dp)
478  a12=max(ac(i+12),0.01_dp)
479 
480  g1= (x-ac(i+7))/a8
481  g3= (-ac(i+7))/a8
482  g2= (x-ac(i+11))/a12
483  g4= (-ac(i+11))/a12
484  pcurr = pcurr + ac(i+4) * ac(i+0) *
485  1 ( tanh( 2*ac(i+2)/ac(i+3) )
486  2 -tanh( 2*(ac(i+2)-sqrt(x))/ac(i+3) ) )
487  3 +ac(i+5)*( tanh(g1) - tanh(g3) )
488  4 +ac(i+9)*( tanh(g2) - tanh(g4) )
489 
490 
491  CASE('rational') !
492 ! I(s) - not I-prime(s). No need for Gaussian quadrature.
493 ! Ratio of polynomials.
494 ! Numerator coefficients: ac(LBOUND) - ac(9)
495 ! Denominator coefficients: ac(10) - ac(UBOUND(ac,1))
496  temp_num = zero
497  temp_denom = zero
498  DO i = 9, ioff, -1
499  temp_num = x * temp_num + ac(i)
500  END DO
501  DO i = ubound(ac,1), 10, -1
502  temp_denom = x * temp_denom + ac(i)
503  END DO
504  IF (temp_denom .ne. zero) THEN
505  pcurr = temp_num / temp_denom
506  ELSE
507  pcurr = huge(pcurr)
508  ENDIF
509 
510  CASE('line_segment_Ip', 'line_segment_ip')
511 ! I(s) from I-prime(s)-values. Integrated line segments to determine the plasma
512 ! current.
513  i = minloc(ac_aux_s(2:),dim=1)
514  CALL line_seg_int(x,pcurr,ac_aux_s,ac_aux_f,i)
515 
516  CASE('line_segment_I', 'line_segment_i')
517 ! I(s) values. Linearly interpolated line segments to determine the plasma
518 ! current.
519  i = minloc(ac_aux_s(2:),dim=1)
520  CALL line_seg(x,pcurr,ac_aux_s,ac_aux_f,i)
521 
522  CASE DEFAULT
523 ! Power series
524 ! I-prime(s) = Sum(i,0,-)[ac(i) * s ** i]
525 ! Analytic integration to get I(s)
526  DO i = ubound(ac,1), ioff, -1
527  pcurr = x*pcurr + ac(i)/(i-ioff+1)
528  END DO
529  IF (trim(pcurr_type_lc) .ne. 'power_series') THEN
530  WRITE(*,*) 'Unrecognized pcurr_type:', pcurr_type
531  WRITE(*,*) ' *** CHECK YOUR INPUT ***'
532  WRITE(*,*) 'Changing pcurr_type from ''', &
533  & trim(pcurr_type), '''to ''power_series''.'
534  pcurr_type = 'power_series'
535  END IF
536  pcurr = x*pcurr
537  END SELECT
538 
539  END FUNCTION pcurr
540 
541 !******************************************************************************
542 
543  FUNCTION piota (x)
544 ! Function to compute the iota/q profile.
545 
546 ! Variables declared in vmec_input:
547 ! ai array (0:20) of coefficients
548 ! ai_aux_s Auxiliary array s, s-values used for splines
549 ! ai_aux_f Auxiliary array f, function-values used for splines
550 ! piota_type character, specifies the parameterization of the profile
551 ! sum_atan Sum of atan functions plus offset.
552 ! Akima_spline Akima spline
553 ! cubic_spline cubic spline
554 ! rational rational (ratio of polynomials)
555 ! nice_quadratic quadratic with rerranged coefficients.
556 ! line_segment Linearly interpolated line segments
557 ! power_series Power Series (Default)
558 ! lRFP logical - Reversed Field Pinch
559 ! True - ai parameterization specifies q-profile (= 1 / iota)
560 ! False - ai parameterization is for iota-profile
561 
562 ! Local Variables
563 ! i integer counter
564 ! iflag error flag for spline call
565 ! x real argument
566 ! x constrained to be between 0 and 1
567 ! piota_type_lc character, piota_type -> lower case
568 
569  USE stel_kinds
570  USE stel_constants, ONLY: zero, one, pi
571  USE vmec_input, ONLY: ai, piota_type, ai_aux_s, ai_aux_f, lrfp
572  USE line_segment
573  IMPLICIT NONE
574 ! ai assumed to be dimensioned (0:n), with n >= 20
575 !-----------------------------------------------
576  INTEGER :: i, iflag, ioff
577  REAL(rprec), INTENT(IN) :: x
578  REAL(rprec) :: piota, temp_num, temp_denom
579  CHARACTER(len=20) :: piota_type_lc
580 !-----------------------------------------------
581  piota = 0
582  ioff = lbound(ai,1) ! Expected to be zero.
583 
584 ! Convert to Lower Case, to avoid typo issues
585  piota_type_lc = piota_type
586  CALL tolower(piota_type_lc)
587  SELECT CASE(trim(piota_type_lc))
588 
589  CASE('sum_atan')
590 ! Sum atan functions mapped to [0:1], with an ai(0) offset
591 ! piota = ai(0) + &
592 ! & ai(1) * (2/pi) * atan(ai(2)*x**ai(3)/(1-x)**ai(4)) + &
593 ! & ai(5) * (2/pi) * atan(ai(6)*x**ai(7)/(1-x)**ai(8)) + &
594 ! & ai(9) * (2/pi) * atan(ai(10)*x**ai(11)/(1-x)**ai(12)) + &
595 ! & ai(13) * (2/pi) * atan(ai(14)*x**ai(15)/(1-x)**ai(16)) + &
596 ! & ai(17) * (2/pi) * atan(ai(18)*x**ai(19)/(1-x)**ai(20))
597  IF (x .ge. one) THEN
598  piota = ai(0) + ai(1) + ai(5) + ai(9) + ai(13) + ai(17)
599  ELSE
600  piota = ai(0) + &
601  & (2/pi) * (ai(1) * atan(ai(2)*x**ai(3)/(1-x)**ai(4)) + &
602  & ai(5) * atan(ai(6)*x**ai(7)/(1-x)**ai(8)) + &
603  & ai(9) * atan(ai(10)*x**ai(11)/(1-x)**ai(12)) + &
604  & ai(13) * atan(ai(14)*x**ai(15)/(1-x)**ai(16)) + &
605  & ai(17) * atan(ai(18)*x**ai(19)/(1-x)**ai(20)))
606  ENDIF
607 
608  CASE('Akima_spline','akima_spline') ! former ipiota=11
609 ! Akima-spline (iogrid + siogrid)
610 ! i = minloc(siogrid(2:ndatafmax),dim=1) ! this is where zeros are again
611  i = minloc(ai_aux_s(2:),dim=1) ! this is where zeros are again
612  if(i < 4) stop 'piota: check s-grid for iota-grid values!'
613 ! call spline_akima(x,piota,siogrid,iogrid,i,iflag)
614  call spline_akima(x,piota,ai_aux_s,ai_aux_f,i,iflag)
615  if(iflag < 0) &
616  & stop'piota: outside value from spline_akima requested'
617 
618  CASE('cubic_spline') ! former ipiota=13
619 ! cubic-spline (iogrid + siogrid)
620 ! i = minloc(siogrid(2:ndatafmax),dim=1) ! this is where zeros are again
621  i = minloc(ai_aux_s(2:),dim=1) ! this is where zeros are again
622  if(i < 4) stop 'piota: check s-grid for iota-grid values!'
623 ! call spline_cubic(x,piota,siogrid,iogrid,i,iflag)
624  call spline_cubic(x,piota,ai_aux_s,ai_aux_f,i,iflag)
625  if(iflag < 0)then
626  if(iflag == -1) then
627  stop'piota: outside value from spline_cubic requested'
628  elseif(iflag == -2) then
629  stop'piota: identical siogrid-values in spline_cubic'
630  else
631  stop 'piota: unknown error from spline_cubic'
632  endif
633  endif
634 
635  CASE('rational') !
636 ! Ratio of polynomials.
637 ! Numerator coefficients: ai(LBOUND) - ai(9)
638 ! Denominator coefficients: ai(10) - ai(UBOUND)
639  temp_num = zero
640  temp_denom = zero
641  DO i = 9, ioff, -1
642  temp_num = x * temp_num + ai(i)
643  END DO
644  DO i = ubound(ai,1), 10, -1
645  temp_denom = x * temp_denom + ai(i)
646  END DO
647  IF (temp_denom .ne. zero) THEN
648  piota = temp_num / temp_denom
649  ELSE
650  piota = huge(piota)
651  ENDIF
652 
653  CASE('nice_quadratic') !
654 ! Quadratic with slightly rearranged coefficients.
655 ! iota(s) = a0(1-s) + a1 s + 4 a2 s (1 - s)
656 ! a0 is the iota value at s = 0
657 ! a1 is the iota value at s = 1
658 ! a2 is the shift of iota(1/2) from the straight line value
659 ! (Thus, a2 = 0 gives a linear iota profile.
660  piota = ai(0) * (one - x) + ai(1) * x + 4 * ai(2) * &
661  & x * (one - x)
662 
663  CASE('line_segment')
664 ! Linearly interpolated line segments to determine the iotabar profile.
665  i = minloc(ai_aux_s(2:),dim=1)
666  CALL line_seg(x,piota,ai_aux_s,ai_aux_f,i)
667 
668  CASE default
669 ! Power Series is the default
670  DO i = ubound(ai,1), lbound(ai,1), -1
671  piota = x*piota + ai(i)
672  END DO
673  IF (trim(piota_type_lc) .ne. 'power_series') THEN
674  WRITE(*,*) 'Unrecognized piota_type:', piota_type
675  WRITE(*,*) ' *** CHECK YOUR INPUT ***'
676  WRITE(*,*) 'Changing piota_type from ''', &
677  & trim(piota_type), '''to ''power_series''.'
678  piota_type = 'power_series'
679  END IF
680  END SELECT
681 
682  IF (lrfp) THEN
683 ! RFP/TOKAMAK: ai are coefficients of q_safety, NOT iota
684  IF (piota .ne. zero) THEN
685  piota = one / piota
686  ELSE
687  piota = huge(piota)
688  END IF
689  END IF
690 
691  END FUNCTION piota
692 
693 !******************************************************************************
694 
695  FUNCTION pmass (xx)
696 ! Function to compute the mass (pressure) profile
697 
698 ! Variables declared in vmec_input:
699 ! am array (0:20) of coefficients
700 ! am_aux_s Auxiliary array s, s-values used for splines
701 ! am_aux_f Auxiliary array f, function-values used for splines
702 ! bloat used to expand the profile
703 ! pmass_type character, specifies the parameterization of the profile
704 ! gauss_trunc Truncated Gaussian
705 ! two_power am(0) * [1 - s**am(1)]**am(2)
706 ! two_power_gs am(0) * [1 - s**am(1)]**am(2)*
707 ! (1 + Sum[am(i)*Exp(-((x - am(i+1))/am(i+2))**2)])
708 ! two_Lorentz two Lorentz-type functions, mapped to [0:1]
709 ! Akima_spline Akima-spline (magrid[-> am_aux_f] + smagrid[-> am_aux_s])
710 ! cubic_spline Cubic-spline (magrid[-> am_aux_f] + smagrid[-> am_aux_s])
711 ! pedestal Pedestal profile
712 ! rational rational (ratio of polynomials)
713 ! line_segment Linearly interpolated line segments
714 ! power_series Power Series (Default)
715 
716 ! Local Variables
717 ! i integer counter
718 ! iflag error flag for spline call
719 ! xx real argument
720 ! x constrained to be between 0 and 1
721 ! pmass_type_lc character, pmass_type -> lower case
722 
723  USE stel_kinds
724  USE stel_constants, ONLY: zero, one
725  USE vmec_input, ONLY: am, bloat, pres_scale, pmass_type, &
726  & am_aux_s, am_aux_f
727 ! am is assumed to be dimensioned starting at zero.
728  USE vparams, ONLY: mu0
729  USE line_segment
730  USE functions
731  IMPLICIT NONE
732 !-----------------------------------------------
733  INTEGER :: i, iflag, ioff
734  REAL(rprec) :: xx, pmass, x, temp_num, temp_denom
735  CHARACTER(len=20) :: pmass_type_lc
736 !-----------------------------------------------
737 ! NOTE: On entry, am is in pascals. pmass internal units are mu0*pascals (B**2 units)
738  x = min(abs(xx * bloat), 1._dp)
739 
740  pmass = 0
741  ioff = lbound(am,1) ! Expected to be zero.
742 
743 ! Convert to Lower Case, to avoid typo issues
744  pmass_type_lc = pmass_type
745  CALL tolower(pmass_type_lc)
746  SELECT CASE(trim(pmass_type_lc))
747 
748  CASE ('gauss_trunc')
749 ! Truncated Gaussian
750 ! p(s) = am(0) * (exp(-(s/am(1))**2) - exp(-(1/am(1))**2)
751 ! Above Old as of 2010-05-26
752 ! Below current as of 2010-05-26
753 ! p(s) = (am(0) / c) * (exp(-(s/am(1))**2) - exp(-(1/am(1))**2)
754 ! c = (1 - exp(-(1/am(1))**2) ! so that p(0) = am(0).
755  pmass = (am(0)/(one - exp(-(one / am(1)) ** 2))) * &
756  & (exp(-(x / am(1)) ** 2) - exp(-(one / am(1)) ** 2))
757 
758  CASE ('two_power')
759 ! Two power profile
760 ! p(s) = [1 - s**am(0)]**am(1) !! Old as of 2010-05-26
761 ! p(s) = am(0) * [1 - s**am(1)]**am(2) !! Current as of 2010-05-26
762  pmass = two_power(x, am)
763 
764  CASE ('two_power_gs')
765 ! Two power profile with guassian peaks.
766  pmass = two_power_gs(x, am)
767 
768  CASE ('two_Lorentz','two_lorentz') ! former ipmass=1
769  pmass = am(0)*(am(1)*(one/(one+( x/am(2)**2)**am(3))**am(4) &
770  & -one/(one+(one/am(2)**2)**am(3))**am(4))/ &
771  & (one-one/(one+(one/am(2)**2)**am(3))**am(4))+ &
772  & (one-am(1))*(one/(one+( x/am(5)**2)**am(6))**am(7) &
773  & -one/(one+(one/am(5)**2)**am(6))**am(7))/ &
774  & (one-one/(one+(one/am(5)**2)**am(6))**am(7))) &
775 
776  CASE ('Akima_spline','akima_spline') ! former ipmass=11
777 ! i = minloc(smagrid(2:ndatafmax),dim=1) ! this is where zeros are again
778  i = minloc(am_aux_s(2:),dim=1) ! this is where zeros are again
779  IF(i < 4) stop 'pmass: check s-grid for mass-grid values!'
780 ! call spline_akima(x,pmass,smagrid,magrid,i,iflag)
781  CALL spline_akima(x,pmass,am_aux_s,am_aux_f,i,iflag)
782  IF(iflag < 0) &
783  & stop 'pmass: outside value from spline_akima requested'
784 
785  CASE ('cubic_spline') ! former ipmass=13
786 ! i = minloc(smagrid(2:ndatafmax),dim=1) ! this is where zeros are again
787  i = minloc(am_aux_s(2:),dim=1) ! this is where zeros are again
788  IF(i < 4) stop 'pmass: check s-grid for mass-grid values!'
789 ! call spline_cubic(x,pmass,smagrid,magrid,i,iflag)
790  CALL spline_cubic(x,pmass,am_aux_s,am_aux_f,i,iflag)
791  IF (iflag < 0) THEN
792  IF (iflag == -1) THEN
793  stop'pmass: outside value from spline_cubic requested'
794  ELSEIF (iflag == -2) THEN
795 ! STOP'pmass: identical smagrid-values in spline_cubic'
796  stop'pmass: identical am_aux_s-values in spline_cubic'
797  ELSE
798  stop 'pmass: unknown error from spline_cubic'
799  ENDIF
800  ENDIF
801 
802  CASE ('pedestal') !PMV Texas group
803 ! Added by SPH 2010-05-26
804 
805  DO i = 15, lbound(am,1), -1
806  pmass = x*pmass + am(i)
807  END DO
808 
809  i = 16
810  IF (am(i+3) .le. 0._dp ) THEN
811  am(i:i+4) = 0
812  am(i+3) = 1.0e30_dp
813  ELSE
814  am(i+4) = 1.0_dp/
815  1 (tanh(2*am(i+2)/am(i+3))-tanh(2*(am(i+2)-1)/am(i+3)))
816  END IF
817 
818  pmass = pmass + am(i+4) * am(i+1) *
819  1 ( tanh( 2*(am(i+2)-sqrt(x))/am(i+3) )
820  2 -tanh( 2*(am(i+2)-1._dp) /am(i+3) ) )
821 
822  CASE('rational') !
823 ! Ratio of polynomials.
824 ! Numerator coefficients: am(LBOUND) - am(9)
825 ! Denominator coefficients: am(10) - am(UBOUND)
826  temp_num = zero
827  temp_denom = zero
828  DO i = 9, ioff, -1
829  temp_num = x * temp_num + am(i)
830  END DO
831  DO i = ubound(am,1), 10, -1
832  temp_denom = x * temp_denom + am(i)
833  END DO
834  IF (temp_denom .ne. zero) THEN
835  pmass = temp_num / temp_denom
836  ELSE
837  pmass = huge(pmass)
838  ENDIF
839 
840  CASE('line_segment')
841 ! Linearly interpolated line segments to determine the pressure profile.
842  i = minloc(am_aux_s(2:),dim=1)
843  CALL line_seg(x,pmass,am_aux_s,am_aux_f,i)
844 
845  CASE DEFAULT
846  DO i = ubound(am,1), lbound(am,1), -1
847  pmass = x*pmass + am(i)
848  END DO
849  IF (trim(pmass_type_lc) .ne. 'power_series') THEN
850  WRITE(*,*) 'Unrecognized pmass_type:', pmass_type
851  WRITE(*,*) ' *** CHECK YOUR INPUT ***'
852  WRITE(*,*) 'Changing pmass_type from ''', &
853  & trim(pmass_type), '''to ''power_series''.'
854  pmass_type = 'power_series'
855  END IF
856  END SELECT
857 
858  pmass = mu0*pres_scale*pmass
859 
860  END FUNCTION pmass
861 
862 #ifdef _ANIMEC
863  FUNCTION photp (xx)
864  USE stel_kinds
865  USE vmec_input, ONLY: ah, bloat
866 C-----------------------------------------------
867  INTEGER :: i
868  REAL(rprec) :: xx, photp, x
869 C-----------------------------------------------
870 ! NOTE: On entry, ah is dimensionless
871 ! HOT PARTICLE PRESSURE (RATIO TO ISOTROPIC PRESSURE)
872 
873  x = min(abs(xx * bloat), 1._dp)
874 
875  photp = 0
876 
877  DO i = ubound(ah,1), lbound(ah,1), -1
878  photp = x*photp + ah(i)
879  END DO
880 
881  photp = photp
882 
883  END FUNCTION photp
884 
885  FUNCTION ptrat (xx)
886  USE stel_kinds
887  USE vmec_input, ONLY: at, bloat
888 C-----------------------------------------------
889  INTEGER :: i
890  REAL(rprec) :: xx, ptrat, x
891 C-----------------------------------------------
892 ! NOTE: On entry, at is dimensionless
893 ! HOT PARTICLE T-perp/T-par
894 
895  x = min(abs(xx * bloat), 1._dp)
896 
897  ptrat = 0
898 
899  DO i = ubound(at,1), lbound(at,1), -1
900  ptrat = x*ptrat + at(i)
901  END DO
902 
903  ptrat = ptrat
904 
905  END FUNCTION ptrat
906 #endif
line_segment
Module is part of the LIBSTELL. This module contains code to create a profile constructed of line sig...
Definition: line_segment.f:13
functions
This module containes functions used by the profiles.
Definition: functions.f:10
line_segment::line_seg
subroutine, public line_seg(x, y, xx, yy, n)
Interpolate a point on a line.
Definition: line_segment.f:41
line_segment::line_seg_int
subroutine, public line_seg_int(x, y, xx, yy, n)
Integrate to a point on a line.
Definition: line_segment.f:88