ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
equilibrium_spider.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
2 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3 SUBROUTINE equilibrium_spider (EQUILIBRIUM_IN, COREPROF_IN, EQUILIBRIUM_OUT, code_parameters)
4 
5 !-------------------------------------------------------!
6 ! This routine provides interface between ETS !
7 ! and SPIDER equilibrium solver. !
8 !-------------------------------------------------------!
9 ! Source: --- !
10 ! Developers: D.Kalupin, G.Pereverzev !
11 ! Kontacts: D.Kalupin@fz-juelich.de !
12 ! !
13 ! Comments: created for ETS V&V !
14 ! !
15 ! for questions on SPIDER contact: !
16 ! Sergei Medvedev !
17 ! (medvedev@a5.kiam.ru) !
18 ! Andrei Ivanov !
19 ! (aai@a5.kiam.ru) !
20 ! !
21 !-------------------------------------------------------!
22 
23 
24  USE euitm_schemas
25  USE copy_structures
26 
27  USE itm_constants
28  USE euitm_xml_parser
29  USE deallocate_structures
30 
31  USE mod_spider
32  USE integration
33 
34  IMPLICIT NONE
35 
36 
37 ! +++ CPO derived types:
38  TYPE (type_equilibrium), POINTER :: equilibrium_in(:) !input CPO with geometry quantities from previous time
39  TYPE (type_equilibrium), POINTER :: equilibrium_out(:) !output CPO with geometry quantities from previous iteration
40  TYPE (type_coreprof), POINTER :: coreprof_in(:) !input CPO with plasma profiles
41 
42  TYPE (type_param) :: code_parameters
43 
44 
45 
46 ! +++ Dimensions:
47  INTEGER :: npsi !number of PSI points
48  INTEGER :: ndim1, ndim2 !number of points for 2-D parameters
49  INTEGER :: idim1, idim2 !index of points for 2-D parameters
50  INTEGER :: max_npoints !index of boundary points
51 
52  INTEGER :: icontrol,iinput !control sitch for SPIDER configuration
53  INTEGER :: key_iter !control sitch for SPIDER iterrations
54  INTEGER :: nrho, ntheta, neq !number of points for SPIDER I/O
55  INTEGER :: nbnd !number of points for SPIDER I/O
56  INTEGER :: ntransp !number of points for COREPROF
57  INTEGER :: ieq !index of points for SPIDER I/O grid
58  INTEGER :: ibnd !index of boundary points
59  INTEGER :: itime !time step number for SPIDER
60  INTEGER :: iter !iterration number
61  INTEGER :: itheta, irho !
62  INTEGER, PARAMETER :: spid_out = 0
63 
64 
65 ! +++ Internal parameters(0-D parameters):
66  REAL (R8) :: r0, z0 !magnetic axis (on first iteration equal to geom.axis)
67  REAL (R8) :: b0 !magnetic field on axis
68  REAL (R8) :: pc, pc_int !total plasma current
69  REAL (R8) :: amin !minor radius
70  REAL (R8) :: rhon !RHO (boundary value)
71  REAL (R8) :: rhonnew !RHO (new boundary value)
72  REAL (R8), ALLOCATABLE :: btmin(:), btmax(:)
73  REAL (R8) :: elongation
74  REAL (R8) :: tria_upper
75  REAL (R8) :: tria_lower
76 
77 ! Input:
78  REAL (R8), ALLOCATABLE :: rho(:)
79  REAL (R8), ALLOCATABLE :: r_inboard(:), r_outboard(:)
80  REAL (R8), ALLOCATABLE :: rr(:), rr_in(:)
81  REAL (R8), ALLOCATABLE :: jpar(:), intjpar(:)
82  REAL (R8), ALLOCATABLE :: j_bs(:), j_cd(:)
83  REAL (R8), ALLOCATABLE :: qsf(:)
84  REAL (R8), ALLOCATABLE :: mu(:)
85  REAL (R8), ALLOCATABLE :: pr(:), dpr(:)
86  REAL (R8), ALLOCATABLE :: te(:)
87  REAL (R8), ALLOCATABLE :: psi(:)
88  REAL (R8), ALLOCATABLE :: sigma(:)
89  REAL (R8), ALLOCATABLE :: pprime(:)
90  REAL (R8), ALLOCATABLE :: ffprime(:)
91  REAL (R8), ALLOCATABLE :: rzbnd(:)
92  REAL (R8), ALLOCATABLE :: rout(:,:), zout(:,:)
93 
94 ! Output from SPIDER:
95  REAL (R8), ALLOCATABLE :: rhonew(:)
96  REAL (R8), ALLOCATABLE :: rhonew_s(:)
97  REAL (R8), ALLOCATABLE :: gm11(:), gm22(:), gm33(:)
98  REAL (R8), ALLOCATABLE :: b2b02(:)
99  REAL (R8), ALLOCATABLE :: bb0(:)
100  REAL (R8), ALLOCATABLE :: b02b2(:)
101  REAL (R8), ALLOCATABLE :: fdia(:)
102  REAL (R8), ALLOCATABLE :: jdia(:), jdias(:)
103  REAL (R8), ALLOCATABLE :: vol(:), vprime(:), vprimes(:), vprime_psi(:)
104  REAL (R8), ALLOCATABLE :: surf(:)
105  REAL (R8), ALLOCATABLE :: gradrho(:)
106  REAL (R8), ALLOCATABLE :: drhoda(:)
107  REAL (R8), ALLOCATABLE :: shafr_shift(:)
108  REAL (R8), ALLOCATABLE :: elong(:)
109  REAL (R8), ALLOCATABLE :: elong_u(:)
110  REAL (R8), ALLOCATABLE :: elong_l(:)
111  REAL (R8), ALLOCATABLE :: triang_l(:)
112  REAL (R8), ALLOCATABLE :: triang_u(:)
113  REAL (R8), ALLOCATABLE :: ftrap(:)
114 
115 ! Output to ETS:
116  REAL (R8), ALLOCATABLE :: g1(:), g2(:), g3(:), g4(:)
117  REAL (R8), ALLOCATABLE :: g5(:), g6(:), g7(:)
118  REAL (R8), ALLOCATABLE :: area(:)
119 
120 ! Convergence parameters:
121  REAL (R8) :: conv
122  REAL (R8) :: tolerance = 1.0e-4_r8
123  REAL (R8) :: err_vprime, err_g1, err_fdia
124  REAL (R8), ALLOCATABLE :: vprime_old(:), g1_old(:), fdia_old(:)
125 
126 ! Work:
127  REAL (R8) :: theta, tria
128  REAL (R8), ALLOCATABLE :: amid(:)
129  REAL (R8), ALLOCATABLE :: fun(:), intfun(:), fun_in(:)
130 
131  REAL (R8), ALLOCATABLE :: amid_2d(:), rho_2d(:)
132  REAL (R8), ALLOCATABLE :: shafr_shift_2d(:), elong_2d(:), elong_u2d(:), elong_l2d(:), triang_l2d(:), triang_u2d(:)
133  REAL (R8), ALLOCATABLE :: psi_2d(:), jpar_2d(:), pr_2d(:), g11_2d(:), g22_2d(:), g33_2d(:)
134 
135  INTEGER :: return_status
136  LOGICAL, SAVE :: first=.true.
137  INTEGER :: i,j
138 
139  CHARACTER*18 :: str
140 !-------------------------------------------------------!
141 
142 
143 
144 ! +++ Assign SPIDER parameters:
145  CALL assign_spider_parameters(code_parameters, return_status)
146 
147  IF (return_status /= 0) THEN
148  WRITE(*,*) 'ERROR: Could not assign SPIDER parametrs.'
149  END IF
150 
151 
152 
153 ! +++ Set dimensions and allocate parameters:
154  npsi = SIZE (equilibrium_in(1)%profiles_1d%rho_tor, dim=1)
155  ndim1 = SIZE (equilibrium_in(1)%coord_sys%position%R, dim=1)
156  ndim2 = SIZE (equilibrium_in(1)%coord_sys%position%R, dim=2)
157  max_npoints = SIZE (equilibrium_in(1)%eqgeometry%boundary(1)%r, dim=1)
158  ntransp = SIZE (coreprof_in(1)%rho_tor, dim=1)
159 
160 
161  nbnd = max_npoints
162 
163 
164 
165 
166 ! +++ Set time index for SPIDER:
167  itime = 0
168 
169 
170 
171 ! +++ Allocate internal parameters:
172  ALLOCATE (rho(npsi))
173  ALLOCATE (r_inboard(npsi))
174  ALLOCATE (r_outboard(npsi))
175  ALLOCATE (rr(npsi))
176  ALLOCATE (rr_in(npsi))
177  ALLOCATE (jpar(npsi))
178  ALLOCATE (intjpar(npsi))
179  ALLOCATE (j_bs(npsi))
180  ALLOCATE (j_cd(npsi))
181  ALLOCATE (qsf(npsi))
182  ALLOCATE (mu(npsi))
183  ALLOCATE (pr(npsi))
184  ALLOCATE (dpr(npsi))
185  ALLOCATE (te(npsi))
186  ALLOCATE (psi(npsi))
187  ALLOCATE (sigma(npsi))
188  ALLOCATE (pprime(npsi))
189  ALLOCATE (ffprime(npsi))
190  ALLOCATE (rzbnd(nbnd * 2))
191  ALLOCATE (rout(npsi,ntheta))
192  ALLOCATE (zout(npsi,ntheta))
193 
194  ALLOCATE (b2b02(npsi))
195  ALLOCATE (bb0(npsi))
196  ALLOCATE (b02b2(npsi))
197  ALLOCATE (rhonew(npsi))
198  ALLOCATE (rhonew_s(npsi))
199  ALLOCATE (gm11(npsi))
200  ALLOCATE (gm22(npsi))
201  ALLOCATE (gm33(npsi))
202  ALLOCATE (fdia(npsi))
203  ALLOCATE (jdia(npsi))
204  ALLOCATE (jdias(npsi))
205  ALLOCATE (btmin(npsi))
206  ALLOCATE (btmax(npsi))
207  ALLOCATE (vol(npsi))
208  ALLOCATE (vprime(npsi))
209  ALLOCATE (vprimes(npsi))
210  ALLOCATE (vprime_psi(npsi))
211  ALLOCATE (surf(npsi))
212  ALLOCATE (gradrho(npsi))
213  ALLOCATE (drhoda(npsi))
214  ALLOCATE (shafr_shift(npsi))
215  ALLOCATE (elong(npsi))
216  ALLOCATE (elong_u(npsi))
217  ALLOCATE (elong_l(npsi))
218  ALLOCATE (triang_l(npsi))
219  ALLOCATE (triang_u(npsi))
220  ALLOCATE (amid(npsi))
221  ALLOCATE (ftrap(npsi))
222 
223  ALLOCATE (g1(npsi))
224  ALLOCATE (g2(npsi))
225  ALLOCATE (g3(npsi))
226  ALLOCATE (g4(npsi))
227  ALLOCATE (g5(npsi))
228  ALLOCATE (g6(npsi))
229  ALLOCATE (g7(npsi))
230  ALLOCATE (area(npsi))
231 
232  ALLOCATE (fun(npsi))
233  ALLOCATE (fun_in(npsi))
234  ALLOCATE (intfun(npsi))
235 
236  ALLOCATE (vprime_old(npsi))
237  ALLOCATE (g1_old(npsi))
238  ALLOCATE (fdia_old(npsi))
239 
240  ALLOCATE (amid_2d(ndim1))
241  ALLOCATE (rho_2d(ndim1))
242  ALLOCATE (elong_2d(ndim1))
243  ALLOCATE (elong_u2d(ndim1))
244  ALLOCATE (elong_l2d(ndim1))
245  ALLOCATE (triang_l2d(ndim1))
246  ALLOCATE (triang_u2d(ndim1))
247  ALLOCATE (shafr_shift_2d(ndim1))
248  ALLOCATE (psi_2d(ndim1))
249  ALLOCATE (jpar_2d(ndim1))
250  ALLOCATE (pr_2d(ndim1))
251  ALLOCATE (g11_2d(ndim1))
252  ALLOCATE (g22_2d(ndim1))
253  ALLOCATE (g33_2d(ndim1))
254 
255  rho = 0.0_r8
256  r_inboard = 0.0_r8
257  r_outboard = 0.0_r8
258  rr = 0.0_r8
259  rr_in = 0.0_r8
260  jpar = 0.0_r8
261  intjpar = 0.0_r8
262  j_bs = 0.0_r8
263  j_cd = 0.0_r8
264  qsf = 0.0_r8
265  mu = 0.0_r8
266  pr = 0.0_r8
267  dpr = 0.0_r8
268  te = 0.0_r8
269  psi = 0.0_r8
270  sigma = 0.0_r8
271  pprime = 0.0_r8
272  ffprime = 0.0_r8
273  rzbnd = 0.0_r8
274  rout = 0.0_r8
275  zout = 0.0_r8
276 
277  b2b02 = 0.0_r8
278  bb0 = 0.0_r8
279  b02b2 = 0.0_r8
280  rhonew = 0.0_r8
281  rhonew_s = 0.0_r8
282  gm11 = 0.0_r8
283  gm22 = 0.0_r8
284  gm33 = 0.0_r8
285  fdia = 0.0_r8
286  jdia = 0.0_r8
287  btmin = 0.0_r8
288  btmax = 0.0_r8
289  vol = 0.0_r8
290  vprime = 0.0_r8
291  vprimes = 0.0_r8
292  vprime_psi = 0.0_r8
293  surf = 0.0_r8
294  gradrho = 0.0_r8
295  drhoda = 0.0_r8
296  shafr_shift = 0.0_r8
297  elong = 0.0_r8
298  elong_u = 0.0_r8
299  elong_l = 0.0_r8
300  triang_l = 0.0_r8
301  triang_u = 0.0_r8
302  amid = 0.0_r8
303  ftrap = 0.0_r8
304 
305  g1 = 0.0_r8
306  g2 = 0.0_r8
307  g3 = 0.0_r8
308  g4 = 0.0_r8
309  g5 = 0.0_r8
310  g6 = 0.0_r8
311  g7 = 0.0_r8
312  area = 0.0_r8
313 
314  fun = 0.0_r8
315 
316  vprime_old = 1.0_r8
317  g1_old = 1.0_r8
318  fdia_old = 1.0_r8
319 
320 
321 
322 
323 
324 ! +++ equidistant equilibrium grid:
325  rhon = equilibrium_in(1)%profiles_1d%rho_tor(npsi)
326  DO ieq=1,npsi
327  rho(ieq) = rhon / (npsi-1) * (ieq-1)
328  END DO
329 
330 
331 
332 ! +++ Equilibrium boundary:
333  DO ibnd=1,nbnd*2
334  IF (ibnd.LE.nbnd) THEN
335  rzbnd(ibnd) = equilibrium_in(1)%eqgeometry%boundary(1)%r(ibnd)
336  ELSE IF (ibnd.GT.nbnd) THEN
337  rzbnd(ibnd) = equilibrium_in(1)%eqgeometry%boundary(1)%z(ibnd-nbnd)
338  END IF
339  END DO
340 
341 
342 
343 
344 ! +++ Set parameters:
345  r0 = equilibrium_in(1)%global_param%mag_axis%position%r
346  z0 = equilibrium_in(1)%global_param%mag_axis%position%z
347  b0 = equilibrium_in(1)%global_param%toroid_field%b0
348  amin = equilibrium_in(1)%eqgeometry%a_minor
349  elongation = equilibrium_in(1)%eqgeometry%elongation
350  tria_upper = equilibrium_in(1)%eqgeometry%tria_upper
351  tria_lower = equilibrium_in(1)%eqgeometry%tria_lower
352 
353 
354 
355 ! From EQUILIBRIUM:
356  CALL l3interp(-equilibrium_in(1)%profiles_1d%psi, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
357  psi, rho, npsi)
358  CALL l3interp(equilibrium_in(1)%profiles_1d%jparallel, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
359  jpar, rho, npsi)
360  CALL l3interp(-equilibrium_in(1)%profiles_1d%q, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
361  qsf, rho, npsi)
362  CALL l3interp(equilibrium_in(1)%profiles_1d%pressure, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
363  pr, rho, npsi)
364  CALL l3interp(equilibrium_in(1)%profiles_1d%gm5 / b0**2, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
365  b2b02, rho, npsi)
366  CALL l3interp(equilibrium_in(1)%profiles_1d%F_dia/r0/b0, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
367  jdia, rho, npsi)
368  CALL l3interp(equilibrium_in(1)%profiles_1d%F_dia, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
369  fdia, rho, npsi)
370  CALL l3interp(equilibrium_in(1)%profiles_1d%volume, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
371  vol, rho, npsi)
372  CALL l3interp(equilibrium_in(1)%profiles_1d%gm1, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
373  g1, rho, npsi)
374  CALL l3interp(equilibrium_in(1)%profiles_1d%gm2, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
375  g2, rho, npsi)
376  CALL l3interp(equilibrium_in(1)%profiles_1d%gm3, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
377  g3, rho, npsi)
378 
379  CALL l3deriv(vol, rho, npsi, vprime, rho, npsi)
380 
381 
382 
383 
384 ! From COREPROF:
385  CALL l3interp(coreprof_in(1)%te%value, coreprof_in(1)%rho_tor, ntransp, &
386  te, rho, npsi)
387  IF (ASSOCIATED(coreprof_in(1)%psi%sigma_par%value)) THEN
388  CALL l3interp(coreprof_in(1)%psi%sigma_par%value, coreprof_in(1)%rho_tor, ntransp, &
389  sigma, rho, npsi)
390  ELSE
391  sigma = 1.0e7_r8
392  END IF
393  IF (ASSOCIATED(coreprof_in(1)%psi%jni%value)) THEN
394  CALL l3interp(coreprof_in(1)%psi%jni%value, coreprof_in(1)%rho_tor, ntransp, &
395  j_cd, rho, npsi)
396  ELSE
397  j_cd = 0.0e0_r8
398  END IF
399 
400 
401 
402 
403 ! +++ Calculate major radius, RR:
404  IF(ASSOCIATED(equilibrium_in(1)%profiles_1d%r_inboard).AND. &
405  ASSOCIATED(equilibrium_in(1)%profiles_1d%r_outboard)) THEN
406  rr_in = (equilibrium_in(1)%profiles_1d%r_inboard + equilibrium_in(1)%profiles_1d%r_outboard) / 2.0_r8
407  CALL l3interp(rr_in, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
408  rr, rho, npsi)
409  ELSE
410  rr = r0
411  END IF
412 
413 
414 ! +++ Set parameters in SPIDER units:
415  mu = 1.0_r8 / qsf
416  j_bs = 0.0_r8 * 1.0e-6_r8 ! [MA/m^2]
417  j_cd = j_cd * 1.0e-6_r8 ! [MA/m^2]
418  jpar = jpar * 1.0e-6_r8 ! [MA/m^2]
419  pr = pr * 1.0e-6_r8 ! [MPa]
420  sigma = sigma * 1.0e-6_r8 ! [1/ (mkOhm*m)]
421 
422 
423 
424 
425 ! +++ Iterrations on equilibrium and current profile:
426  iter = 0
427 
428 7 CONTINUE
429 
430  iter = iter +1
431 
432  IF (iter.GT.100) THEN
433  WRITE(*,*) '### ERROR! Equilibrium loop in EQUILIBRIUM_SPIDER not converging'
434  stop 'Error - no convergence in EQUILIBRIUM_SPIDER'
435  RETURN
436 
437  END IF
438 
439  vprime_old = vprime !copy of equilibrium quantities for iterations
440  g1_old = g1
441  fdia_old = fdia
442 
443 
444 
445 
446 ! +++ Calculate P' and FF':
447  IF (iinput.EQ.0) THEN
448  CALL l3deriv(pr, rho, npsi, dpr, rho, npsi)
449 
450  pprime = - rr * qsf / b0 / rho * dpr
451 
452  IF (rho(1).EQ.0.0_r8) pprime(1) = pprime(2)
453 
454  ffprime = jdia / b2b02 * (jpar * rr / r0 - pprime)
455 
456  ELSE IF (iinput.EQ.1) THEN
457 
458  CALL l3interp(equilibrium_in(1)%profiles_1d%pprime, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
459  pprime, rho, npsi)
460  CALL l3interp(equilibrium_in(1)%profiles_1d%ffprime, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
461  ffprime, rho, npsi)
462  END IF
463 
464 
465 ! +++ Check current:
466  fun = 0.0e0_r8
467  fun = vprime * jpar / 2.0e0_r8 / itm_pi * b0 / fdia**2
468  CALL integral_value(npsi, rho, fun, intjpar)
469  pc_int = intjpar(npsi) * fdia(npsi) / 1.0e6_r8
470 
471  IF(equilibrium_in(1)%global_param%i_plasma.GT.-1.0e+40_r8) THEN
472  pc = equilibrium_in(1)%global_param%i_plasma / 1.0e6_r8
473 
474  ELSE
475  pc = pc_int
476 
477  END IF
478 
479 
480 
481 
482 ! +++ Save 'spidat.out' file:
483  IF (spid_out.EQ.1) THEN
484  WRITE(*,*) ''
485  WRITE(*,*) '-------------------------------------------------------'
486  WRITE(*,*) 'Set Ip (MA) ',pc
487  WRITE(*,*) 'Calculated Ip (MA) ',pc_int
488  WRITE(*,*) 'R0, RR, B0 ', r0, rr(npsi), b0
489  WRITE(*,*) '-------------------------------------------------------'
490  WRITE(*,*) ''
491 
492 
493  OPEN (2,file='spidat.out')
494  WRITE(2,100) ' jNPSIL,jnteta,NBND'
495  WRITE(2,102) nrho,ntheta,nbnd
496  WRITE(2,100) ' (rzbnd(j),j=1,2*NBND)'
497  WRITE(2,101) (rzbnd(i),i=1,nbnd*2)
498  WRITE(2,100) ' na1'
499  WRITE(2,102) npsi
500  WRITE(2,100) ' (eqpf(j),j=1,na1)'
501  WRITE(2,101) (pprime(i),i=1,npsi)
502  WRITE(2,100) ' (eqff(j),j=1,na1)'
503  WRITE(2,101) (ffprime(i),i=1,npsi)
504  WRITE(2,100) ' (fp(j),j=1,na1)'
505  WRITE(2,101) (psi(i),i=1,npsi)
506  WRITE(2,100) ' (rho(j),j=1,na1)'
507  WRITE(2,101) (rho(i),i=1,npsi)
508  WRITE(2,100) ' ipl,rtor,btor,roc,jnstep,'
509  WRITE(2,103) pc,r0,b0,rhon,itime
510  WRITE(2,100) ' (fp(j),j=1,na1)'
511  WRITE(2,101) (psi(i),i=1,npsi)
512  WRITE(2,100) ' (g11(j),j=1,na1)'
513  WRITE(2,101) (gm11(i),i=1,npsi)
514  WRITE(2,100) ' (g22(j),j=1,na1)'
515  WRITE(2,101) (gm22(i),i=1,npsi)
516  WRITE(2,100) ' (g33(j),j=1,na1)'
517  WRITE(2,101) (gm33(i),i=1,npsi)
518  WRITE(2,100) ' (vr(j),j=1,na1)'
519  WRITE(2,101) (vprime(i),i=1,npsi)
520  WRITE(2,100) ' (vrs(j),j=1,na1)'
521  WRITE(2,101) (vprimes(i),i=1,npsi)
522  WRITE(2,100) ' (slat(j),j=1,na1)'
523  WRITE(2,101) (surf(i),i=1,npsi)
524  WRITE(2,100) ' (gradro(j),j=1,na1)'
525  WRITE(2,101) (gradrho(i),i=1,npsi)
526  WRITE(2,100) ' (mu(j),j=1,na1)'
527  WRITE(2,101) (mu(i),i=1,npsi)
528  WRITE(2,100) ' (ipol(j),j=1,na1)'
529  WRITE(2,101) (jdia(j),j=1,npsi)
530  WRITE(2,100) ' (cc(j),j=1,na1)'
531  WRITE(2,101) (sigma(j),j=1,npsi)
532  WRITE(2,100) ' (Te(j),j=1,na1)'
533  WRITE(2,101) (te(j),j=1,npsi)
534  WRITE(2,100) ' (cubs(j),j=1,na1)'
535  WRITE(2,101) (j_bs(j),j=1,npsi)
536  WRITE(2,100) ' (cd(j),j=1,na1)'
537  WRITE(2,101) (j_cd(j),j=1,npsi)
538  WRITE(2,100) ' (pres(j),j=1,na1)'
539  WRITE(2,101) (pr(j),j=1,npsi)
540  WRITE(2,100) ' (cu(j),j=1,na1)'
541  WRITE(2,101) (jpar(j),j=1,npsi)
542  CLOSE(2)
543 ! +++
544  100 FORMAT(a)
545  101 FORMAT(1p,3g24.15e3)
546  102 FORMAT(3i12)
547  103 FORMAT(1p,3g24.15e3/g24.15e3,i12)
548  END IF
549 
550 
551 
552 
553 ! +++ Call SPIDER equilibrium solver:
554  CALL int_spider( &
555 ! Control:
556  icontrol &! key_equil in SPIDER (control key for spider input)
557 ! ICONTROL = 0 - equilibrium with prescribed p' and ff' profiles
558 ! + toroidal plasma current, Ipl;
559 ! The following inputs are used:
560 ! fp - psi 1D array,fp size is na1,
561 ! eqpf - p'*rtor 1D array,yeqpf size is na1,
562 ! eqff - ff'/rtor 1D array,yeqff size is na1,
563 ! ipl - toroidal plasma current.
564 ! rho, cu, pres are not used.
565 ! ICONTROL =-10 - equilibrium with prescribed p and <(j,B)>;
566 ! The following inputs are used:
567 ! rtor, btor - btor is vacuum toroidal field at r=rtor
568 ! rho - toroidal flux 1D array, sqrt(Phi/btor) rho size is na1
569 ! cu - <(j,B)>/btor 1D array, cu size is na1
570 ! pres - pressure 1D array,pres size is na1
571 ! fp, eqpf, eqpp, ipl are not used.
572 ! Input:
573  ,nrho &! radial grid in SPIDER
574  ,ntheta &! poloidal grid in SPIDER
575  ,nbnd &! boundary points
576  ,rzbnd &! array(2*nbnd) of boundary points
577  ,npsi &! radial grid for input (dimension of all arrays below)
578  ,pprime &! p'
579  ,ffprime &! FF'
580  ,rho &! rho
581  ,pc &! Ip [MA]
582  ,r0 &! R0 [m]
583  ,z0 &! Z0 [m]
584  ,b0 &! b0 [T]
585  ,rhon &! rho_bnd [m]
586 ! Output:
587  ,rhonnew &! new rho_bnd [m]
588 ! Input:
589  ,itime &! time step number
590 ! Input (and output if key_dmf = -2, -3)
591  ,psi &! psi [Wb] !defined on the main grid
592 ! Output:
593  ,gm11 &! <|nabla_rho|^2> * V' !defined on intermediate grid
594  ,gm22 &! V' R0 / (4*pi^2*Jdia) * <|nabla_rho/r|^2> !defined on intermediate grid
595  ,gm33 &! <R0^2/r^2> !defined on the main grid
596  ,vprime &! V'=dV/drho on the na1 grid !defined on the main grid
597  ,vprimes &! V'=dV/drho on the shifted grid !defined on intermediate grid
598  ,surf &! magnetic surface area [m^2] !defined on intermediate grid
599  ,gradrho &! |nabla_rho| !defined on intermediate grid
600  ,mu &! mu = 1/q !defined on intermediate grid
601  ,jdia &! F_dia/R0/B0 !defined on the main grid
602  ,btmax &! max of Bt on magn.axis !defined on intermediate grid
603  ,btmin &! min of Bt on magn.axis !defined on intermediate grid
604  ,b2b02 &! <B**2/B0**2> !defined on intermediate grid
605  ,bb0 &! <B/B0> !defined on intermediate grid
606  ,b02b2 &! <B0**2/B**2> !defined on intermediate grid
607  ,drhoda &! d_rho/d_a - translation from a to rho grid
608  ,rout &! flux surface coordinate [m]
609  ,zout &! flux surface coordinate [m]
610  ,shafr_shift &!
611  ,elong &!
612  ,elong_u &!
613  ,elong_l &!
614  ,triang_l &!
615  ,triang_u &!
616  ,amid &!
617  ,ftrap &!
618 ! Input:
619  ,sigma &! sigma [10^6*Ohm^-1*m^-1] !defined on the main grid
620  ,te &! Te [eV] !defined on the main grid
621  ,j_bs &! j_bootstrap [MA/m^2] !defined on the main grid
622  ,j_cd &! j_CD [MA/m^2] !defined on the main grid
623 ! Optional input:
624  ,pr &! pressure [MJ/m^3] !defined on the main grid
625  ,jpar) ! j_parallel [MA/m^2] !defined on the main grid
626 
627 
628 
629 
630 ! +++ Output equilibrium grids from SPIDER output:
631  DO ieq=1,npsi-1
632  rhonew(ieq) = rho(ieq) !main grid
633  END DO
634  rhonew(npsi) = rhonnew !last point of main and intermediate grids coinside
635 
636  DO ieq=1,npsi-1
637  rhonew_s(ieq) = 0.5_r8 * (rhonew(ieq) + rhonew(ieq+1)) !intermediate grid
638  END DO
639  rhonew_s(npsi-1) = rhonew_s(npsi-2) + rhonew(2) - rhonew(1) !NPSI-1 point of main and intermediate grids coinside
640  rhonew_s(npsi) = rhonnew !last point of main and intermediate grids coinside
641 
642 
643 
644 
645 ! +++ Calculate quantities not included in SPIDER output:
646 ! Metric coefficients defined on the main grid:
647  g1 = gm33 / r0**2
648  fdia = jdia * r0 * b0
649 
650  CALL integral_value(npsi, rhonew, vprime, vol) !volume [m^3]
651 
652 
653 
654 
655 ! Metric coefficients defined on intermediate grid:
656  g2 = gm22 / vprimes * 4.0_r8 * itm_pi**2 / r0
657  g3 = gm11 / vprimes
658  g4 = b02b2 / b0**2
659  g5 = b2b02 * b0**2
660  g6 = gm11 / b0**2
661  g7 = gradrho
662 
663  area = vprimes * gradrho !poloidal cross section [m^2]
664 
665  g4(npsi) = g4(npsi-1) + (g4(npsi-1)-g4(npsi-2)) & !linear interpolation, because SPIDER outputs two
666  / (rhonew_s(npsi-1)-rhonew_s(npsi-2)) & !equal last points
667  * (rhonew_s(npsi)-rhonew_s(npsi-1))
668  g5(npsi) = g5(npsi-1) + (g5(neq-1)-g5(neq-2)) & !linear interpolation, because SPIDER outputs two
669  / (rhonew_s(neq-1)-rhonew_s(neq-2)) & !equal last points
670  * (rhonew_s(neq)-rhonew_s(neq-1))
671 
672  CALL l3interp(g2, rhonew_s, neq, fun, rhonew, neq)
673  g2 = fun * jdia
674  CALL l3interp(g3, rhonew_s, neq, fun, rhonew, neq)
675  g3 = fun
676  CALL l3interp(g4, rhonew_s, neq, fun, rhonew, neq)
677  g4 = fun
678  CALL l3interp(g5, rhonew_s, neq, fun, rhonew, neq)
679  g5 = fun
680  CALL l3interp(g6, rhonew_s, neq, fun, rhonew, neq)
681  g6 = fun
682  CALL l3interp(g7, rhonew_s, neq, fun, rhonew, neq)
683  g7 = fun
684  CALL l3interp(mu, rhonew_s, neq, fun, rhonew, neq)
685  mu = fun
686  CALL l3interp(area, rhonew_s, neq, fun, rhonew, neq)
687  area = fun
688  CALL l3interp(surf, rhonew_s, neq, fun, rhonew, neq)
689  surf = fun
690 
691  rho = rhonew
692  rhon = rhonnew
693 
694  amid(neq) = amin
695 
696 ! +++ recalculation of the safety factor:
697  qsf = 1.0_r8 / mu
698 
699 
700 
701 
702 
703 
704 ! +++ CONVERGENCY CHECK:
705  conv = 0.0_r8
706 
707  err_vprime = maxval(abs(1.0_r8 - vprime_old(2:neq-1) / vprime(2:neq-1)))
708  err_g1 = maxval(abs(1.0_r8 - g1_old(2:neq-1) / g1(2:neq-1)))
709  err_fdia = maxval(abs(1.0_r8 - fdia_old(2:neq-1) / fdia(2:neq-1)))
710 
711 
712  conv = max( err_vprime, err_g1, err_fdia)
713 
714  write(*,*) ' CONV SPIDER =', conv, iter
715  write(*,*) ' err_VPRIME =', err_vprime
716  write(*,*) ' err_G1 =', err_g1
717  write(*,*) ' err_FDIA =', err_fdia
718 
719 ! IF(CONV.GT.TOLERANCE.AND.KEY_ITER.NE.0) GOTO 7
720 
721 
722 
723 
724 ! +++ Save output to EQUILIBRIUM CPO:
725 ! +++ Allocate output CPO:
726 ! IF(.NOT.ASSOCIATED(EQUILIBRIUM_OUT)) THEN
727  WRITE(*,*) 'Allocating EQUILIBRIUM_OUT'
728  ALLOCATE(equilibrium_out(1))
729 ! END IF
730 
731 
732 
733 
734 ! +++ Fill new EQUILIBRIUM CPO with quantities calculated SPIDER solver:
735 
736 
737 
738 
739 ! 1-D profiles:
740  CALL l3deriv(vol, psi, neq, vprime_psi, psi, neq)
741 
742 
743 ! Constraint:
744  call deallocate_cpo(equilibrium_out(1)%eqconstraint)
745  CALL copy_cpo(equilibrium_in(1)%eqconstraint, equilibrium_out(1)%eqconstraint)
746 
747 
748 ! Geometry:
749  call deallocate_cpo(equilibrium_out(1)%eqgeometry)
750  CALL copy_cpo(equilibrium_in(1)%eqgeometry, equilibrium_out(1)%eqgeometry)
751 
752 
753 ! Global parameters:
754  call deallocate_cpo(equilibrium_out(1)%global_param)
755  CALL copy_cpo(equilibrium_in(1)%global_param, equilibrium_out(1)%global_param)
756 
757 
758 
759 
760 ! GLOBAL_PARAM
761 
762  equilibrium_out(1)%global_param%area = area(neq)
763  equilibrium_out(1)%global_param%volume = vol(neq)
764  equilibrium_out(1)%global_param%i_plasma = pc * 1.0e6_r8
765  equilibrium_out(1)%global_param%psi_ax = -psi(1)
766  equilibrium_out(1)%global_param%psi_bound = -psi(neq)
767 
768 
769 
770 ! PROFILES_1D
771 ! Allocate 1-D profiles:
772  ALLOCATE (equilibrium_out(1)%profiles_1d%psi(npsi))
773  ALLOCATE (equilibrium_out(1)%profiles_1d%phi(npsi))
774  ALLOCATE (equilibrium_out(1)%profiles_1d%pressure(npsi))
775  ALLOCATE (equilibrium_out(1)%profiles_1d%F_dia(npsi))
776  ALLOCATE (equilibrium_out(1)%profiles_1d%pprime(npsi))
777  ALLOCATE (equilibrium_out(1)%profiles_1d%ffprime(npsi))
778  ALLOCATE (equilibrium_out(1)%profiles_1d%jphi(npsi))
779  ALLOCATE (equilibrium_out(1)%profiles_1d%jparallel(npsi))
780  ALLOCATE (equilibrium_out(1)%profiles_1d%q(npsi))
781  ALLOCATE (equilibrium_out(1)%profiles_1d%r_inboard(npsi))
782  ALLOCATE (equilibrium_out(1)%profiles_1d%r_outboard(npsi))
783  ALLOCATE (equilibrium_out(1)%profiles_1d%rho_tor(npsi))
784  ALLOCATE (equilibrium_out(1)%profiles_1d%rho_vol(npsi))
785  ALLOCATE (equilibrium_out(1)%profiles_1d%beta_pol(npsi))
786  ALLOCATE (equilibrium_out(1)%profiles_1d%li(npsi))
787  ALLOCATE (equilibrium_out(1)%profiles_1d%elongation(npsi))
788  ALLOCATE (equilibrium_out(1)%profiles_1d%tria_upper(npsi))
789  ALLOCATE (equilibrium_out(1)%profiles_1d%tria_lower(npsi))
790  ALLOCATE (equilibrium_out(1)%profiles_1d%volume(npsi))
791  ALLOCATE (equilibrium_out(1)%profiles_1d%vprime(npsi))
792  ALLOCATE (equilibrium_out(1)%profiles_1d%area(npsi))
793  ALLOCATE (equilibrium_out(1)%profiles_1d%aprime(npsi))
794  ALLOCATE (equilibrium_out(1)%profiles_1d%surface(npsi))
795  ALLOCATE (equilibrium_out(1)%profiles_1d%ftrap(npsi))
796  ALLOCATE (equilibrium_out(1)%profiles_1d%dpsidrho_tor(npsi))
797  ALLOCATE (equilibrium_out(1)%profiles_1d%b_av(npsi))
798 
799  ALLOCATE (equilibrium_out(1)%profiles_1d%gm1(npsi))
800  ALLOCATE (equilibrium_out(1)%profiles_1d%gm2(npsi))
801  ALLOCATE (equilibrium_out(1)%profiles_1d%gm3(npsi))
802  ALLOCATE (equilibrium_out(1)%profiles_1d%gm4(npsi))
803  ALLOCATE (equilibrium_out(1)%profiles_1d%gm5(npsi))
804  ALLOCATE (equilibrium_out(1)%profiles_1d%gm6(npsi))
805  ALLOCATE (equilibrium_out(1)%profiles_1d%gm7(npsi))
806  ALLOCATE (equilibrium_out(1)%profiles_1d%gm8(npsi))
807  ALLOCATE (equilibrium_out(1)%profiles_1d%gm9(npsi))
808 
809 
810 
811  DEALLOCATE (rhonew)
812  DEALLOCATE (fun)
813  ALLOCATE (rhonew(npsi))
814  ALLOCATE (fun(npsi))
815 
816  DO i = 1, npsi
817  rhonew(i) = rho(neq)*(i-1)/(npsi-1)
818  END DO
819 
820 
821  equilibrium_out(1)%profiles_1d%rho_tor = rhonew
822 
823 
824  fun_in = sqrt(vol/vol(neq))
825  CALL l3interp(fun_in, rho, neq, fun, rhonew, npsi)
826  equilibrium_out(1)%profiles_1d%rho_vol = fun
827  CALL l3interp(g1, rho, neq, fun, rhonew, npsi)
828  equilibrium_out(1)%profiles_1d%gm1 = fun
829  CALL l3interp(g2, rho, neq, fun, rhonew, npsi)
830  equilibrium_out(1)%profiles_1d%gm2 = fun
831  CALL l3interp(g3, rho, neq, fun, rhonew, npsi)
832  equilibrium_out(1)%profiles_1d%gm3 = fun
833  CALL l3interp(g4, rho, neq, fun, rhonew, npsi)
834  equilibrium_out(1)%profiles_1d%gm4 = fun
835  CALL l3interp(g5, rho, neq, fun, rhonew, npsi)
836  equilibrium_out(1)%profiles_1d%gm5 = fun
837  CALL l3interp(g6, rho, neq, fun, rhonew, npsi)
838  equilibrium_out(1)%profiles_1d%gm6 = fun
839  CALL l3interp(g7, rho, neq, fun, rhonew, npsi)
840  equilibrium_out(1)%profiles_1d%gm7 = g7
841  fun_in = 1._r8/g1**0.5
842  CALL l3interp(fun_in, rho, neq, fun, rhonew, npsi)!this is approximation
843  equilibrium_out(1)%profiles_1d%gm8 = fun
844  fun_in = g1**0.5
845  CALL l3interp(fun_in, rho, neq, fun, rhonew, npsi)!this is approximation
846  equilibrium_out(1)%profiles_1d%gm9 = fun
847  CALL l3interp(fdia, rho, neq, fun, rhonew, npsi)
848  equilibrium_out(1)%profiles_1d%F_dia = fun
849  CALL l3interp(vol, rho, neq, fun, rhonew, npsi)
850  equilibrium_out(1)%profiles_1d%volume = fun
851  CALL l3interp(vprime_psi, rho, neq, fun, rhonew, npsi)
852  equilibrium_out(1)%profiles_1d%vprime = -fun
853  CALL l3interp(psi, rho, neq, fun, rhonew, npsi)
854  equilibrium_out(1)%profiles_1d%psi = -fun
855  fun_in = pr * 1.0e6_r8
856  CALL l3interp(fun_in, rho, neq, fun, rhonew, npsi)
857  equilibrium_out(1)%profiles_1d%pressure = fun
858  fun_in = jpar * 1.0e6_r8
859  CALL l3interp(fun_in, rho, neq, fun, rhonew, npsi)
860  equilibrium_out(1)%profiles_1d%jparallel = fun
861  CALL l3interp(qsf, rho, neq, fun, rhonew, npsi)
862  equilibrium_out(1)%profiles_1d%q = -fun
863  CALL l3interp(surf, rho, neq, fun, rhonew, npsi)
864  equilibrium_out(1)%profiles_1d%surface = fun
865  CALL l3interp(area, rho, neq, fun, rhonew, npsi)
866  equilibrium_out(1)%profiles_1d%area = fun
867  CALL l3interp(elong, rho, neq, fun, rhonew, npsi)
868  equilibrium_out(1)%profiles_1d%elongation = fun
869  CALL l3interp(triang_u, rho, neq, fun, rhonew, npsi)
870  equilibrium_out(1)%profiles_1d%tria_upper = fun
871  CALL l3interp(triang_l, rho, neq, fun, rhonew, npsi)
872  equilibrium_out(1)%profiles_1d%tria_lower = fun
873  CALL l3interp(amid, rho, neq, fun, rhonew, npsi)
874  equilibrium_out(1)%profiles_1d%r_inboard = r0 - fun
875  equilibrium_out(1)%profiles_1d%r_outboard = r0 + fun
876  CALL l3interp(ftrap, rho, neq, fun, rhonew, npsi)
877  equilibrium_out(1)%profiles_1d%ftrap = 1.0_r8 - fun
878  CALL l3deriv(psi, rho, neq, fun, rhonew, npsi)
879  equilibrium_out(1)%profiles_1d%dpsidrho_tor = -fun
880  CALL l3interp(pprime, rho, neq, fun, rhonew, npsi)
881  equilibrium_out(1)%profiles_1d%pprime = -fun
882  CALL l3interp(ffprime, rho, neq, fun, rhonew, npsi)
883  equilibrium_out(1)%profiles_1d%ffprime = fun
884  fun_in = bb0* b0
885  CALL l3interp(fun_in, rho, neq, fun, rhonew, npsi)
886  equilibrium_out(1)%profiles_1d%b_av = fun
887 
888 
889 ! Stabilization of current oscillations should be removed later
890  CALL l3interp(g2, rho, neq, fun, rhonew, npsi)
891  CALL l3interp(equilibrium_in(1)%profiles_1d%gm2, &
892  equilibrium_in(1)%profiles_1d%rho_tor, &
893  SIZE (equilibrium_in(1)%profiles_1d%rho_tor), &
894  equilibrium_out(1)%profiles_1d%gm2, &
895  equilibrium_out(1)%profiles_1d%rho_tor, &
896  SIZE (equilibrium_out(1)%profiles_1d%rho_tor))
897  equilibrium_out(1)%profiles_1d%gm2 = equilibrium_out(1)%profiles_1d%gm2 * 0.9_r8 + fun * 0.1_r8
898 
899 
900 
901 
902 
903 
904 ! COORD_SYS
905  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%grid%dim1)) DEALLOCATE(equilibrium_out(1)%coord_sys%grid%dim1)
906  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%grid%dim2)) DEALLOCATE(equilibrium_out(1)%coord_sys%grid%dim2)
907  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%position%R)) DEALLOCATE(equilibrium_out(1)%coord_sys%position%R)
908  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%position%Z)) DEALLOCATE(equilibrium_out(1)%coord_sys%position%Z)
909  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%g_11)) DEALLOCATE(equilibrium_out(1)%coord_sys%g_11)
910  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%g_22)) DEALLOCATE(equilibrium_out(1)%coord_sys%g_22)
911  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%g_33)) DEALLOCATE(equilibrium_out(1)%coord_sys%g_33)
912  ALLOCATE(equilibrium_out(1)%coord_sys%grid_type(4))
913  ALLOCATE(equilibrium_out(1)%coord_sys%grid%dim1(ndim1))
914  ALLOCATE(equilibrium_out(1)%coord_sys%grid%dim2(ndim2))
915  ALLOCATE(equilibrium_out(1)%coord_sys%position%R(ndim1,ndim2))
916  ALLOCATE(equilibrium_out(1)%coord_sys%position%Z(ndim1,ndim2))
917  ALLOCATE(equilibrium_out(1)%coord_sys%g_11(ndim1,ndim2))
918  ALLOCATE(equilibrium_out(1)%coord_sys%g_22(ndim1,ndim2))
919  ALLOCATE(equilibrium_out(1)%coord_sys%g_33(ndim1,ndim2))
920 
921  rho_loop2: DO i=1,ndim1
922  amid_2d(i) = amin/(ndim1-1)*(i-1)
923  END DO rho_loop2
924  CALL l3interp(rho, amid, neq, rho_2d, amid_2d, ndim1)
925  CALL l3interp(elong, amid, neq, elong_2d, amid_2d, ndim1)
926  CALL l3interp(elong_u, amid, neq, elong_u2d, amid_2d, ndim1)
927  CALL l3interp(elong_l, amid, neq, elong_l2d, amid_2d, ndim1)
928  CALL l3interp(triang_l, amid, neq, triang_l2d, amid_2d, ndim1)
929  CALL l3interp(triang_u, amid, neq, triang_u2d, amid_2d, ndim1)
930  CALL l3interp(shafr_shift, amid, neq, shafr_shift_2d, amid_2d, ndim1)
931  CALL l3interp(psi, amid, neq, psi_2d, amid_2d, ndim1)
932  fun_in = jpar * 1.0e6_r8
933  CALL l3interp(fun_in, amid, neq, jpar_2d, amid_2d, ndim1)
934  fun_in = pr * 1.0e6_r8
935  CALL l3interp(fun_in, amid, neq, pr_2d, amid_2d, ndim1)
936  CALL l3interp(gm11, amid, neq, g11_2d, amid_2d, ndim1)
937  CALL l3interp(gm22, amid, neq, g22_2d, amid_2d, ndim1)
938  CALL l3interp(gm33, amid, neq, g33_2d, amid_2d, ndim1)
939 
940 
941 
942  DO itheta=1, ndim2
943  theta = REAL(itheta-1,r8)/REAL(ndim2-1,r8)*2.0_r8*itm_pi
944 
945  DO irho=1, ndim1
946  equilibrium_out(1)%coord_sys%grid_type(1) = '2'
947  equilibrium_out(1)%coord_sys%grid_type(2) = 'inverse'
948  equilibrium_out(1)%coord_sys%grid_type(3) = '3'
949  equilibrium_out(1)%coord_sys%grid_type(4) = 'polar'
950  equilibrium_out(1)%coord_sys%grid%dim1(irho) = -psi_2d(irho)
951  equilibrium_out(1)%coord_sys%grid%dim2(itheta) = theta
952  IF (theta.LE.itm_pi) tria = triang_u2d(irho)
953  IF (theta.GT.itm_pi) tria = triang_l2d(irho)
954  IF (theta.LE.itm_pi) elongation = elong_u2d(irho)
955  IF (theta.GT.itm_pi) elongation = elong_l2d(irho)
956 
957  equilibrium_out(1)%coord_sys%position%R(irho, itheta) = r0 + shafr_shift_2d(irho) + amid_2d(irho) &
958  * (cos(theta)-tria*(sin(theta))**2)
959  equilibrium_out(1)%coord_sys%position%Z(irho, itheta) = z0 + amid_2d(irho) * elongation * sin(theta)
960  equilibrium_out(1)%coord_sys%g_11(irho, itheta) = g11_2d(irho)
961  equilibrium_out(1)%coord_sys%g_22(irho, itheta) = g22_2d(irho)
962  equilibrium_out(1)%coord_sys%g_33(irho, itheta) = g33_2d(irho)
963  END DO
964  END DO
965 
966 
967 
968 
969 
970 
971 ! Allocate 2-D profiles:
972  ALLOCATE (equilibrium_out(1)%profiles_2d(1))
973  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%grid_type(4))
974  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%grid%dim1(ndim1))
975  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%grid%dim2(ndim2))
976  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%R(ndim1, ndim2))
977  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%Z(ndim1, ndim2))
978  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%psi(ndim1, ndim2))
979  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%theta(ndim1, ndim2))
980  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%jphi(ndim1, ndim2))
981  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%jpar(ndim1, ndim2))
982  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%br(ndim1, ndim2))
983  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%bz(ndim1, ndim2))
984  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%bphi(ndim1, ndim2))
985  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%vphi(ndim1, ndim2))
986  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%vtheta(ndim1, ndim2))
987  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%rho_mass(ndim1, ndim2))
988  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%pressure(ndim1, ndim2))
989  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%temperature(ndim1, ndim2))
990 
991 
992  DO itheta=1, ndim2
993  theta = REAL(itheta-1,r8)/REAL(ndim2-1,r8)*2.0_r8*itm_pi
994 
995  DO irho=1, ndim1
996  equilibrium_out(1)%profiles_2d(1)%grid_type(1) = '2'
997  equilibrium_out(1)%profiles_2d(1)%grid_type(2) = 'inverse'
998  equilibrium_out(1)%profiles_2d(1)%grid_type(3) = '3'
999  equilibrium_out(1)%profiles_2d(1)%grid_type(4) = 'polar'
1000  equilibrium_out(1)%profiles_2d(1)%grid%dim1(irho) = -psi_2d(irho)
1001  equilibrium_out(1)%profiles_2d(1)%grid%dim2(itheta) = theta
1002  IF (theta.LE.itm_pi) tria = triang_u2d(irho)
1003  IF (theta.GT.itm_pi) tria = triang_l2d(irho)
1004  IF (theta.LE.itm_pi) elongation = elong_u2d(irho)
1005  IF (theta.GT.itm_pi) elongation = elong_l2d(irho)
1006 
1007  equilibrium_out(1)%profiles_2d(1)%R(irho, itheta) = r0 + shafr_shift_2d(irho) + amid_2d(irho) &
1008  * (cos(theta)-tria*(sin(theta))**2)
1009  equilibrium_out(1)%profiles_2d(1)%Z(irho, itheta) = z0 + amid_2d(irho) * elongation * sin(theta)
1010  equilibrium_out(1)%profiles_2d(1)%theta(irho, itheta) = theta
1011 
1012  equilibrium_out(1)%profiles_2d(1)%psi(irho, itheta) = -psi_2d(irho)
1013  equilibrium_out(1)%profiles_2d(1)%jpar(irho, itheta) = jpar_2d(irho)
1014  equilibrium_out(1)%profiles_2d(1)%pressure(irho, itheta) = pr_2d(irho)
1015  END DO
1016  END DO
1017 
1018 
1019 
1020 
1021 
1022 
1023 ! CODE_PARAM
1024  IF(.NOT.ASSOCIATED(equilibrium_out(1)%codeparam%codename)) THEN
1025  ALLOCATE(equilibrium_out(1)%codeparam%codename(1)) ! For a string of 132 characters max.
1026  equilibrium_out(1)%codeparam%codename(1) = 'SPIDER'
1027  END IF
1028  IF(.NOT.ASSOCIATED(equilibrium_out(1)%codeparam%codeversion)) THEN
1029  ALLOCATE(equilibrium_out(1)%codeparam%codeversion(1)) ! For a string of 132 characters max.
1030  equilibrium_out(1)%codeparam%codeversion(1) = '14.12.2012'
1031  END IF
1032  IF(.NOT.ASSOCIATED(equilibrium_out(1)%codeparam%parameters)) THEN
1033  ALLOCATE(equilibrium_out(1)%codeparam%parameters(1)) ! For a string of 132 characters max.
1034  equilibrium_out(1)%codeparam%parameters(1) = 'my_code_specific_parameters'
1035  END IF
1036  IF(.NOT.ASSOCIATED(equilibrium_out(1)%codeparam%output_diag)) THEN
1037  ALLOCATE(equilibrium_out(1)%codeparam%output_diag(1)) ! For a string of 132 characters max.
1038  equilibrium_out(1)%codeparam%output_diag(1) = 'my_output_diag'
1039  END IF
1040  equilibrium_out(1)%codeparam%output_flag = 0 ! Integer output flag, 0 means the run was successful and can be used
1041 
1042 
1043 
1044 
1045 
1046 
1047 ! DATA_INFO
1048  IF(.NOT.ASSOCIATED(equilibrium_out(1)%datainfo%dataprovider)) THEN
1049  ALLOCATE(equilibrium_out(1)%datainfo%dataprovider(1)) ! For a string of 132 characters max.
1050  END IF
1051  equilibrium_out(1)%datainfo%dataprovider(1) = 'Denis Kalupin'
1052  equilibrium_out(1)%datainfo%cocos = 13!!!!12
1053 
1054 
1055 
1056 
1057 
1058 ! +++ Deallocate parameters:
1059  DEALLOCATE (rho)
1060  DEALLOCATE (r_inboard)
1061  DEALLOCATE (r_outboard)
1062  DEALLOCATE (rr)
1063  DEALLOCATE (rr_in)
1064  DEALLOCATE (jpar)
1065  DEALLOCATE (intjpar)
1066  DEALLOCATE (j_bs)
1067  DEALLOCATE (j_cd)
1068  DEALLOCATE (qsf)
1069  DEALLOCATE (mu)
1070  DEALLOCATE (pr)
1071  DEALLOCATE (dpr)
1072  DEALLOCATE (te)
1073  DEALLOCATE (psi)
1074  DEALLOCATE (sigma)
1075  DEALLOCATE (pprime)
1076  DEALLOCATE (ffprime)
1077  DEALLOCATE (rzbnd)
1078  DEALLOCATE (rout)
1079  DEALLOCATE (zout)
1080 
1081  DEALLOCATE (b2b02)
1082  DEALLOCATE (bb0)
1083  DEALLOCATE (b02b2)
1084  DEALLOCATE (rhonew)
1085  DEALLOCATE (rhonew_s)
1086  DEALLOCATE (gm11)
1087  DEALLOCATE (gm22)
1088  DEALLOCATE (gm33)
1089  DEALLOCATE (fdia)
1090  DEALLOCATE (jdia)
1091  DEALLOCATE (jdias)
1092  DEALLOCATE (vol)
1093  DEALLOCATE (vprime)
1094  DEALLOCATE (vprimes)
1095  DEALLOCATE (vprime_psi)
1096  DEALLOCATE (surf)
1097  DEALLOCATE (gradrho)
1098  DEALLOCATE (drhoda)
1099  DEALLOCATE (shafr_shift)
1100  DEALLOCATE (elong)
1101  DEALLOCATE (elong_u)
1102  DEALLOCATE (elong_l)
1103  DEALLOCATE (triang_l)
1104  DEALLOCATE (triang_u)
1105  DEALLOCATE (amid)
1106  DEALLOCATE (ftrap)
1107 
1108  DEALLOCATE (g1)
1109  DEALLOCATE (g2)
1110  DEALLOCATE (g3)
1111  DEALLOCATE (g4)
1112  DEALLOCATE (g5)
1113  DEALLOCATE (g6)
1114  DEALLOCATE (g7)
1115  DEALLOCATE (area)
1116 
1117  DEALLOCATE (fun)
1118  DEALLOCATE (fun_in)
1119  DEALLOCATE (intfun)
1120 
1121  DEALLOCATE (vprime_old)
1122  DEALLOCATE (g1_old)
1123  DEALLOCATE (fdia_old)
1124 
1125  DEALLOCATE (amid_2d)
1126  DEALLOCATE (rho_2d)
1127  DEALLOCATE (elong_2d)
1128  DEALLOCATE (elong_u2d)
1129  DEALLOCATE (elong_l2d)
1130  DEALLOCATE (triang_l2d)
1131  DEALLOCATE (triang_u2d)
1132  DEALLOCATE (shafr_shift_2d)
1133  DEALLOCATE (psi_2d)
1134  DEALLOCATE (jpar_2d)
1135  DEALLOCATE (pr_2d)
1136  DEALLOCATE (g11_2d)
1137  DEALLOCATE (g22_2d)
1138  DEALLOCATE (g33_2d)
1139 
1140 
1141 
1142 !+++ Remove SPIDER files from the directory:
1143  str = 'rm -rf *.dat *.wr'
1144  CALL system(str)
1145 
1146  RETURN
1147 
1148 
1149 
1150 
1151 
1152 
1153 
1154 
1155 
1156 
1157 
1158  CONTAINS
1159 
1160 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1161 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1162  SUBROUTINE assign_spider_parameters(codeparameters, return_status)
1163 
1164  !-------------------------------------------------------!
1165  ! This subroutine calls the XML parser for !
1166  ! the SPIDER parameters and assign the !
1167  ! resulting values to the corresponding variables !
1168  !-------------------------------------------------------!
1169  ! Source: --- !
1170  ! Developers: D.Kalupin !
1171  ! Kontacts: Denis.Kalupin@efda.org !
1172  ! !
1173  ! Comments: created for V&V between ETS and !
1174  ! ASTRA !
1175  ! !
1176  !-------------------------------------------------------!
1177 
1178  USE itm_types
1179  USE euitm_schemas
1180  USE euitm_xml_parser
1181 
1182  IMPLICIT NONE
1183 
1184 
1185  TYPE(type_param) :: codeparameters
1186  TYPE(tree) :: parameter_list
1187  TYPE(element), POINTER :: temp_pointer
1188 
1189  INTEGER(ITM_I4) :: return_status
1190  INTEGER(ITM_I4) :: i, nparm, n_values
1191  INTEGER :: n_data
1192 
1193  CHARACTER(len = 132) :: cname
1194  !-------------------------------------------------------!
1195 
1196 
1197 
1198  return_status = 0 ! no error
1199 
1200 
1201 ! +++ parse xml-string codeparameters%parameters
1202  WRITE(6,*) 'Calling euitm_xml_parse'
1203  CALL euitm_xml_parse(codeparameters, nparm, parameter_list)
1204  WRITE(6,*) 'Called euitm_xml_parse'
1205 
1206 
1207 ! +++ assign variables
1208  temp_pointer => parameter_list%first
1209 
1210  outer: do
1211  cname = char2str(temp_pointer%cname) ! necessary for AIX
1212  select case (cname)
1213 
1214 
1215 ! +++ parameters overall
1216  case ("parameters")
1217  temp_pointer => temp_pointer%child
1218  cycle
1219 
1220 
1221 ! +++ SPIDER settings
1222  case ("SPIDER")
1223  temp_pointer => temp_pointer%child
1224  cycle
1225 
1226  case ("icontrol")
1227  if (allocated(temp_pointer%cvalue)) &
1228  call char2num(temp_pointer%cvalue, icontrol)
1229 
1230  case ("key_iter")
1231  if (allocated(temp_pointer%cvalue)) &
1232  call char2num(temp_pointer%cvalue, key_iter)
1233 
1234  case ("nrho")
1235  if (allocated(temp_pointer%cvalue)) &
1236  call char2num(temp_pointer%cvalue, nrho)
1237 
1238  case ("ntheta")
1239  if (allocated(temp_pointer%cvalue)) &
1240  call char2num(temp_pointer%cvalue, ntheta)
1241 
1242  case ("neq")
1243  if (allocated(temp_pointer%cvalue)) &
1244  call char2num(temp_pointer%cvalue, neq)
1245 
1246  case ("tolerance")
1247  if (allocated(temp_pointer%cvalue)) &
1248  call char2num(temp_pointer%cvalue, tolerance)
1249 
1250  case ("iinput")
1251  if (allocated(temp_pointer%cvalue)) &
1252  call char2num(temp_pointer%cvalue, iinput)
1253 
1254  case default
1255  write(*, *) 'ERROR: invalid parameter', cname
1256  return_status = 1
1257  exit
1258  end select
1259 
1260 
1261  do
1262  if (associated(temp_pointer%sibling)) then
1263  temp_pointer => temp_pointer%sibling
1264  exit
1265  end if
1266  if (associated(temp_pointer%parent, parameter_list%first )) &
1267  exit outer
1268  if (associated(temp_pointer%parent)) then
1269  temp_pointer => temp_pointer%parent
1270  else
1271  write(*, *) 'ERROR: broken list.'
1272  return
1273  end if
1274  end do
1275  end do outer
1276 
1277 ! +++ destroy tree
1278  CALL destroy_xml_tree(parameter_list)
1279 
1280 
1281  RETURN
1282 
1283  END SUBROUTINE assign_spider_parameters
1284 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1285 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1286 
1287 
1288 END SUBROUTINE equilibrium_spider
1289 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1290 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1291 
1292 
1293 
1294 
1295 
l3deriv
subroutine l3deriv(y_in, x_in, nr_in, dydx_out, x_out, nr_out)
Definition: l3interp.f90:59
mod_spider
Definition: mod_spider.f90:1
fun
subroutine fun(X, F)
Definition: Ev2.f:10
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
r_inboard
subroutine r_inboard
Definition: mapping_mod.f90:1397
equilibrium_spider
subroutine equilibrium_spider(EQUILIBRIUM_IN, COREPROF_IN, EQUILIBRIUM_OUT, code_parameters)
Definition: equilibrium_spider.f90:3
assign_spider_parameters
subroutine assign_spider_parameters(codeparameters, return_status)
Definition: equilibrium_spider.f90:1162
helena21_mod::fdia
real(r8) function fdia(psi_n)
Definition: helena21_mod.f90:139
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
integral_value
subroutine integral_value(N, X, Y, INTY)
Definition: prepare_input_cpos.F90:1897
mod_spider::int_spider
subroutine int_spider(
Definition: mod_spider.f90:9