56 USE stel_constants,
ONLY: zero, one, pi
57 USE vmec_input,
ONLY: ac, bloat, pcurr_type, ac_aux_s, ac_aux_f
63 INTEGER :: i, ioff, iflag
64 REAL(rprec) :: xx, pcurr, x, xp, temp_num, temp_denom
65 CHARACTER(len=20) :: pcurr_type_lc
67 INTEGER,
PARAMETER :: gln = 10
69 REAL(rprec),
DIMENSION(gln),
PARAMETER :: glx = (/
74 REAL(rprec),
DIMENSION(gln),
PARAMETER :: glw = (/
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
92 x = min(abs(xx * bloat), one)
98 pcurr_type_lc = pcurr_type
99 CALL tolower(pcurr_type_lc)
100 SELECT CASE(trim(pcurr_type_lc))
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!"
143 delx=1.0_dp/(ncssq-1.0_dp)
148 bsta(i)=max(0.0_dp,xi(i)-delx)
149 bend(i)=min(xi(i)+delx,1.0_dp)
153 if( x .le. bend(i) .and. x .gt. bsta(i))
then
159 if( x .gt. bsta(i) .and. x .le. bend(i) )
then
161 pcurr=pcurr + 0.5*ac(i)*(x+
162 & delx*sin(pi*(x-xi(i))/delx)/pi)
164 pcurr=pcurr + 0.5*ac(i)*(x-xi(i)+delx+
165 & delx*sin(pi*(x-xi(i))/delx)/pi)
169 pcurr=pcurr+0.5*ac(i)*delx
171 pcurr=pcurr+ac(i)*delx
176 CASE (
'sum_cossq_sqrts')
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!"
221 delx=1.0_dp/(ncssq-1.0_dp)
228 bsta(i)=max(0.0_dp,xi(i)-delx)
229 bend(i)=min(xi(i)+delx,1.0_dp)
233 if( sqx .le. bend(i) .and. sqx .gt. bsta(i))
then
239 if( sqx .gt. bsta(i) .and. sqx .le. bend(i) )
then
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))
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))
254 pcurr=pcurr+ac(i)*(0.25_dp-1/pisq)*delxsq
256 pcurr=pcurr+ac(i)*delx*xi(i)
261 CASE (
'sum_cossq_s_free')
298 bsta(i)=max(0.0_dp,xi(i)-wd(i))
299 bend(i)=min(xi(i)+wd(i),1.0_dp)
308 if(ac(3*(i-1)) .ne. 0.0_dp)
then
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) ) ) )
329 pcurr = pcurr + glw(gli) * ac(0) * (exp(-(xp / ac(1)) ** 2)
330 & - exp(-(1 / ac(1)) ** 2))
341 pcurr = pcurr + glw(gli) * two_power(xp, ac)
345 CASE (
'two_power_gs')
351 pcurr = pcurr + glw(gli) * two_power_gs(xp, ac)
364 pcurr = ac(0) + ac(1) + ac(5) + ac(9) + ac(13) + ac(17)
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)))
374 CASE(
'power_series_I',
'power_series_i')
380 DO i = ubound(ac,1), ioff, -1
381 pcurr = (pcurr + ac(i))*x
384 CASE(
'Akima_spline_I',
'akima_spline_i')
388 i = minloc(ac_aux_s(2:),dim=1)
389 IF (i < 4) stop
'pcurr: check s-grid for curr-grid values!'
391 CALL spline_akima(x,pcurr,ac_aux_s,ac_aux_f,i,iflag)
393 & stop
'pcurr: outside value from spline_akima requested'
395 CASE(
'Akima_spline_Ip',
'akima_spline_ip')
399 i = minloc(ac_aux_s(2:),dim=1)
400 IF(i < 4)stop
'pcurr: check s-grid for curr.dens.-grid values!'
402 CALL spline_akima_int(x,pcurr,ac_aux_s,ac_aux_f,i,iflag)
405 stop
'pcurr: outside value from spline_akima_int requested'
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'
416 CASE(
'cubic_spline_I',
'cubic_spline_i')
420 i = minloc(ac_aux_s(2:),dim=1)
421 IF(i < 4) stop
'pcurr: check s-grid for curr-grid values!'
423 CALL spline_cubic(x,pcurr,ac_aux_s,ac_aux_f,i,iflag)
426 stop
'pcurr: outside value from spline_cubic requested'
427 ELSEIF(iflag == -2)
THEN
428 stop
'pcurr: identical scugrid-values in spline_cubic'
430 stop
'pcurr: unknown error from spline_cubic'
434 CASE(
'cubic_spline_Ip',
'cubic_spline_ip')
438 i = minloc(ac_aux_s(2:),dim=1)
439 IF(i < 4) stop
'pcurr: check s-grid for curr.dens.-grid values!'
441 CALL spline_cubic_int(x,pcurr,ac_aux_s,ac_aux_f,i,iflag)
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'
448 stop
'pcurr: unknown error from spline_cubic_int'
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'
464 pcurr = x*pcurr + ac(i)/(i-ioff+1)
469 IF (ac(i+3) .le. 0._dp )
THEN
474 1 (tanh(2*ac(i+2)/ac(i+3))-tanh(2*(ac(i+2)-1)/ac(i+3)))
477 a8 =max(ac(i+8),0.01_dp)
478 a12=max(ac(i+12),0.01_dp)
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) )
499 temp_num = x * temp_num + ac(i)
501 DO i = ubound(ac,1), 10, -1
502 temp_denom = x * temp_denom + ac(i)
504 IF (temp_denom .ne. zero)
THEN
505 pcurr = temp_num / temp_denom
510 CASE(
'line_segment_Ip',
'line_segment_ip')
513 i = minloc(ac_aux_s(2:),dim=1)
516 CASE(
'line_segment_I',
'line_segment_i')
519 i = minloc(ac_aux_s(2:),dim=1)
520 CALL line_seg(x,pcurr,ac_aux_s,ac_aux_f,i)
526 DO i = ubound(ac,1), ioff, -1
527 pcurr = x*pcurr + ac(i)/(i-ioff+1)
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'
570 USE stel_constants,
ONLY: zero, one, pi
571 USE vmec_input,
ONLY: ai, piota_type, ai_aux_s, ai_aux_f, lrfp
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
585 piota_type_lc = piota_type
586 CALL tolower(piota_type_lc)
587 SELECT CASE(trim(piota_type_lc))
598 piota = ai(0) + ai(1) + ai(5) + ai(9) + ai(13) + ai(17)
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)))
608 CASE(
'Akima_spline',
'akima_spline')
611 i = minloc(ai_aux_s(2:),dim=1)
612 if(i < 4) stop
'piota: check s-grid for iota-grid values!'
614 call spline_akima(x,piota,ai_aux_s,ai_aux_f,i,iflag)
616 & stop
'piota: outside value from spline_akima requested'
621 i = minloc(ai_aux_s(2:),dim=1)
622 if(i < 4) stop
'piota: check s-grid for iota-grid values!'
624 call spline_cubic(x,piota,ai_aux_s,ai_aux_f,i,iflag)
627 stop
'piota: outside value from spline_cubic requested'
628 elseif(iflag == -2)
then
629 stop
'piota: identical siogrid-values in spline_cubic'
631 stop
'piota: unknown error from spline_cubic'
642 temp_num = x * temp_num + ai(i)
644 DO i = ubound(ai,1), 10, -1
645 temp_denom = x * temp_denom + ai(i)
647 IF (temp_denom .ne. zero)
THEN
648 piota = temp_num / temp_denom
653 CASE(
'nice_quadratic')
660 piota = ai(0) * (one - x) + ai(1) * x + 4 * ai(2) *
665 i = minloc(ai_aux_s(2:),dim=1)
666 CALL line_seg(x,piota,ai_aux_s,ai_aux_f,i)
670 DO i = ubound(ai,1), lbound(ai,1), -1
671 piota = x*piota + ai(i)
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'
684 IF (piota .ne. zero)
THEN
724 USE stel_constants,
ONLY: zero, one
725 USE vmec_input,
ONLY: am, bloat, pres_scale, pmass_type, &
728 USE vparams,
ONLY: mu0
733 INTEGER :: i, iflag, ioff
734 REAL(rprec) :: xx, pmass, x, temp_num, temp_denom
735 CHARACTER(len=20) :: pmass_type_lc
738 x = min(abs(xx * bloat), 1._dp)
744 pmass_type_lc = pmass_type
745 CALL tolower(pmass_type_lc)
746 SELECT CASE(trim(pmass_type_lc))
755 pmass = (am(0)/(one - exp(-(one / am(1)) ** 2))) *
756 & (exp(-(x / am(1)) ** 2) - exp(-(one / am(1)) ** 2))
762 pmass = two_power(x, am)
764 CASE (
'two_power_gs')
766 pmass = two_power_gs(x, am)
768 CASE (
'two_Lorentz',
'two_lorentz')
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)))
776 CASE (
'Akima_spline',
'akima_spline')
778 i = minloc(am_aux_s(2:),dim=1)
779 IF(i < 4) stop
'pmass: check s-grid for mass-grid values!'
781 CALL spline_akima(x,pmass,am_aux_s,am_aux_f,i,iflag)
783 & stop
'pmass: outside value from spline_akima requested'
785 CASE (
'cubic_spline')
787 i = minloc(am_aux_s(2:),dim=1)
788 IF(i < 4) stop
'pmass: check s-grid for mass-grid values!'
790 CALL spline_cubic(x,pmass,am_aux_s,am_aux_f,i,iflag)
792 IF (iflag == -1)
THEN
793 stop
'pmass: outside value from spline_cubic requested'
794 ELSEIF (iflag == -2)
THEN
796 stop
'pmass: identical am_aux_s-values in spline_cubic'
798 stop
'pmass: unknown error from spline_cubic'
805 DO i = 15, lbound(am,1), -1
806 pmass = x*pmass + am(i)
810 IF (am(i+3) .le. 0._dp )
THEN
815 1 (tanh(2*am(i+2)/am(i+3))-tanh(2*(am(i+2)-1)/am(i+3)))
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) ) )
829 temp_num = x * temp_num + am(i)
831 DO i = ubound(am,1), 10, -1
832 temp_denom = x * temp_denom + am(i)
834 IF (temp_denom .ne. zero)
THEN
835 pmass = temp_num / temp_denom
842 i = minloc(am_aux_s(2:),dim=1)
843 CALL line_seg(x,pmass,am_aux_s,am_aux_f,i)
846 DO i = ubound(am,1), lbound(am,1), -1
847 pmass = x*pmass + am(i)
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'
858 pmass = mu0*pres_scale*pmass
865 USE vmec_input,
ONLY: ah, bloat
868 REAL(rprec) :: xx, photp, x
873 x = min(abs(xx * bloat), 1._dp)
877 DO i = ubound(ah,1), lbound(ah,1), -1
878 photp = x*photp + ah(i)
887 USE vmec_input,
ONLY: at, bloat
890 REAL(rprec) :: xx, ptrat, x
895 x = min(abs(xx * bloat), 1._dp)
899 DO i = ubound(at,1), lbound(at,1), -1
900 ptrat = x*ptrat + at(i)