-
Notifications
You must be signed in to change notification settings - Fork 0
/
neo_magfie_perturbation.f90
696 lines (663 loc) · 28.3 KB
/
neo_magfie_perturbation.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
MODULE neo_magfie_perturbation
!
! module containing numerical constants
USE neo_precision
! module containing switches from the input file (neo.in)
USE neo_control, ONLY: lab_swi, inp_swi
! interface to spline routines
USE inter_interfaces, ONLY: splinecof3_hi_driv,&
tf, tfp, tfpp, tfppp, splint_horner3
! get boozer_iota (from the axiymmetric file - is the same)
USE neo_magfie_mod, ONLY: boozer_iota
! magfie
USE magfie_mod, ONLY : magfie
! used for normalization (hat-quantities)
USE partpa_mod, ONLY : bmod0
!
IMPLICIT NONE
!
! define kind of double complex, quad, quad complex
PRIVATE dcp, qp, qcp
INTEGER, PARAMETER :: dcp=KIND(1.0_dp)
INTEGER, PARAMETER :: qp=SELECTED_REAL_KIND(33, 4931)
INTEGER, PARAMETER :: qcp=KIND(1.0_qp)
!
! internal (private) storage arrays for the
! Fourier spectra from the input file
PRIVATE mnmax_pert, ns_pert
INTEGER :: mnmax_pert, ns_pert
PROTECTED ixm_pert, ixn_pert
INTEGER, DIMENSION(:), ALLOCATABLE :: ixm_pert, ixn_pert
! arrays for s-values and Fourier modes
PRIVATE es_pert
PRIVATE rmnc_pert, zmnc_pert, lmnc_pert, bmnc_pert
REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: es_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: rmnc_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: zmnc_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: lmnc_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: bmnc_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: rmns_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: zmns_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: lmns_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: bmns_pert
!
! storage arrays for the splines of the 2d
! Fourier arrays
!PRIVATE a_rmnc_pert, b_rmnc_pert, c_rmnc_pert, d_rmnc_pert
!PRIVATE a_zmnc_pert, b_zmnc_pert, c_zmnc_pert, d_zmnc_pert
!PRIVATE a_lmnc_pert, b_lmnc_pert, c_lmnc_pert, d_lmnc_pert
PRIVATE a_bmnc_pert, b_bmnc_pert, c_bmnc_pert, d_bmnc_pert
!PRIVATE a_rmns_pert, b_rmns_pert, c_rmns_pert, d_rmns_pert
!PRIVATE a_zmns_pert, b_zmns_pert, c_zmns_pert, d_zmns_pert
!PRIVATE a_lmns_pert, b_lmns_pert, c_lmns_pert, d_lmns_pert
PRIVATE a_bmns_pert, b_bmns_pert, c_bmns_pert, d_bmns_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: a_rmnc_pert, b_rmnc_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: c_rmnc_pert, d_rmnc_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: a_zmnc_pert, b_zmnc_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: c_zmnc_pert, d_zmnc_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: a_lmnc_pert, b_lmnc_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: c_lmnc_pert, d_lmnc_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: a_bmnc_pert, b_bmnc_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: c_bmnc_pert, d_bmnc_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: a_rmns_pert, b_rmns_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: c_rmns_pert, d_rmns_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: a_zmns_pert, b_zmns_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: c_zmns_pert, d_zmns_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: a_lmns_pert, b_lmns_pert
!REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: c_lmns_pert, d_lmns_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: a_bmns_pert, b_bmns_pert
REAL(kind=dp), DIMENSION(:,:), ALLOCATABLE :: c_bmns_pert, d_bmns_pert
!
! storage arrays for r_mhalf and sp_index
! requested by the spline routines
PRIVATE r_m_pert, r_mhalf_pert, sp_index_pert
REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: r_m_pert, r_mhalf_pert
INTEGER(I4B), DIMENSION(:), ALLOCATABLE :: sp_index_pert
CHARACTER(len=1024), PUBLIC :: in_file_pert
INTEGER, PARAMETER :: m_phi_sgn = 1 ! positive or negative m_phi mode
INTEGER, PUBLIC :: m_phi
!
PUBLIC neo_magfie_pert
PRIVATE neo_magfie_pert_a, neo_magfie_pert_b
INTERFACE neo_magfie_pert
MODULE PROCEDURE neo_magfie_pert_a, neo_magfie_pert_b
END INTERFACE neo_magfie_pert
!
CONTAINS
SUBROUTINE neo_read_pert_control()
INTEGER :: r_un
character(1) :: dummy
r_un = 13123124
open(unit=r_un,file='neo_pert.in',status='old',form='formatted')
read (r_un,*) dummy
read (r_un,*) dummy
read (r_un,*) dummy
read (r_un,*) in_file_pert
close(unit=r_un)
END SUBROUTINE neo_read_pert_control
!
! Read Boozer files for the perturbation field
SUBROUTINE neo_read_pert()
!
! local definitions
!
! specify mode numbers and number of flux surfaces
INTEGER :: m0b_pert, n0b_pert
INTEGER :: m_max_pert, n_max_pert
INTEGER :: nfp_pert
! loop counters, status variables
INTEGER :: i,j
INTEGER :: i_alloc, r_un
! dummy variables
REAL(kind=dp) :: real_dummy
CHARACTER(5) :: char_dummy
!
! open Boozer file and read first quantities
r_un=27052014
in_file_pert=TRIM(ADJUSTL(in_file_pert))
OPEN(unit=r_un,file=in_file_pert,status='old',form='formatted')
!
IF (inp_swi .EQ. 8) THEN ! NEW IPP TOKAMAK
READ (r_un,*) char_dummy
READ (r_un,*) char_dummy
READ (r_un,*) char_dummy
READ (r_un,*) char_dummy
READ (r_un,*) char_dummy
READ (r_un,*) m0b_pert,n0b_pert,ns_pert,&
nfp_pert,real_dummy,real_dummy,real_dummy
m_max_pert = m0b_pert+1
n_max_pert = n0b_pert+1
mnmax_pert = m_max_pert*n_max_pert
!
! allocate storage arrays
ALLOCATE(ixm_pert(mnmax_pert), ixn_pert(mnmax_pert), stat = i_alloc)
IF(i_alloc /= 0) THEN
STOP "Allocation for the arrays containing the mode numbers&
& of the perturbation field failed!"
END IF
!
ALLOCATE(es_pert(ns_pert), stat = i_alloc)
IF(i_alloc /= 0) THEN
STOP "Allocation for the real array containing&
& the s-values of the perturbation field failed!"
END IF
!
ALLOCATE(rmnc_pert(ns_pert,mnmax_pert), zmnc_pert(ns_pert,mnmax_pert),&
lmnc_pert(ns_pert,mnmax_pert), bmnc_pert(ns_pert,mnmax_pert), stat = i_alloc)
IF(i_alloc /= 0) THEN
STOP "Allocation for the Fourier arrays for the perturbation field failed!"
END IF
!
! read input arrays
DO i =1, ns_pert
READ(r_un,*) char_dummy
READ(r_un,*) char_dummy
READ(r_un,*) es_pert(i),real_dummy,real_dummy,real_dummy,real_dummy,real_dummy
READ(r_un,*) char_dummy
DO j=1,mnmax_pert
!print *, 'j: ',j
READ(r_un,*) ixm_pert(j),ixn_pert(j),&
rmnc_pert(i,j),zmnc_pert(i,j),lmnc_pert(i,j),&
bmnc_pert(i,j)
!PRINT *, 'ixm,ixn,bmnc: ',ixm_pert(j),ixn_pert(j),bmnc_pert(i,j)
END DO
END DO
ELSEIF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
READ (r_un,*) char_dummy
READ (r_un,*) char_dummy
READ (r_un,*) char_dummy
READ (r_un,*) char_dummy
READ (r_un,*) char_dummy
READ (r_un,*) m0b_pert,n0b_pert,ns_pert,&
nfp_pert,real_dummy,real_dummy,real_dummy
m_max_pert = m0b_pert+1
n_max_pert = n0b_pert+1
mnmax_pert = m_max_pert*n_max_pert
!
! allocate storage arrays
ALLOCATE(ixm_pert(mnmax_pert), ixn_pert(mnmax_pert), stat = i_alloc)
IF(i_alloc /= 0) THEN
STOP "Allocation for the arrays containing the mode numbers&
& of the perturbation field failed!"
END IF
!
ALLOCATE(es_pert(ns_pert), stat = i_alloc)
IF(i_alloc /= 0) THEN
STOP "Allocation for the real array containing&
& the s-values of the perturbation field failed!"
END IF
!
ALLOCATE(rmnc_pert(ns_pert,mnmax_pert), zmnc_pert(ns_pert,mnmax_pert),&
lmnc_pert(ns_pert,mnmax_pert), bmnc_pert(ns_pert,mnmax_pert), stat = i_alloc)
IF(i_alloc /= 0) THEN
STOP "Allocation for the Fourier arrays for the perturbation field failed!"
END IF
ALLOCATE(rmns_pert(ns_pert,mnmax_pert), zmns_pert(ns_pert,mnmax_pert),&
lmns_pert(ns_pert,mnmax_pert), bmns_pert(ns_pert,mnmax_pert), stat = i_alloc)
IF(i_alloc /= 0) THEN
STOP "Allocation for the Fourier arrays for the perturbation field failed!"
END IF
!
! read input arrays
DO i =1, ns_pert
READ(r_un,*) char_dummy
READ(r_un,*) char_dummy
READ(r_un,*) es_pert(i),real_dummy,real_dummy,real_dummy,real_dummy,real_dummy
READ(r_un,*) char_dummy
DO j=1,mnmax_pert
!print *, 'j: ',j
READ(r_un,*) ixm_pert(j),ixn_pert(j),&
rmnc_pert(i,j),rmns_pert(i,j),zmnc_pert(i,j),zmns_pert(i,j),&
lmnc_pert(i,j),lmns_pert(i,j),bmnc_pert(i,j),bmns_pert(i,j)
!PRINT *, 'ixm,ixn,bmnc: ',ixm_pert(j),ixn_pert(j),bmnc_pert(i,j),bmns_pert(i,j),es_pert(i)
END DO
END DO
ELSE
PRINT *,'FATAL: There is yet no other input type for the perturbed field defined'
STOP
END IF
!
IF (lab_swi .EQ. 8) THEN ! NEW IPP TOKAMAK
ixn_pert = ixn_pert * nfp_pert
m_phi = ixn_pert(1)*m_phi_sgn
ixm_pert = ixm_pert
ELSEIF (lab_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
ixn_pert = ixn_pert * nfp_pert
m_phi = ixn_pert(1)*m_phi_sgn
ixm_pert = ixm_pert
ELSE
PRINT *,'FATAL: There is yet no other Laboratory for the perturbed field defined!'
STOP
END IF
! close Boozer file
CLOSE (unit=r_un)
RETURN
END SUBROUTINE neo_read_pert
!
! Initialization for splines along s
SUBROUTINE neo_init_spline_pert()
!
! local definitions
!
! loop variables, poloidal mode number
INTEGER :: i
INTEGER(I4B) :: m
! some maximum value of poloidal mode numbers
! (used for the computation of r_mhalf)
INTEGER, PARAMETER :: m_max_sp = 12
!
! testing
!
!!$ ! define switch and select a flux surface
!!$ INTEGER(I4B) :: swd
!!$ INTEGER :: k
!!$ REAL(dp) :: m0, f_es_pert, f_bmnc_pert, dummy
!
! allocate 2d arrays for the splines of the Fourier coefficients
!at the moment unused stuff
!ALLOCATE ( a_rmnc_pert(ns_pert,mnmax_pert), b_rmnc_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( c_rmnc_pert(ns_pert,mnmax_pert), d_rmnc_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( a_zmnc_pert(ns_pert,mnmax_pert), b_zmnc_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( c_zmnc_pert(ns_pert,mnmax_pert), d_zmnc_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( a_lmnc_pert(ns_pert,mnmax_pert), b_lmnc_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( c_lmnc_pert(ns_pert,mnmax_pert), d_lmnc_pert(ns_pert,mnmax_pert) )
ALLOCATE ( a_bmnc_pert(ns_pert,mnmax_pert), b_bmnc_pert(ns_pert,mnmax_pert) )
ALLOCATE ( c_bmnc_pert(ns_pert,mnmax_pert), d_bmnc_pert(ns_pert,mnmax_pert) )
! Additional data from Boozer files without Stellarator symmetry
IF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
!at the moment unused stuff
!ALLOCATE ( a_rmns_pert(ns_pert,mnmax_pert), b_rmns_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( c_rmns_pert(ns_pert,mnmax_pert), d_rmns_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( a_zmns_pert(ns_pert,mnmax_pert), b_zmns_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( c_zmns_pert(ns_pert,mnmax_pert), d_zmns_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( a_lmns_pert(ns_pert,mnmax_pert), b_lmns_pert(ns_pert,mnmax_pert) )
!ALLOCATE ( c_lmns_pert(ns_pert,mnmax_pert), d_lmns_pert(ns_pert,mnmax_pert) )
ALLOCATE ( a_bmns_pert(ns_pert,mnmax_pert), b_bmns_pert(ns_pert,mnmax_pert) )
ALLOCATE ( c_bmns_pert(ns_pert,mnmax_pert), d_bmns_pert(ns_pert,mnmax_pert) )
END IF
! allocate arrays requested by the spline routines
ALLOCATE ( r_m_pert(mnmax_pert), r_mhalf_pert(mnmax_pert) )
ALLOCATE ( sp_index_pert(ns_pert) )
!
! compute r_mhalf_pert and sp_index_pert requested by the spline routines
! (r_mhalf_pert enters the test function tf, which is used within the spline
! routines - tf(x) = x^r_mhalf_pert, x is here boozer_s. Since mode number m
! has to be symmetric around 0 for the representation of the perturbation
! field, r_mhalf_pert should be also symmetric -> m = abs(ixm_pert(i)).
! Without abs() the value of the test function grows unbounded for
! boozer_s << 1 and r_mhalf_pert < 0)
DO i = 1,mnmax_pert
! abs() - because of negative m-values
! (not necessary for Boozer files with m>=0, see neo_init_spline/neo_sub.f90)
m = ABS(ixm_pert(i))
IF (m .LE. m_max_sp) THEN
r_m_pert(i) = DBLE(m)
ELSE
IF (MODULO(m,2) .EQ. 1) THEN
r_m_pert(i) = DBLE(m_max_sp+1)
ELSE
r_m_pert(i) = DBLE(m_max_sp)
END IF
END IF
r_mhalf_pert(i) = r_m_pert(i) / 2._dp
END DO
sp_index_pert = (/ (i, i=1,ns_pert) /)
!
! 1-d splines of 2-d arrays
!at the moment unused stuff
!CALL splinecof3_hi_driv(es_pert, rmnc_pert, r_mhalf_pert,&
! a_rmnc_pert, b_rmnc_pert, c_rmnc_pert, d_rmnc_pert, sp_index_pert, tf)
!CALL splinecof3_hi_driv(es_pert, zmnc_pert, r_mhalf_pert,&
! a_zmnc_pert, b_zmnc_pert, c_zmnc_pert, d_zmnc_pert, sp_index_pert, tf)
!CALL splinecof3_hi_driv(es_pert, lmnc_pert, r_mhalf_pert,&
! a_lmnc_pert, b_lmnc_pert, c_lmnc_pert, d_lmnc_pert, sp_index_pert, tf)
CALL splinecof3_hi_driv(es_pert, bmnc_pert, r_mhalf_pert,&
a_bmnc_pert, b_bmnc_pert, c_bmnc_pert, d_bmnc_pert, sp_index_pert, tf)
! Additional data from Boozer files without Stellarator symmetry
IF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
!at the moment unused stuff
!CALL splinecof3_hi_driv(es_pert, rmns_pert, r_mhalf_pert,&
! a_rmns_pert, b_rmns_pert, c_rmns_pert, d_rmns_pert, sp_index_pert, tf)
!CALL splinecof3_hi_driv(es_pert, zmns_pert, r_mhalf_pert,&
! a_zmns_pert, b_zmns_pert, c_zmns_pert, d_zmns_pert, sp_index_pert, tf)
!CALL splinecof3_hi_driv(es_pert, lmns_pert, r_mhalf_pert,&
! a_lmns_pert, b_lmns_pert, c_lmns_pert, d_lmns_pert, sp_index_pert, tf)
CALL splinecof3_hi_driv(es_pert, bmns_pert, r_mhalf_pert,&
a_bmns_pert, b_bmns_pert, c_bmns_pert, d_bmns_pert, sp_index_pert, tf)
END IF
!
! Testing
!
!!$ swd = 0 ! no derivatives
!!$ k = 1
!!$ f_es_pert = es_pert(k)
!!$ DO i=1,mnmax_pert
!!$ m0 = r_mhalf_pert(i)
!!$ CALL splint_horner3(es_pert,a_bmnc_pert(:,i), b_bmnc_pert(:,i), c_bmnc_pert(:,i),&
!!$ d_bmnc_pert(:,i), swd, m0, f_es_pert, tf, tfp, tfpp, tfppp, f_bmnc_pert,&
!!$ dummy,dummy,dummy)
!!$ PRINT *, ixm_pert(i),ixn_pert(i),bmnc_pert(k,i),f_bmnc_pert
!!$ END DO
!
RETURN
END SUBROUTINE neo_init_spline_pert
!
! Compute the perturbation field for a certain x-value
SUBROUTINE neo_magfie_pert_a( x, bn_hat_pert )
!
! input / output
!
REAL(kind=dp), DIMENSION(:), INTENT(in) :: x
COMPLEX(kind=dcp), INTENT(out) :: bn_hat_pert
!
! local definitions
!
COMPLEX(kind=dcp), PARAMETER :: imun=(0.0_dp,1.0_dp)
INTEGER(i4b) :: swd
INTEGER :: i, m, n
REAL(kind=dp) :: yp, ypp, yppp
REAL(kind=dp) :: bmnc_pert_val, bmns_pert_val
COMPLEX(kind=dcp) :: expv
!
! read Boozer file and prepare spline routines (1st call)
!
IF (.NOT. ALLOCATED(es_pert)) THEN
CALL neo_read_pert()
CALL neo_init_spline_pert()
END IF
!
! direct summation of Fourier components
!
bn_hat_pert = (0.0_dp,0.0_dp)
!PRINT *,'mnmax_pert: ', mnmax_pert
DO i = 1, mnmax_pert
swd = 1
CALL splint_horner3(es_pert, a_bmnc_pert(:,i), b_bmnc_pert(:,i),&
c_bmnc_pert(:,i), d_bmnc_pert(:,i), swd, r_mhalf_pert(i),&
x(1), tf, tfp, tfpp, tfppp, bmnc_pert_val, yp, ypp, yppp)
! Additional data from Boozer files without Stellarator symmetry
IF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
CALL splint_horner3(es_pert, a_bmns_pert(:,i), b_bmns_pert(:,i),&
c_bmns_pert(:,i), d_bmns_pert(:,i), swd, r_mhalf_pert(i),&
x(1), tf, tfp, tfpp, tfppp, bmns_pert_val, yp, ypp, yppp)
END IF
IF (inp_swi .EQ. 8) THEN ! NEW IPP TOKAMAK
! requested representation of the perturbation field
! $B_n = \sum_{m>-\infty} \tilde{b}_{mn} \exp{i(m\vartheta_B+n\varphi_B)}$
n = ixn_pert(i)
!PRINT *,'n: ',n
! minus sign due to spectrum of the form
! $B = \sum_{mn} b_{mn} cos(m\vartheta_B-n\varphi_B)$
m = (-1)*ixm_pert(i)
!PRINT *,'m: ',m
! perturbation field is represented in terms of an expansion
! over the periodic toroidal Boozer angle times a phase
! (Note that x(3) is here a periodic angle - see mag_interface.f90 -
! while x(3) is the coordinate along the field line! Make sure that
! the different periodicities are treated correctly)
if (m_phi_sgn > 0) then
expv = EXP(imun*(m*x(3)+n*x(2)))
!PRINT *,'exponential (expv): ',expv
! part corresponding to the expansion over the periodic angle
bn_hat_pert = bn_hat_pert + bmnc_pert_val * expv
else
expv = EXP(-imun*(m*x(3)+n*x(2)))
!PRINT *,'exponential (expv): ',expv
! part corresponding to the expansion over the periodic angle
bn_hat_pert = bn_hat_pert + bmnc_pert_val * expv
end if
! alternative definition (for large aspect ratios, there is a factor 2 difference)
!m = ixm_pert(i)
!bn_hat_pert = bn_hat_pert + bmnc_pert_val * COS((m*boozer_iota-n)*x(2)) - &
! imun * bmnc_pert_val * SIN((m*boozer_iota-n)*x(2))
ELSEIF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
! requested representation of the perturbation field
! $B_n = \sum_{m>-\infty} \tilde{b}_{mn} \exp{i(m\vartheta_B+n\varphi_B)}$
n = ixn_pert(i)
m = ixm_pert(i)
!PRINT *,'m: ',m
! perturbation field is represented in terms of an expansion
! over the periodic toroidal Boozer angle times a phase
! (Note that x(3) is here a periodic angle - see mag_interface.f90 -
! while x(3) is the coordinate along the field line! Make sure that
! the different periodicities are treated correctly)
if (m_phi_sgn > 0) then
expv = EXP(imun*(m*x(3)+n*x(2)))
!PRINT *,'exponential (expv): ',expv
! part corresponding to the expansion over the periodic angle
bn_hat_pert = bn_hat_pert + (bmnc_pert_val - imun * bmns_pert_val) * expv
else
expv = EXP(-imun*(m*x(3)+n*x(2)))
!PRINT *,'exponential (expv): ',expv
! part corresponding to the expansion over the periodic angle
bn_hat_pert = bn_hat_pert + (bmnc_pert_val + imun * bmns_pert_val) * expv
end if
! alternative definition
!bn_hat_pert = bn_hat_pert + bmnc_pert_val * COS((m*boozer_iota+n)*x(2)) + &
! bmns_pert_val * SIN((m*boozer_iota+n)*x(2)) + imun * ( bmnc_pert_val * &
! SIN((m*boozer_iota+n)*x(2)) - bmns_pert_val * COS((m*boozer_iota+n)*x(2)) )
END IF
END DO
bn_hat_pert = bn_hat_pert / bmod0
!PRINT *,'bhat, bn_hat_pert: ',(bn_hat_pert/bhat),(1.0d-3*EXP(imun*m_phi*x(2)))
!PRINT *,'bn_hat_pert: ',bn_hat_pert
!PRINT *,'bmod0: ',bmod0
END SUBROUTINE neo_magfie_pert_a
!
! Compute the perturbation field for a certain x-value and its
! derivative over theta
SUBROUTINE neo_magfie_pert_b( x, bn_hat_pert, dbn_hat_pert_dtheta )
!
! input / output
!
REAL(kind=dp), DIMENSION(:), INTENT(in) :: x
COMPLEX(kind=dcp), INTENT(out) :: bn_hat_pert
COMPLEX(kind=dcp), INTENT(out) :: dbn_hat_pert_dtheta
!
! local definitions
!
COMPLEX(kind=dcp), PARAMETER :: imun=(0.0_dp,1.0_dp)
INTEGER(i4b) :: swd
INTEGER :: i, m, n
REAL(kind=dp) :: yp, ypp, yppp
REAL(kind=dp) :: bmnc_pert_val, bmns_pert_val
COMPLEX(kind=dcp) :: expv
!
! read Boozer file and prepare spline routines (1st call)
!
IF (.NOT. ALLOCATED(es_pert)) THEN
CALL neo_read_pert()
CALL neo_init_spline_pert()
END IF
!
! direct summation of Fourier components
!
bn_hat_pert = (0.0_dp,0.0_dp)
dbn_hat_pert_dtheta = (0.0_dp,0.0_dp)
!PRINT *,'mnmax_pert: ', mnmax_pert
DO i = 1, mnmax_pert
swd = 1
CALL splint_horner3(es_pert, a_bmnc_pert(:,i), b_bmnc_pert(:,i),&
c_bmnc_pert(:,i), d_bmnc_pert(:,i), swd, r_mhalf_pert(i),&
x(1), tf, tfp, tfpp, tfppp, bmnc_pert_val, yp, ypp, yppp)
! Additional data from Boozer files without Stellarator symmetry
IF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
CALL splint_horner3(es_pert, a_bmns_pert(:,i), b_bmns_pert(:,i),&
c_bmns_pert(:,i), d_bmns_pert(:,i), swd, r_mhalf_pert(i),&
x(1), tf, tfp, tfpp, tfppp, bmns_pert_val, yp, ypp, yppp)
END IF
IF (inp_swi .EQ. 8) THEN ! NEW IPP TOKAMAK
! requested representation of the perturbation field
! $B_n = \sum_{m>-\infty} \tilde{b}_{mn} \exp{i(m\vartheta_B+n\varphi_B)}$
n = ixn_pert(i)
!PRINT *,'n: ',n
! minus sign due to spectrum of the form
! $B = \sum_{mn} b_{mn} cos(m\vartheta_B-n\varphi_B)$
m = (-1)*ixm_pert(i)
!PRINT *,'m: ',m
! perturbation field is represented in terms of an expansion
! over the periodic toroidal Boozer angle times a phase
! (Note that x(3) is here a periodic angle - see mag_interface.f90 -
! while x(3) is the coordinate along the field line! Make sure that
! the different periodicities are treated correctly)
if (m_phi_sgn > 0) then
expv = EXP(imun*(m*x(3)+n*x(2)))
!PRINT *,'exponential (expv): ',expv
! part corresponding to the expansion over the periodic angle
bn_hat_pert = bn_hat_pert + bmnc_pert_val * expv
dbn_hat_pert_dtheta = dbn_hat_pert_dtheta + imun * m * bmnc_pert_val * expv
else
expv = EXP(-imun*(m*x(3)+n*x(2)))
!PRINT *,'exponential (expv): ',expv
! part corresponding to the expansion over the periodic angle
bn_hat_pert = bn_hat_pert + bmnc_pert_val * expv
dbn_hat_pert_dtheta = dbn_hat_pert_dtheta - imun * m * bmnc_pert_val * expv
end if
ELSEIF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
! requested representation of the perturbation field
! $B_n = \sum_{m>-\infty} \tilde{b}_{mn} \exp{i(m\vartheta_B+n\varphi_B)}$
n = ixn_pert(i)
m = ixm_pert(i)
!PRINT *,'m: ',m
! perturbation field is represented in terms of an expansion
! over the periodic toroidal Boozer angle times a phase
! (Note that x(3) is here a periodic angle - see mag_interface.f90 -
! while x(3) is the coordinate along the field line! Make sure that
! the different periodicities are treated correctly)
if (m_phi_sgn > 0) then
expv = EXP(imun*(m*x(3)+n*x(2)))
!PRINT *,'exponential (expv): ',expv
! part corresponding to the expansion over the periodic angle
bn_hat_pert = bn_hat_pert + (bmnc_pert_val - imun * bmns_pert_val) * expv
!PRINT *,'m, n, expv, bn_hat_pert, bmnc, bmns: ',m,n,expv,bn_hat_pert,bmnc_pert_val,bmns_pert_val
dbn_hat_pert_dtheta = dbn_hat_pert_dtheta +&
imun * m * (bmnc_pert_val - imun * bmns_pert_val) * expv
else
expv = EXP(-imun*(m*x(3)+n*x(2)))
!PRINT *,'exponential (expv): ',expv
! part corresponding to the expansion over the periodic angle
bn_hat_pert = bn_hat_pert + (bmnc_pert_val + imun * bmns_pert_val) * expv
!PRINT *,'m, n, expv, bn_hat_pert, bmnc, bmns: ',m,n,expv,bn_hat_pert,bmnc_pert_val,bmns_pert_val
dbn_hat_pert_dtheta = dbn_hat_pert_dtheta -&
imun * m * (bmnc_pert_val + imun * bmns_pert_val) * expv
end if
END IF
END DO
bn_hat_pert = bn_hat_pert / bmod0
dbn_hat_pert_dtheta = dbn_hat_pert_dtheta / bmod0
!PRINT *,'bhat, bn_hat_pert: ',(bn_hat_pert/bhat),(1.0d-3*EXP(imun*m_phi*x(2)))
!PRINT *,'bn_hat_pert: ',bn_hat_pert
!PRINT *,'bmod0: ',bmod0
END SUBROUTINE neo_magfie_pert_b
!
!
! compute the flux surface average of $g_{\varphi\varphi}$
! for symmetry flux coordinates
SUBROUTINE calc_av_gphph(x_start,phi_arr,ibeg,iend,av_gphph)
!
! input / output
!
REAL(kind=dp), DIMENSION(:), INTENT(in) :: x_start
REAL(kind=dp), DIMENSION(ibeg:iend), INTENT(in) :: phi_arr
INTEGER, INTENT(in) :: ibeg, iend
REAL(kind=dp), INTENT(out) :: av_gphph
!
! local definitions:
!
INTEGER :: k, w_un
REAL(kind=dp) :: bmoda, sqrtg, R_val, Z_val, fac1, fac2
REAL(kind=dp), DIMENSION(3) :: x, bder, hcovar, hctrvr
REAL(kind=dp), DIMENSION(3) :: hcurl, bcovar_s_hat_der
REAL(kind=dp), DIMENSION(ibeg:iend) :: R_arr, bmod_arr
!
! compute R, bmod as a function of phi
w_un=14112014
OPEN(unit=w_un,file='fluxsurface.dat',status='replace',form='formatted')
DO k = ibeg,iend
x(1)=x_start(1)
x(2)=phi_arr(k)
x(3)=x_start(3)+(x(2)-x_start(2))*boozer_iota
CALL magfie( x, bmoda, sqrtg, bder, hcovar, hctrvr,&
hcurl, bcovar_s_hat_der, R_val, Z_val )
R_arr(k) = R_val
bmod_arr(k) = bmoda
WRITE(w_un,*) x(2), R_val, Z_val, bmoda
END DO
CLOSE(unit=w_un)
!
! evaluate flux surface average
fac1 = 0.0_dp
fac2 = 0.0_dp
DO k = ibeg,iend-1
fac1 = fac1 + (phi_arr(k+1)-phi_arr(k)) * &
( (R_arr(k+1)/bmod_arr(k+1))**2.0_dp + (R_arr(k)/bmod_arr(k))**2.0_dp )
fac2 = fac2 + (phi_arr(k+1)-phi_arr(k)) * &
( (1.0_dp/bmod_arr(k+1))**2.0_dp + (1.0_dp/bmod_arr(k))**2.0_dp )
END DO
fac1 = 0.5_dp * fac1
fac2 = 0.5_dp * fac2
av_gphph = fac1 / fac2
!PRINT *,'av_gphph: ',av_gphph
!
END SUBROUTINE calc_av_gphph
!
!
! Compute the perturbation field amplitude
! for constant toroidal modenumber n = ixm_pert(1)
SUBROUTINE neo_magfie_pert_amp( x, bn_hat_pert )
!
! input / output
!
REAL(kind=dp), DIMENSION(:), INTENT(in) :: x
COMPLEX(kind=dcp), INTENT(out) :: bn_hat_pert
!
! local definitions
!
COMPLEX(kind=dcp), PARAMETER :: imun=(0.0_dp,1.0_dp)
INTEGER(i4b) :: swd
INTEGER :: i, m
REAL(kind=dp) :: yp, ypp, yppp
REAL(kind=dp) :: bmnc_pert_val, bmns_pert_val
COMPLEX(kind=dcp) :: expv
!
! read Boozer file and prepare spline routines (1st call)
!
IF (.NOT. ALLOCATED(es_pert)) THEN
CALL neo_read_pert()
CALL neo_init_spline_pert()
END IF
bn_hat_pert = (0.0_dp,0.0_dp)
DO i = 1, mnmax_pert
swd = 1
CALL splint_horner3(es_pert, a_bmnc_pert(:,i), b_bmnc_pert(:,i),&
c_bmnc_pert(:,i), d_bmnc_pert(:,i), swd, r_mhalf_pert(i),&
x(1), tf, tfp, tfpp, tfppp, bmnc_pert_val, yp, ypp, yppp)
! Additional data from Boozer files without Stellarator symmetry
IF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
CALL splint_horner3(es_pert, a_bmns_pert(:,i), b_bmns_pert(:,i),&
c_bmns_pert(:,i), d_bmns_pert(:,i), swd, r_mhalf_pert(i),&
x(1), tf, tfp, tfpp, tfppp, bmns_pert_val, yp, ypp, yppp)
END IF
IF (inp_swi .EQ. 8) THEN ! NEW IPP TOKAMAK
m = (-1)*ixm_pert(i)
if (m_phi_sgn > 0) then
expv = EXP(imun*m*x(3))
bn_hat_pert = bn_hat_pert + bmnc_pert_val * expv
else
expv = EXP(-imun*m*x(3))
bn_hat_pert = bn_hat_pert + bmnc_pert_val * expv
end if
ELSEIF (inp_swi .EQ. 9) THEN ! ASDEX-U (E. Strumberger)
m = ixm_pert(i)
if (m_phi_sgn > 0) then
expv = EXP(imun*m*x(3))
bn_hat_pert = bn_hat_pert + (bmnc_pert_val-imun*bmns_pert_val)*expv
else
expv = EXP(-imun*m*x(3))
bn_hat_pert = bn_hat_pert + (bmnc_pert_val+imun*bmns_pert_val)*expv
end if
END IF
END DO
bn_hat_pert = bn_hat_pert / bmod0
END SUBROUTINE neo_magfie_pert_amp
END MODULE neo_magfie_perturbation