V3FIT
vmec_info.f90
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
9 !*******************************************************************************
10  MODULE vmec_info
11  USE stel_kinds
12  USE read_wout_mod, ns_vmec=>ns, mpol_vmec=>mpol, ntor_vmec=>ntor, &
13  rmnc_vmec=>rmnc, zmns_vmec=>zmns, lmns_vmec=>lmns, &
14  xm_vmec=>xm, xn_vmec=>xn, chipf_vmec=>chipf, & ! MRC 4/1/2016
15  rmns_vmec=>rmns, zmnc_vmec=>zmnc, lmnc_vmec=>lmnc, & ! MRC 12/1/2016
16  phipf_vmec=>phipf, presf_vmec=>presf, nfp_vmec=>nfp, &
17  wb_vmec=>wb, wp_vmec=>wp, gamma_vmec=>gamma, &
18  volume_vmec=>volume, raxis_vmec=>raxis, &
19  lasym_vmec=>lasym, iasym_vmec=>iasym, &
20  vmec_curtor=>itor
21  USE fourier, ONLY: f_cos, f_sin
22 
23  IMPLICIT NONE
24 
25 !*******************************************************************************
26 ! vmec_info module variables
27 !*******************************************************************************
29  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: rmnc_spline
31  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: zmns_spline
33  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: lmns_spline
35  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: rmns_spline
37  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: zmnc_spline
39  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: lmnc_spline
41  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: currumnc_spline
43  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: currvmnc_spline
45  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: currumns_spline
47  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: currvmns_spline
48 
50  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: rmnc_i
52  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: zmns_i
54  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: lmns_i
56  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: rmns_i
58  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: zmnc_i
60  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: lmnc_i
62  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: jcurrumnc
64  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: jcurrvmnc
66  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: jcurrumns
68  REAL (dp), DIMENSION(:,:,:), ALLOCATABLE :: jcurrvmns
69 
70 !*******************************************************************************
71 ! vmec_info module parameters
72 !*******************************************************************************
74  INTEGER, PARAMETER :: iflipj = 1
76  REAL (dp), PARAMETER :: lower_b = -1.e30_dp
78  REAL (dp), PARAMETER :: upper_b = -1.e30_dp
79 
80  CONTAINS
81 !*******************************************************************************
82 ! CONSTRUCTION SUBROUTINES
83 !*******************************************************************************
84 !-------------------------------------------------------------------------------
94 !-------------------------------------------------------------------------------
95  SUBROUTINE vmec_info_construct_spline(mnmax, ns, is_asym)
96 
97  IMPLICIT NONE
98 
99 ! Declare Arguments
100  INTEGER, INTENT(in) :: mnmax
101  INTEGER, INTENT(in) :: ns
102  LOGICAL, INTENT(in) :: is_asym
103 
104 ! Start of executable code
105 
106 ! Stellarator symmetric quantities.
107  ALLOCATE(rmnc_spline(mnmax,ns))
108  ALLOCATE(zmns_spline(mnmax,ns))
109  ALLOCATE(lmns_spline(mnmax,ns + 1))
110  ALLOCATE(currumnc_spline(mnmax,ns))
111  ALLOCATE(currvmnc_spline(mnmax,ns))
112 
113 ! Asymmetric quantities.
114  IF (is_asym) THEN
115  ALLOCATE(rmns_spline(mnmax,ns))
116  ALLOCATE(zmnc_spline(mnmax,ns))
117  ALLOCATE(lmnc_spline(mnmax,ns + 1))
118  ALLOCATE(currumns_spline(mnmax,ns))
119  ALLOCATE(currvmns_spline(mnmax,ns))
120  END IF
121 
122  END SUBROUTINE
123 
124 !-------------------------------------------------------------------------------
135 !-------------------------------------------------------------------------------
136  SUBROUTINE vmec_info_construct_island(mpol, ntor, ns, is_asym)
137 
138  IMPLICIT NONE
139 
140 ! Declare Arguments
141  INTEGER, INTENT(in) :: mpol
142  INTEGER, INTENT(in) :: ntor
143  INTEGER, INTENT(in) :: ns
144  LOGICAL, INTENT(in) :: is_asym
145 
146 ! Start of executable code
147 
148 ! Stellarator symmetric quantities.
149  ALLOCATE(rmnc_i(0:mpol,-ntor:ntor,ns))
150  ALLOCATE(zmns_i(0:mpol,-ntor:ntor,ns))
151  ALLOCATE(lmns_i(0:mpol,-ntor:ntor,ns + 1))
152  ALLOCATE(jcurrumnc(0:mpol,-ntor:ntor,ns))
153  ALLOCATE(jcurrvmnc(0:mpol,-ntor:ntor,ns))
154 
155  rmnc_i = 0
156  zmns_i = 0
157  lmns_i = 0
158  jcurrumnc = 0
159  jcurrvmnc = 0
160 
161 ! Asymmetric quantities.
162  IF (is_asym) THEN
163  ALLOCATE(rmns_i(0:mpol,-ntor:ntor,ns))
164  ALLOCATE(zmnc_i(0:mpol,-ntor:ntor,ns))
165  ALLOCATE(lmnc_i(0:mpol,-ntor:ntor,ns + 1))
166  ALLOCATE(jcurrumns(0:mpol,-ntor:ntor,ns))
167  ALLOCATE(jcurrvmns(0:mpol,-ntor:ntor,ns))
168 
169  rmns_i = 0
170  zmnc_i = 0
171  lmnc_i = 0
172  jcurrumns = 0
173  jcurrvmns = 0
174  END IF
175 
176  END SUBROUTINE
177 
178 !*******************************************************************************
179 ! DESTRUCTION SUBROUTINES
180 !*******************************************************************************
181 !-------------------------------------------------------------------------------
185 !-------------------------------------------------------------------------------
186  SUBROUTINE vmec_info_destruct_spline
187 
188  IMPLICIT NONE
189 
190 ! Start of executable code
191 
192 ! Stellarator symmetric quantities.
193  DEALLOCATE(rmnc_spline)
194  DEALLOCATE(zmns_spline)
195  DEALLOCATE(lmns_spline)
196  DEALLOCATE(currumnc_spline)
197  DEALLOCATE(currvmnc_spline)
198 
199 ! Asymmetric quantities.
200  IF (ALLOCATED(rmns_spline)) THEN
201  DEALLOCATE(rmns_spline)
202  DEALLOCATE(zmnc_spline)
203  DEALLOCATE(lmnc_spline)
204  DEALLOCATE(currumns_spline)
205  DEALLOCATE(currvmns_spline)
206  END IF
207 
208  END SUBROUTINE
209 
210 !-------------------------------------------------------------------------------
214 !-------------------------------------------------------------------------------
215  SUBROUTINE vmec_info_destruct_island
216 
217  IMPLICIT NONE
218 
219 ! Start of executable code
220 
221 ! Stellarator symmetric quantities.
222  DEALLOCATE(rmnc_i)
223  DEALLOCATE(zmns_i)
224  DEALLOCATE(lmns_i)
225  DEALLOCATE(jcurrumnc)
226  DEALLOCATE(jcurrvmnc)
227 
228 ! Asymmetric quantities.
229  IF (ALLOCATED(rmns_i)) THEN
230  DEALLOCATE(rmns_i)
231  DEALLOCATE(zmnc_i)
232  DEALLOCATE(lmnc_i)
233  DEALLOCATE(jcurrumns)
234  DEALLOCATE(jcurrvmns)
235  END IF
236 
237  END SUBROUTINE
238 
239 !*******************************************************************************
240 ! SETTER SUBROUTINES
241 !*******************************************************************************
242 !-------------------------------------------------------------------------------
250 !-------------------------------------------------------------------------------
251  SUBROUTINE vmec_info_set_wout(wout_file, ns_in, mpol_in, ntor_in, nfp_in)
252  USE descriptor_mod, ONLY: iam
253  USE v3_utilities, ONLY: assert_eq, assert
254  USE island_params
255  USE stel_constants, ONLY: twopi, zero, mu0
256  USE island_params, hs=>hs_i, ns=>ns_i
257  USE shared_data, ONLY: lasym, torflux, polflux
258 
259  IMPLICIT NONE
260 
261 ! Declare Arguments
262  CHARACTER(len=*), INTENT(in) :: wout_file
263  INTEGER, INTENT(IN) :: ns_in
264  INTEGER, INTENT(IN) :: mpol_in
265  INTEGER, INTENT(IN) :: ntor_in
266  INTEGER, INTENT(IN) :: nfp_in
267 
268 ! Local variables
269  INTEGER :: istat
270  INTEGER :: js
271  REAL (dp) :: temp
272 
273 ! Start of executable code
274 
275 ! Load wout file.
276  CALL read_wout_file(wout_file, istat)
277  CALL assert_eq(0, istat, 'Read-wout error in vmec_info_set_wout')
278 
279  IF (nfp_in .lt. 1) THEN
280  nfp_i = nfp_vmec
281  ELSE
282  nfp_i = nfp_in
283  END IF
284 
285  wb_i = (twopi*twopi)*wb_vmec
286  wp_i = (twopi*twopi)*wp_vmec
287  volume_i = volume_vmec
288  rmajor_i = raxis_vmec(0,1)
289 
290  CALL assert(wb_i .gt. zero, 'wb_vmec = 0!')
291 
292  ns = ns_in
293  nsh = ns_in - 1
294 
295  fourier_context => fourier_class(mpol_in, ntor_in, nu_i, nv_i, nfp_i, &
296  lasym)
297 
298 ! Spline r, z, and l Fourier components in s from original vmec mesh (s ~ phi,
299 ! ns_vmec points) to a "polar" mesh s ~ sqrt(phi) for better axis resolution.
300  CALL vmec_info_construct_spline(mnmax, ns, lasym)
301  CALL vmec_info_set_rzl_splines(rmnc_vmec, zmns_vmec, lmns_vmec, &
302  currumnc, currvmnc, &
303  rmnc_spline, zmns_spline, lmns_spline, &
304  currumnc_spline, currvmnc_spline, f_cos)
305  IF (lasym) THEN
306  CALL vmec_info_set_rzl_splines(rmns_vmec, zmnc_vmec, lmnc_vmec, &
307  currumns, currvmns, &
308  rmns_spline, zmnc_spline, lmnc_spline, &
309  currumns_spline, currvmns_spline, &
310  f_sin)
311  END IF
312 
313 ! Spline 1-D arrays: careful -> convert phipf VMEC and multiply by
314 ! ds-vmec/ds-island, since phipf_i = d(PHI)/ds-island
315  ALLOCATE(phipf_i(ns), chipf_i(ns), presf_i(ns), stat=istat)
316  CALL vmec_info_spline_oned_array(chipf_vmec, chipf_i, istat)
317  CALL vmec_info_spline_oned_array(phipf_vmec, phipf_i, istat)
318  presf_vmec = mu0*presf_vmec
319  CALL vmec_info_spline_oned_array(presf_vmec, presf_i, istat)
320 
321 ! Pessure should never be negative.
322  WHERE (presf_i .lt. 0)
323  presf_i = 0
324  END WHERE
325 
326 ! Scale phipf_i and convert to sqrt(flux) mesh by multiplying by 2*s
327  phipf_i = phipf_i/twopi
328  chipf_i = chipf_i/twopi
329 
330 ! Mapping s_vmec = (s_siesta)^2, so d(s_vmec)/d(s_siesta) = 2 s_siesta, where
331 ! s_siesta = hs_i*(js-1)
332  DO js = 1, ns
333  phipf_i(js) = 2.0*hs*(js - 1)*phipf_i(js)
334  chipf_i(js) = 2.0*hs*(js - 1)*chipf_i(js)
335  END DO
336 
337  IF (iam .EQ. 0) THEN
338  temp = twopi*hs*(sum(phipf_i(2:ns)) + sum(phipf_i(1:ns-1)))/2
339  WRITE (33, 1000) volume_i, 1.e-6_dp*(wb_i+wp_i/(gamma-1))/mu0, &
340  1.e-6_dp*wp_i/mu0, temp, wp_i/wb_i, &
341  chipf_i(2)/phipf_i(2), chipf_i(ns)/phipf_i(ns), &
342  gamma
343  END IF
344 
345  ALLOCATE(torflux(ns), polflux(ns))
346  torflux(1) = 0
347  polflux(1) = 0;
348  DO js = 2, ns
349  polflux(js) = polflux(js - 1) + hs*(chipf_i(js) + chipf_i(js - 1))/2
350  torflux(js) = torflux(js - 1) + hs*(phipf_i(js) + phipf_i(js - 1))/2
351  END DO
352 
353  CALL vmec_info_set_rzl
354 
355 ! Construct R
356 1000 FORMAT(/,' INITIAL PARAMETERS (FROM VMEC)',/, &
357  ' PHYSICS QUANTITIES',/, &
358  ' PLASMA VOLUME (M^3): ',1pe12.4, &
359  ' TOTAL MHD ENERGY (MJ): ', 1pe16.8,/, &
360  ' THERMAL ENERGY (MJ): ', 1pe12.4, &
361  ' EDGE TOROIDAL FLUX (Wb): ', 1pe12.4,/, &
362  ' DIMENSIONLESS QUANTITIES',/, &
363  ' <BETA>: ',1pe12.4, ' IOTA(0): ',1pe12.4, &
364  ' IOTA(1) ',1pe12.4, ' GAMMA: ',1pe12.4,/, 21('-'),/)
365 
366  END SUBROUTINE
367 
368 !-------------------------------------------------------------------------------
384 !-------------------------------------------------------------------------------
385  SUBROUTINE vmec_info_set_rzl_splines(rmn_vmec, zmn_vmec, lmn_vmec, &
386  currumn_vmec, currvmn_vmec, &
387  rmn_spline, zmn_spline, lmn_spline, &
388  currumn_spline, currvmn_spline, &
389  parity)
390  USE shared_data, ONLY: jsupvdota
391  USE stel_constants, ONLY: mu0
392  USE v3_utilities, ONLY: assert, assert_eq
393 
394  IMPLICIT NONE
395 
396 ! Declare Arguments
397  REAL (dp), DIMENSION(:,:), INTENT(in) :: rmn_vmec
398  REAL (dp), DIMENSION(:,:), INTENT(in) :: zmn_vmec
399  REAL (dp), DIMENSION(:,:), INTENT(in) :: lmn_vmec
400  REAL (dp), DIMENSION(:,:), INTENT(in) :: currumn_vmec
401  REAL (dp), DIMENSION(:,:), INTENT(in) :: currvmn_vmec
402  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(inout) :: rmn_spline
403  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(inout) :: zmn_spline
404  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(inout) :: lmn_spline
405  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(inout) :: currumn_spline
406  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(inout) :: currvmn_spline
407  INTEGER, INTENT(in) :: parity
408 
409 ! Local variables
410  INTEGER :: istat
411  INTEGER :: modes
412  INTEGER :: modes_nyq
413  INTEGER :: n_off
414  INTEGER :: n_max
415  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: tempmn1
416  REAL (dp), DIMENSION(:,:), ALLOCATABLE :: tempmn2
417 
418 ! Start of executable code
419  IF (parity .eq. f_sin .and. .not.lasym_vmec) THEN
420  rmn_spline = 0
421  zmn_spline = 0
422  lmn_spline = 0
423  currumn_spline = 0
424  currvmn_spline = 0
425  RETURN
426  END IF
427 
428  istat = 0
429 
430  CALL vmec_info_spline_modes(rmn_vmec, rmn_spline, xm_vmec, istat)
431  CALL assert_eq(0, istat, 'Error splining rmn')
432 
433  CALL vmec_info_spline_modes(zmn_vmec, zmn_spline, xm_vmec, istat)
434  CALL assert_eq(0, istat, 'Error splining zmn')
435 
436  ALLOCATE(tempmn1(mnmax,ns_vmec + 1), stat=istat)
437  CALL vmec_info_set_ghost_points(lmn_vmec, tempmn1)
438  CALL vmec_info_spline_modes(tempmn1, lmn_spline, xm_vmec, istat)
439  CALL assert_eq(0, istat, 'Error splining lmn')
440  DEALLOCATE(tempmn1)
441 
442 ! The currents use the _nyq spectrum so we need to truncae down to just the
443 ! non-nyquest modes.
444  ALLOCATE (tempmn1(mnmax,ns_vmec), &
445  tempmn2(mnmax,ns_vmec), stat=istat)
446  tempmn1 = 0
447  tempmn2 = 0
448 
449 ! Assuming here that _nyq modes will be larger than the regular modes.
450  CALL assert(mnmax_nyq .ge. mnmax, 'ERROR: More modes than _nyq modes')
451 
452  n_max = int(maxval(xn_vmec))
453  n_off = 2*int((maxval(xn_nyq) - maxval(xn_vmec))/nfp_vmec) + 1
454 
455  modes_nyq = 1
456  DO modes = 1, mnmax
457  CALL assert(xm_vmec(modes) .eq. xm_nyq(modes_nyq), 'm mode mismatch')
458  CALL assert(xn_vmec(modes) .eq. xn_nyq(modes_nyq), 'n mode mismatch')
459 
460  tempmn1(modes,:) = mu0*currumn_vmec(modes_nyq,:)
461  tempmn2(modes,:) = mu0*currvmn_vmec(modes_nyq,:)
462 
463  IF (int(xn_vmec(modes)) .eq. n_max) THEN
464  modes_nyq = modes_nyq + n_off
465  ELSE
466  modes_nyq = modes_nyq + 1
467  END IF
468  END DO
469  IF (parity .eq. f_cos) THEN
470  jsupvdota = mu0*(sum(tempmn2(1,2:ns_vmec - 1)) &
471  + (tempmn2(1, 1) + tempmn2(1, ns_vmec))/2)
472  jsupvdota = jsupvdota/(ns_vmec - 1)
473  END IF
474 
475  CALL vmec_info_spline_modes(tempmn1, currumn_spline, xm_vmec, istat)
476  CALL vmec_info_spline_modes(tempmn2, currvmn_spline, xm_vmec, istat)
477 
478  DEALLOCATE(tempmn1, tempmn2)
479 
480  END SUBROUTINE
481 
482 !-------------------------------------------------------------------------------
489 !-------------------------------------------------------------------------------
490  SUBROUTINE vmec_info_set_ghost_points(ymn_vmec, ymn_ghost)
491 
492  IMPLICIT NONE
493 
494 ! Declare Arguments
495  REAL (dp), DIMENSION(:,:), INTENT(in) :: ymn_vmec
496  REAL (dp), DIMENSION(:,:), INTENT(out) :: ymn_ghost
497 
498 ! Local variables
499  INTEGER :: modes
500 
501 ! Start of executable code
502  ymn_ghost(:,2:ns_vmec) = ymn_vmec(:,2:ns_vmec)
503  ymn_ghost(:,ns_vmec + 1) = 2*ymn_vmec(:,ns_vmec) &
504  & - ymn_vmec(:,ns_vmec - 1)
505 
506  WHERE (xm_vmec .eq. 0)
507  ymn_ghost(:,1) = 2*ymn_vmec(:,2) - ymn_vmec(:,3)
508  ELSE WHERE
509  ymn_ghost(:,1) = 0
510  END WHERE
511 
512  END SUBROUTINE
513 
514 !-------------------------------------------------------------------------------
520 !-------------------------------------------------------------------------------
521  SUBROUTINE vmec_info_set_rzl()
522  USE shared_data, ONLY: lasym
523  USE island_params, ONLY: mpol=>mpol_i, ntor=>ntor_i, ns=>ns_i
524 
525  IMPLICIT NONE
526 
527 ! Start executable code.
528  CALL vmec_info_construct_island(mpol, ntor, ns, lasym)
529  CALL vmec_info_repack(rmnc_i, zmns_i, lmns_i, jcurrumnc, jcurrvmnc, &
530  rmnc_spline, zmns_spline, lmns_spline, &
531  currumnc_spline, currvmnc_spline, f_cos)
532  IF (lasym) THEN
533  CALL vmec_info_repack(rmns_i, zmnc_i, lmnc_i, jcurrumns, jcurrvmns, &
534  rmns_spline, zmnc_spline, lmnc_spline, &
535  currumns_spline, currvmns_spline, f_sin)
536  END IF
537  CALL vmec_info_destruct_spline
538 
539  END SUBROUTINE
540 
541 !*******************************************************************************
542 ! Utilitiy SUBROUTINES
543 !*******************************************************************************
544 !-------------------------------------------------------------------------------
578 !-------------------------------------------------------------------------------
579  SUBROUTINE vmec_info_repack(rmn, zmn, lmn, currumn, currvmn, &
580  rmn_spl, zmn_spl, lmn_spl, &
581  currumn_spl, currvmn_spl, parity)
582  USE fourier, ONLY: m0, m1, n0, n1
583  USE island_params, ONLY: mpol=>mpol_i, ntor=>ntor_i, ns=>ns_i, &
584  & nfp=>nfp_i, hs=>hs_i, fourier_context
585  USE descriptor_mod, ONLY: iam
586  USE shared_data, ONLY: lverbose
587 
588  IMPLICIT NONE
589 
590 ! Declare Arguments
591  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: rmn
592  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: zmn
593  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: lmn
594  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: currumn
595  REAL (dp), ALLOCATABLE, DIMENSION(:,:,:), INTENT(inout) :: currvmn
596  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(in) :: rmn_spl
597  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(in) :: zmn_spl
598  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(in) :: lmn_spl
599  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(in) :: currumn_spl
600  REAL (dp), ALLOCATABLE, DIMENSION(:,:), INTENT(in) :: currvmn_spl
601  INTEGER, INTENT(in) :: parity
602 
603 ! Local variables
604  INTEGER :: mn_mode
605  INTEGER :: js
606  INTEGER :: m
607  INTEGER :: n
608  INTEGER :: n_abs
609  INTEGER :: s1
610  INTEGER :: s2
611  INTEGER :: istat
612 
613 ! Start executable code.
614  CALL assert_eq(mnmax, SIZE(rmn_spl,1),'rmn_spl wrong size1 in REPACK')
615  CALL assert_eq(ns, SIZE(rmn_spl,2),'rmn_spl wrong size2 in REPACK')
616  CALL assert_eq(ns + 1, SIZE(lmn_spl,2),'lmn_spl wrong size2 in REPACK')
617 
618  IF (parity .eq. f_cos) THEN
619  s1 = 1
620  ELSE
621  s2 = 1
622  END IF
623 
624  IF (iam .eq. 0 .and. lverbose .and. parity .eq. f_cos) THEN
625  m = int(maxval(xm_vmec))
626  n = int(maxval(abs(xn_vmec)))/nfp
627  IF (m .gt. mpol .or. n .gt. ntor) THEN
628  WRITE(*,1000) m, n
629  END IF
630  CALL assert_eq(0, mod(nfp_vmec, nfp), 'nfpin should be an even ' // &
631  'divisor of the VMEC nfp.')
632  END IF
633 
634 ! Flip vmec n to siesta -n so trig arg is in siesta notation.
635  DO mn_mode = 1,mnmax
636  m = xm_vmec(mn_mode)
637  n = xn_vmec(mn_mode)/nfp
638  n_abs = abs(n)
639  IF (m .le. mpol .and. n_abs .le. ntor) THEN
640 
641 ! For m=0 load only the positive n values.
642  IF (m .eq. 0) THEN
643  IF (parity .eq. f_cos) THEN
644  s2 = -sign(1,n)
645  END IF
646  IF (parity .eq. f_sin) THEN
647  s1 = -sign(1,n)
648  END IF
649  rmn(m,n_abs,:) = rmn(m,n_abs,:) + s1*rmn_spl(mn_mode,:)
650  zmn(m,n_abs,:) = zmn(m,n_abs,:) + s2*zmn_spl(mn_mode,:)
651  lmn(m,n_abs,:) = lmn(m,n_abs,:) + s2*lmn_spl(mn_mode,:)*iflipj
652  currumn(m,n_abs,:) = currumn(m,n_abs,:) &
653  + s1*currumn_spl(mn_mode,:)
654  currvmn(m,n_abs,:) = currvmn(m,n_abs,:) &
655  + s1*currvmn_spl(mn_mode,:)
656  ELSE
657  rmn(m,-n*iflipj,:) = rmn_spl(mn_mode,:)
658  zmn(m,-n*iflipj,:) = zmn_spl(mn_mode,:)*iflipj
659  lmn(m,-n*iflipj,:) = lmn_spl(mn_mode,:)
660  currumn(m,-n*iflipj,:) = currumn_spl(mn_mode,:)
661  currvmn(m,-n*iflipj,:) = currvmn_spl(mn_mode,:)
662  END IF
663  END IF
664  END DO
665 
666 ! Divide out the orthonorm factor used by the Fourier routines. mu0 is already
667 ! multipled into the current in vmec_info_set_RZL_splines.
668  DO js = 1, ns
669  lmn(:,:,js) = lmn(:,:,js)/fourier_context%orthonorm(0:mpol,-ntor:ntor)
670  rmn(:,:,js) = rmn(:,:,js)/fourier_context%orthonorm(0:mpol,-ntor:ntor)
671  zmn(:,:,js) = zmn(:,:,js)/fourier_context%orthonorm(0:mpol,-ntor:ntor)
672 ! IF (js .eq. 1) THEN
673 ! currumn(:,:,js) = currumn(:,:,js) &
674 ! / fourier_context%orthonorm(0:mpol,-ntor:ntor)
675 ! currvmn(:,:,js) = currvmn(:,:,js) &
676 ! / fourier_context%orthonorm(0:mpol,-ntor:ntor)
677 ! ELSE
678  currumn(:,:,js) = 2*hs*(js - 1)*currumn(:,:,js) &
679  / fourier_context%orthonorm(0:mpol,-ntor:ntor)
680  currvmn(:,:,js) = 2*hs*(js - 1)*currvmn(:,:,js) &
681  / fourier_context%orthonorm(0:mpol,-ntor:ntor)
682 ! END IF
683  END DO
684 
685  lmn(:,:,ns + 1) = lmn(:,:,ns + 1) &
686  / fourier_context%orthonorm(0:mpol,-ntor:ntor)
687 
688 ! lmn(m1:,:,1) = 0
689  rmn(m1:,:,1) = 0
690  zmn(m1:,:,1) = 0
691 ! currumn(:,:,1) = 0
692 ! currumn(m1:,:,1) = 0
693 ! currvmn(m1:,:,1) = 0
694 
695  IF (parity .eq. f_cos) THEN
696  lmn(m0,-ntor:n0,:) = 0
697  rmn(m0,-ntor:-n1,:) = 0
698  zmn(m0,-ntor:n0,:) = 0
699  currumn(m0,-ntor:-n1,:) = 0
700  currvmn(m0,-ntor:-n1,:) = 0
701  ELSE
702  lmn(m0,-ntor:-n0,:) = 0
703  rmn(m0,-ntor:n0,:) = 0
704  zmn(m0,-ntor:-n1,:) = 0
705  currumn(m0,-ntor:n0,:) = 0
706  currvmn(m0,-ntor:n0,:) = 0
707  END IF
708 
709  currumn(:,:,ns) = 0
710  currvmn(:,:,ns) = 0
711 
712 1000 FORMAT(' It is recommended to increase the number of modes to at' &
713  ' least VMEC values mpol=',i4,' ntor=',i4)
714 
715  END SUBROUTINE
716 
717 !-------------------------------------------------------------------------------
731 !-------------------------------------------------------------------------------
732  SUBROUTINE vmec_info_spline_modes(ymn_vmec, ymn_spline, xm_in, istat)
733  USE island_params, hs=>hs_i, ohs=>ohs_i, ns=>ns_i
734 
735  IMPLICIT NONE
736 
737 ! Declare Arguments
738  REAL (dp), DIMENSION(:,:), INTENT(in) :: ymn_vmec
739  REAL (dp), DIMENSION(:,:), INTENT(out) :: ymn_spline
740  REAL (dp), DIMENSION(:), INTENT(in) :: xm_in
741  INTEGER, INTENT(inout) :: istat
742 
743 ! Local variables
744  INTEGER :: js
745  INTEGER :: modes
746  INTEGER :: mp
747  REAL (dp), ALLOCATABLE, DIMENSION(:) :: snodes_vmec
748  REAL (dp), ALLOCATABLE, DIMENSION(:) :: y2_vmec
749  REAL (dp), ALLOCATABLE, DIMENSION(:) :: fac1
750  REAL (dp), ALLOCATABLE, DIMENSION(:) :: y_vmec
751  REAL (dp), ALLOCATABLE, DIMENSION(:) :: snodes
752  REAL (dp), ALLOCATABLE, DIMENSION(:) :: fac2
753  REAL (dp), ALLOCATABLE, DIMENSION(:) :: y_spline
754  REAL (dp) :: hsv
755  REAL (dp) :: expm
756  REAL (dp) :: offset
757 
758 ! Local parameters
759  REAL (dp), PARAMETER :: one = 1
760 
761 ! Start of executable code
762 
763 ! Check array sizes.
764  IF (SIZE(ymn_vmec, 2) .le. 1) THEN
765  istat = 1
766  RETURN
767  END IF
768 
769  IF (ns .LE. 1) THEN
770  istat = 2
771  RETURN
772  END IF
773 
774  ALLOCATE(snodes_vmec(SIZE(ymn_vmec, 2)))
775  ALLOCATE(y2_vmec(SIZE(ymn_vmec, 2)))
776  ALLOCATE(fac1(SIZE(ymn_vmec, 2)))
777  ALLOCATE(y_vmec(SIZE(ymn_vmec, 2)))
778 
779  ALLOCATE(snodes(SIZE(ymn_spline, 2)))
780  ALLOCATE(fac2(SIZE(ymn_spline, 2)))
781  ALLOCATE(y_spline(SIZE(ymn_spline, 2)))
782 
783  IF (SIZE(snodes_vmec) .eq. ns_vmec) THEN
784  offset = 1.0_dp
785  ELSE
786  offset = 1.5_dp
787  snodes_vmec(SIZE(snodes_vmec)) = 1
788  fac1(SIZE(snodes_vmec)) = 1;
789  fac2(SIZE(snodes)) = 1
790  snodes(SIZE(snodes)) = 1
791  END IF
792 
793 ! Set up "knots" on initial (vmec, svmec ~ phi) mesh and factor for taking out
794 ! (putting back) [sqrt(s)]**m factor for odd m.
795  hsv = one/(ns_vmec - 1)
796  snodes_vmec(1) = 0; fac1(1) = 1;
797  DO js = 2, ns_vmec
798  snodes_vmec(js) = hsv*(js - offset)
799  fac1(js) = one/sqrt(snodes_vmec(js)) !regularization factor: sqrt(FLUX) for odd modes
800  END DO
801  fac1(2:ns_vmec) = one/sqrt(snodes_vmec(2:ns_vmec))
802 
803 ! Set up s-nodes on final (snodes, splined) mesh [s_siesta(sfinal) ~ sfinal**2]
804 ! if s_siesta ~ sfinal, this is then the original vmec mesh.
805  ohs = ns - 1
806  hs = one/ohs
807  snodes(1) = 0
808  DO js = 2, ns
809  snodes(js) = hs*(js - offset)
810  END DO
811 
812  snodes = snodes*snodes !SIESTA s==fac2 ~ sqrt(s-vmec) mesh
813  fac2 = sqrt(snodes)
814 
815  DO modes = 1, SIZE(xm_in)
816  expm = 0 ! MRC: Check this. Orginal code looked like this would always
817  ! be non zero if the first loop iteration reset it.
818 
819  y_vmec = ymn_vmec(modes,:)
820  mp = xm_in(modes)
821 
822  IF (istat.EQ.0 .and. mp.GT.0) THEN
823  IF (mod(mp,2) .eq. 1) THEN
824  expm = 1
825  ELSE
826  expm = 2
827  END IF
828  y_vmec = y_vmec*(fac1**expm)
829  IF (mp .le. 2) THEN
830  y_vmec(1) = 2*y_vmec(2) - y_vmec(3)
831  END IF
832  END IF
833 
834 ! Initialize spline for each mode amplitude (factor out sqrt(s) factor for
835 ! odd-m) (spline routines in LIBSTELL Miscel folder).
836  CALL spline(snodes_vmec, y_vmec, SIZE(snodes_vmec), &
837  lower_b, upper_b, y2_vmec)
838 
839 ! Interpolate onto snodes mesh.
840  DO js = 1, SIZE(snodes)
841  CALL splint(snodes_vmec, y_vmec, y2_vmec, &
842  SIZE(snodes_vmec), snodes(js), y_spline(js))
843  END DO
844 
845  IF (istat .eq. 0 .and. mp .gt. 0 .and. expm .ne. 0) THEN
846  y_spline = y_spline*(fac2**expm)
847  END IF
848  ymn_spline(modes,:) = y_spline(:)
849  END DO
850 
851  istat = 0
852 
853  DEALLOCATE(snodes_vmec)
854  DEALLOCATE(y2_vmec)
855  DEALLOCATE(fac1)
856  DEALLOCATE(y_vmec)
857 
858  DEALLOCATE(snodes)
859  DEALLOCATE(fac2)
860  DEALLOCATE(y_spline)
861 
862  END SUBROUTINE
863 
864 !-------------------------------------------------------------------------------
873 !-------------------------------------------------------------------------------
874  SUBROUTINE vmec_info_spline_oned_array(y_vmec, y_spline, istat)
875  USE island_params, ONLY: hs=>hs_i, ns=>ns_i
876 
877  IMPLICIT NONE
878 
879 ! Declare Arguments
880  REAL (dp), DIMENSION(ns_vmec), INTENT(in) :: y_vmec
881  REAL (dp), DIMENSION(ns), INTENT(out) :: y_spline
882  INTEGER, INTENT(OUT) :: istat
883 
884 ! Local variables
885  INTEGER :: js
886  REAL (dp), DIMENSION(ns_vmec) :: snodes_vmec
887  REAL (dp), DIMENSION(ns_vmec) :: y2_vmec
888  REAL (dp), DIMENSION(ns) :: snodes
889  REAL (dp) :: hs_vmec
890 
891 ! Local parameters
892  REAL (dp), PARAMETER :: one = 1
893 
894 ! Start of executable code
895 
896 ! Check grid dimensions.
897  IF (ns_vmec .le. 1) THEN
898  istat = 1
899  RETURN
900  END IF
901 
902  IF (ns .LE. 1) THEN
903  istat = 2
904  RETURN
905  END IF
906 
907 ! Set up "knots" on initial (vmec, svmec ~ phi) mesh.
908  hs_vmec = one/(ns_vmec - 1)
909  DO js = 1, ns_vmec
910  snodes_vmec(js) = hs_vmec*(js - 1)
911  END DO
912 
913 ! Set up s-nodes on final (snodes, splined) mesh [s_siesta(sfinal) ~ sfinal**2]
914 ! if s_siesta ~ sfinal, this is then the original vmec mesh.
915 ! ohs = ns - 1
916 ! hs = one/ohs
917  DO js = 1, ns
918  snodes(js) = hs*(js - 1) !vmec mesh s-siesta~s-vmec
919  snodes(js) = snodes(js)*snodes(js) !polar mesh, s-siesta~s-vmec**2
920  END DO
921 
922 ! Initialize spline for each mode amplitude (factor out sqrt(s) factor for
923 ! odd-m)
924  CALL spline(snodes_vmec, y_vmec, ns_vmec, lower_b, upper_b, y2_vmec)
925 
926 ! Interpolate onto snodes mesh
927  DO js = 1, ns
928  CALL splint(snodes_vmec, y_vmec, y2_vmec, ns_vmec, &
929  snodes(js), y_spline(js))
930  END DO
931 
932  END SUBROUTINE
933 
934  END MODULE
shared_data::lverbose
logical lverbose
Use verbose screen output.
Definition: shared_data.f90:234
fourier::n1
integer, parameter n1
n = 1 mode.
Definition: fourier.f90:67
v3_utilities::assert_eq
Definition: v3_utilities.f:62
shared_data::torflux
real(dp), dimension(:), allocatable torflux
Toroidal flux profile.
Definition: shared_data.f90:169
island_params::phipf_i
real(dp), dimension(:), allocatable phipf_i
Radial toroidal flux derivative.
Definition: island_params.f90:68
island_params::fourier_context
type(fourier_class), pointer fourier_context
Fourier transform object.
Definition: island_params.f90:76
fourier
Module contains subroutines for computing FFT on parallel radial intervals. Converts quantities betwe...
Definition: fourier.f90:13
shared_data::jsupvdota
real(dp) jsupvdota
FIXME: UNKNOWN.
Definition: shared_data.f90:167
fourier::m0
integer, parameter m0
m = 0 mode.
Definition: fourier.f90:58
v3_utilities::assert
Definition: v3_utilities.f:55
fourier::f_sin
integer, parameter f_sin
Sine parity.
Definition: fourier.f90:27
island_params::wp_i
real(dp) wp_i
Thermal energy. 2Pi*wp read_wout_mod::wp.
Definition: island_params.f90:60
island_params::rmajor_i
real(dp) rmajor_i
Major radius.
Definition: island_params.f90:64
island_params
This file contains fix parameters related to the computational grids.
Definition: island_params.f90:10
island_params::wb_i
real(dp) wb_i
Magnetic energy. 2Pi*wb read_wout_mod::wb.
Definition: island_params.f90:58
shared_data::polflux
real(dp), dimension(:), allocatable polflux
Poloidal flux profile.
Definition: shared_data.f90:171
fourier::n0
integer, parameter n0
n = 0 mode.
Definition: fourier.f90:65
island_params::gamma
real(dp), parameter gamma
Adiabatic constant.
Definition: island_params.f90:66
shared_data::lasym
logical lasym
Use non-stellarator symmetry.
Definition: shared_data.f90:230
island_params::hs_i
real(dp) hs_i
Radial grid spacing. hs = s_i+1 - s_i.
Definition: island_params.f90:47
island_params::chipf_i
real(dp), dimension(:), allocatable chipf_i
Radial poloidal flux derivative.
Definition: island_params.f90:70
island_params::presf_i
real(dp), dimension(:), allocatable presf_i
Radial pressure. FIXME: Check if this is really needed.
Definition: island_params.f90:72
island_params::volume_i
real(dp) volume_i
Equilibrium volume read_wout_mod::volume.
Definition: island_params.f90:62
fourier::m1
integer, parameter m1
m = 1 mode.
Definition: fourier.f90:60
island_params::ohs_i
real(dp) ohs_i
Radial grid derivative factor. d X/ ds = (x_i+1 - x_i)/(s_i+1 - s_i) (1) Where ohs = 1/(s_i+1 - s_i) ...
Definition: island_params.f90:51
fourier::f_cos
integer, parameter f_cos
Cosine parity.
Definition: fourier.f90:25
shared_data
This file contains variables and parameters used by many modules in SIESTA.
Definition: shared_data.f90:10
fourier::fourier_class
Base class containing fourier memory.
Definition: fourier.f90:77
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17