ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
emeq.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
10 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
11 MODULE emeq
12 
13 CONTAINS
14 
15  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
16  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
17 
18  SUBROUTINE emeq_e3m (EQUILIBRIUM_IN, EQUILIBRIUM_OUT, code_parameters)
19 
20  !-------------------------------------------------------!
21  ! This routine provides interface between ETS !
22  ! and three moment equilibrium solver. !
23  !-------------------------------------------------------!
24  ! Source: --- !
25  ! Developers: D.Kalupin, G.Pereverzev !
26  ! Kontacts: D.Kalupin@fz-juelich.de !
27  ! !
28  ! Comments: created for V&V between ETS and !
29  ! ASTRA !
30  ! !
31  !-------------------------------------------------------!
32 
33  USE itm_types
34  USE itm_constants
35  USE euitm_schemas
36  USE copy_structures
37  USE euitm_xml_parser
38  USE deallocate_structures
39 
40 
41  IMPLICIT NONE
42 
43  INTEGER :: ifail
44 
45 ! +++ CPO derived types:
46  TYPE (type_equilibrium), POINTER :: equilibrium_in(:) !input CPO with geometry quantities from previous time
47  TYPE (type_equilibrium), POINTER :: equilibrium_out(:) !output CPO with geometry quantities from previous iteration
48  TYPE (type_param) :: code_parameters
49  INTEGER :: return_status
50 
51 
52 ! +++ Dimensions:
53  INTEGER :: npsi !number of radial points (input, determined from EQUILIBRIUM_IN CPO)
54  INTEGER :: ndim1, ndim2 !number of points for 2-D parameters
55  INTEGER :: idim1, idim2 !index of points for 2-D parameters
56  INTEGER :: max_npoints !index of boundary points
57  INTEGER :: ipsi, itheta
58 
59 ! +++ Internal parameters(0-D parameters):
60  REAL (R8) :: r0 !magnetic axis (on first iteration equal to geom.axis)
61  REAL (R8) :: b0 !magnetic field on axis
62  REAL (R8) :: pc !total plasma current
63  REAL (R8) :: delta !Shafranov shift = 0 !!!not defined in EQUILIBRIUM CPO
64  REAL (R8) :: el !elongation
65  REAL (R8) :: tr !triangularity
66  REAL (R8) :: amin !minor radius
67  REAL (R8) :: theta
68  REAL (R8) :: rgeo
69  REAL (R8) :: zgeo
70  REAL (R8) :: bgeo
71  REAL (R8) :: rmag
72  REAL (R8) :: bmag
73 
74 
75 ! +++ 1-D parameters:
76  REAL (R8), ALLOCATABLE :: rho(:)
77  REAL (R8), ALLOCATABLE :: jpar(:), intjpar(:)
78  REAL (R8), ALLOCATABLE :: qsf(:)
79  REAL (R8), ALLOCATABLE :: pr(:)
80  REAL (R8), ALLOCATABLE :: psi(:)
81 
82  REAL (R8), ALLOCATABLE :: b2b02(:)
83  REAL (R8), ALLOCATABLE :: g1(:), g2(:), g3(:)
84  REAL (R8), ALLOCATABLE :: g4(:), g5(:), g6(:)
85  REAL (R8), ALLOCATABLE :: g7(:)
86  REAL (R8), ALLOCATABLE :: fdia(:)
87  REAL (R8), ALLOCATABLE :: vol(:), vprime(:), vprime_psi(:)
88  REAL (R8), ALLOCATABLE :: triang(:)
89  REAL (R8), ALLOCATABLE :: area(:), areaprime(:)
90  REAL (R8), ALLOCATABLE :: elong(:), shafr_shift(:)
91  REAL (R8), ALLOCATABLE :: pprime(:), ffprime(:)
92  ! Work
93  REAL (R8), ALLOCATABLE :: amid(:)
94  REAL (R8), ALLOCATABLE :: fun(:)
95 
96 
97 ! +++ 1-D parameters for 2-D output:
98  REAL (R8), ALLOCATABLE :: amid_2d(:), rho_2d(:)
99  REAL (R8), ALLOCATABLE :: psi_2d(:), jpar_2d(:), pr_2d(:)
100  REAL (R8), ALLOCATABLE :: shafr_shift_2d(:)
101  REAL (R8), ALLOCATABLE :: elong_2d(:), triang_l2d(:), triang_u2d(:)
102 
103 
104  REAL (R8), SAVE :: aceqlb
105  REAL (R8), SAVE :: convergence
106  INTEGER, SAVE :: itermax
107  INTEGER, SAVE :: neq !number of equilibrium points (input parameter)
108 
109  LOGICAL, SAVE :: first=.true.
110 
111  REAL (R8) :: s_ip, s_bt, s_q, s_psi, s_jpar
112 
113 
114 
115 !-------------------------------------------------------
116 !-------------------------------------------------------
117 
118  IF(first) THEN
119  aceqlb = 1.0e-3_r8
120  convergence = 1.0e-4_r8
121  itermax = 1000
122 
123 
124  IF(ASSOCIATED(equilibrium_in)) THEN
125  IF(ASSOCIATED(equilibrium_in(1)%profiles_1d%psi)) THEN
126  neq=SIZE(equilibrium_in(1)%profiles_1d%psi)
127  WRITE(*,*) 'Initial NEQ set to ',neq,' based on input equilibrium size of PSI'
128  ELSEIF(ASSOCIATED(equilibrium_in(1)%profiles_1d%rho_vol)) THEN
129  neq=SIZE(equilibrium_in(1)%profiles_1d%rho_vol)
130  WRITE(*,*) 'Initial NEQ set to ',neq,' based on input equilibrium size of RHO_VOL'
131  ELSEIF(ASSOCIATED(equilibrium_in(1)%profiles_1d%rho_tor)) THEN
132  neq=SIZE(equilibrium_in(1)%profiles_1d%rho_tor)
133  WRITE(*,*) 'Initial NEQ set to ',neq,' based on input equilibrium size of RHO_TOR'
134  ELSE
135  neq=100
136  WRITE(*,*) 'Initial NEQ set to ',neq
137  ENDIF
138  ELSE
139  neq=100
140  WRITE(*,*) 'Initial NEQ set to ',neq
141  ENDIF
142 
143  IF (.NOT. ASSOCIATED(code_parameters%parameters)) THEN
144  WRITE(6, *) 'ERROR: code parameters not associated!'
145  stop
146  ELSE
147  IF (.NOT.ASSOCIATED(equilibrium_in(1)%codeparam%parameters)) &
148  ALLOCATE(equilibrium_in(1)%codeparam%parameters(SIZE(code_parameters%parameters)))
149  call deallocate_cpo(equilibrium_in(1)%codeparam%parameters)
150  CALL copy_cpo(code_parameters%parameters, equilibrium_in(1)%codeparam%parameters)
151 
152  CALL assign_code_parameters(code_parameters, return_status)
153 
154  WRITE(*,*) 'Final NEQ set to ',neq
155 
156  IF (return_status /= 0) THEN
157  WRITE(*,*) 'ERROR: Could not assign EMEQ parameters.'
158  END IF
159 
160  END IF
161 
162  first=.false.
163 
164  END IF
165 
166 
167 
168 
169 ! +++ Set dimensions and allocate parameters:
170  npsi = SIZE (equilibrium_in(1)%profiles_1d%rho_tor, dim=1)
171  ndim1 = SIZE (equilibrium_in(1)%coord_sys%position%R, dim=1)
172  ndim2 = SIZE (equilibrium_in(1)%coord_sys%position%R, dim=2)
173  max_npoints = SIZE (equilibrium_in(1)%eqgeometry%boundary(1)%r, dim=1)
174 
175  s_ip = sign(1.0_r8, equilibrium_in(1)%global_param%i_plasma)
176  s_bt = sign(1.0_r8, equilibrium_in(1)%global_param%toroid_field%b0)
177  s_q = - s_ip * s_bt
178  s_psi = - s_ip
179  s_jpar = s_ip
180 
181  write(*,*) 's_ip, s_bt, s_q, s_psi, s_jpar = ', s_ip, s_bt, s_q, s_psi, s_jpar
182 
183  ALLOCATE (rho(npsi))
184  ALLOCATE (jpar(npsi))
185  ALLOCATE (intjpar(npsi))
186  ALLOCATE (qsf(npsi))
187  ALLOCATE (pr(npsi))
188  ALLOCATE (psi(npsi))
189  ALLOCATE (b2b02(npsi))
190  ALLOCATE (g1(npsi))
191  ALLOCATE (g2(npsi))
192  ALLOCATE (g3(npsi))
193  ALLOCATE (g4(npsi))
194  ALLOCATE (g5(npsi))
195  ALLOCATE (g6(npsi))
196  ALLOCATE (g7(npsi))
197  ALLOCATE (fdia(npsi))
198  ALLOCATE (vol(npsi))
199  ALLOCATE (vprime(npsi))
200  ALLOCATE (vprime_psi(npsi))
201  ALLOCATE (area(npsi))
202  ALLOCATE (areaprime(npsi))
203  ALLOCATE (amid(npsi))
204  ALLOCATE (shafr_shift(npsi))
205  ALLOCATE (elong(npsi))
206  ALLOCATE (triang(npsi))
207  ALLOCATE (pprime(npsi))
208  ALLOCATE (ffprime(npsi))
209  ALLOCATE (fun(npsi))
210 
211  ALLOCATE (amid_2d(ndim1))
212  ALLOCATE (rho_2d(ndim1))
213  ALLOCATE (elong_2d(ndim1))
214  ALLOCATE (triang_l2d(ndim1))
215  ALLOCATE (triang_u2d(ndim1))
216  ALLOCATE (shafr_shift_2d(ndim1))
217  ALLOCATE (psi_2d(ndim1))
218  ALLOCATE (jpar_2d(ndim1))
219  ALLOCATE (pr_2d(ndim1))
220 
221 
222  ! +++ Set parameters:
223  r0 = equilibrium_in(1)%global_param%toroid_field%r0
224  b0 = equilibrium_in(1)%global_param%toroid_field%b0 * s_bt
225 
226  rgeo = equilibrium_in(1)%eqgeometry%geom_axis%r
227  zgeo = equilibrium_in(1)%eqgeometry%geom_axis%z
228  bgeo = r0 * b0 / rgeo
229 
230  rmag = equilibrium_in(1)%global_param%mag_axis%position%r
231  bmag = equilibrium_in(1)%global_param%mag_axis%bphi
232 
233  amin = equilibrium_in(1)%eqgeometry%a_minor
234 
235  delta = r0 - rgeo
236  el = equilibrium_in(1)%eqgeometry%elongation
237  tr = (equilibrium_in(1)%eqgeometry%tria_upper + equilibrium_in(1)%eqgeometry%tria_lower) / 2.0_r8
238  amin = equilibrium_in(1)%eqgeometry%a_minor
239 
240  rho = equilibrium_in(1)%profiles_1d%rho_tor
241  jpar = equilibrium_in(1)%profiles_1d%jparallel * s_ip
242  qsf = equilibrium_in(1)%profiles_1d%q * s_q
243  pr = equilibrium_in(1)%profiles_1d%pressure
244  b2b02 = equilibrium_in(1)%profiles_1d%gm5 / b0**2
245  fdia = equilibrium_in(1)%profiles_1d%F_dia * s_bt / r0 / b0
246  psi = equilibrium_in(1)%profiles_1d%psi * s_psi
247 
248  pc = equilibrium_in(1)%global_param%i_plasma / 1.0e6_r8 * s_ip
249 
250  DO ipsi=1,npsi
251  IF (jpar(ipsi).LT.0.0_r8) jpar(ipsi) = 0.0_r8
252  END DO
253 
254  IF (ASSOCIATED(equilibrium_in(1)%profiles_1d%r_inboard).AND. &
255  ASSOCIATED(equilibrium_in(1)%profiles_1d%r_outboard)) THEN
256  amid = 0.5_r8 * (equilibrium_in(1)%profiles_1d%r_outboard - equilibrium_in(1)%profiles_1d%r_inboard)
257 
258  ELSE
259  amid = rho / sqrt(el)
260  ENDIF
261 
262  DO ipsi = 2, npsi
263  IF (amid(ipsi).LT.amid(ipsi-1)) THEN
264  WRITE(*,*) 'ERROR : AMID is non monotonic.'
265  amid = rho / sqrt(el)
266  END IF
267  END DO
268 
269  IF (amid(npsi).NE.amin) THEN
270  WRITE(*,*) 'ERROR : AMID(NPSI) is diffrent from AMIN.',amid(npsi),amin
271  amid = amid * amin / amid(npsi)
272  END IF
273 
274 
275  DO ipsi = 2, npsi
276  IF (amid(ipsi).LT.amid(ipsi-1)) THEN
277  WRITE(*,*) 'ERROR : AMID is non monotonic.'
278  amid = rho / sqrt(el)
279  END IF
280  END DO
281 
282  IF (amid(npsi).NE.amin) THEN
283  WRITE(*,*) 'ERROR : AMID(NPSI) is diffrent from AMIN.',amid(npsi),amin
284  amid = amid * amin / amid(npsi)
285  END IF
286 
287 
288 
289 
290 
291 
292 
293 ! +++ Call three moment equilibrium solver from ASTRA (EMEQ)
294 ! (input/output is on the RHO grid)
295  CALL emeq_interface( &
296  !
297  ! INPUT FROM THE TRANSPORT CODE TO EQUILIBRIUM:
298  ! values:
299  rgeo, & ! - Geometrical major radius [m]
300  bgeo, & ! - Toroidal magnetic field [T]
301  pc, & ! - Plasma current [A]
302  delta, & ! - Shift of plasma column [m]
303  el, & ! - elongation, boundary value
304  tr, & ! - triangularity, boundary value
305  amin, & ! - minor radius [m]
306  npsi, & ! - Number of radial points
307  neq, & ! - Number of equilibrium points
308  aceqlb, & !
309  convergence,& !
310  itermax, & !
311  ! arrays(NPSI):
312  rho, & ! - Rho [m]
313  amid, & ! - Normalized minor radius, [m]
314  jpar, & ! - Parallel current density [A/m^2]
315  qsf, & ! - safety factor
316  pr, & ! - pressure profile [Pa]
317  !
318  ! OUTPUT FROM THE EQUILIBRIUM TO THE TRANSPORT CODE:
319  g1, & ! - <1/R^2> [m^-2]
320  g2, & ! - <grad_rho^2/R^2> [m^-2]
321  g3, & ! - <grad_rho^2> [-]
322  g4, & ! - <1/B^2> [1/T^2]
323  g5, & ! - <B^2> [T^2]
324  g6, & ! - <grad_rho^2/B^2> [m^2/T^2]
325  g7, & ! - <grad_rho> [-]
326  b2b02, & ! - <B**2/B0**2> [-]
327  vol, & ! - volume (V) [m^3]
328  vprime, & ! - volume derivative (V') [m^2]
329  area, & ! - cross-sectional area [m^2]
330  areaprime, & ! - c.-s. area derivative [m]
331  fdia, & ! - diamagnetic function [-]
332  shafr_shift,& !
333  elong, & !
334  triang) !
335 
336 
337 
338 
339 ! +++ Save output to EQUILIBRIUM CPO:
340 ! +++ Allocate output CPO:
341 ! IF(ASSOCIATED(EQUILIBRIUM_OUT)) CALL DEALLOCATE_CPO(EQUILIBRIUM_OUT)
342  WRITE(*,*) 'Allocating EQUILIBRIUM_OUT'
343  ALLOCATE(equilibrium_out(1))
344 
345 
346 
347 
348 ! +++ Fill new EQUILIBRIUM CPO with quantities calculated by EMEQ solver:
349 ! 1-D profiles:
350  CALL l3deriv(vol, psi, neq, vprime_psi, psi, neq)
351  IF(vol(1).LT.0) THEN
352  WRITE(*,*) 'VOL(1) = ',vol(1)
353  vol(1)=0.0_r8
354  ENDIF
355 
356 
357 ! Constraint:
358  call deallocate_cpo(equilibrium_out(1)%eqconstraint)
359  CALL copy_cpo(equilibrium_in(1)%eqconstraint, equilibrium_out(1)%eqconstraint)
360 
361 
362 ! Eq. Geometry:
363  call deallocate_cpo(equilibrium_out(1)%eqgeometry)
364  CALL copy_cpo(equilibrium_in(1)%eqgeometry, equilibrium_out(1)%eqgeometry)
365 
366 
367  equilibrium_out(1)%eqgeometry%tria_upper = triang(npsi)
368  equilibrium_out(1)%eqgeometry%tria_lower = triang(npsi)
369 
370 
371 
372 ! GLOBAL_PARAM
373  call deallocate_cpo(equilibrium_out(1)%global_param)
374  CALL copy_cpo(equilibrium_in(1)%global_param, equilibrium_out(1)%global_param)
375 
376 
377  equilibrium_out(1)%global_param%area = area(npsi)
378  equilibrium_out(1)%global_param%volume = vol(npsi)
379  equilibrium_out(1)%global_param%i_plasma = pc * 1.0e6_r8 * s_ip
380  equilibrium_out(1)%global_param%psi_ax = psi(1) * s_psi
381  equilibrium_out(1)%global_param%psi_bound = psi(npsi) * s_psi
382  equilibrium_out(1)%global_param%mag_axis%position%R = rgeo + shafr_shift(1)
383  equilibrium_out(1)%global_param%mag_axis%position%Z = 0.0_r8
384  equilibrium_out(1)%global_param%mag_axis%bphi = fdia(1) * r0 * b0 * s_bt / (rgeo + shafr_shift(1))
385  equilibrium_out(1)%global_param%mag_axis%q = qsf(1) * s_q
386  equilibrium_out(1)%global_param%q_min = minval(qsf) * s_q
387 
388 
389 
390 ! PROFILES_1D
391 ! Allocate 1-D profiles:
392  ALLOCATE (equilibrium_out(1)%profiles_1d%psi(npsi))
393  ALLOCATE (equilibrium_out(1)%profiles_1d%phi(npsi))
394  ALLOCATE (equilibrium_out(1)%profiles_1d%pressure(npsi))
395  ALLOCATE (equilibrium_out(1)%profiles_1d%F_dia(npsi))
396  ALLOCATE (equilibrium_out(1)%profiles_1d%pprime(npsi))
397  ALLOCATE (equilibrium_out(1)%profiles_1d%ffprime(npsi))
398  ALLOCATE (equilibrium_out(1)%profiles_1d%jphi(npsi))
399  ALLOCATE (equilibrium_out(1)%profiles_1d%jparallel(npsi))
400  ALLOCATE (equilibrium_out(1)%profiles_1d%q(npsi))
401  ALLOCATE (equilibrium_out(1)%profiles_1d%r_inboard(npsi))
402  ALLOCATE (equilibrium_out(1)%profiles_1d%r_outboard(npsi))
403  ALLOCATE (equilibrium_out(1)%profiles_1d%rho_tor(npsi))
404  ALLOCATE (equilibrium_out(1)%profiles_1d%rho_vol(npsi))
405  ALLOCATE (equilibrium_out(1)%profiles_1d%beta_pol(npsi))
406  ALLOCATE (equilibrium_out(1)%profiles_1d%li(npsi))
407  ALLOCATE (equilibrium_out(1)%profiles_1d%elongation(npsi))
408  ALLOCATE (equilibrium_out(1)%profiles_1d%tria_upper(npsi))
409  ALLOCATE (equilibrium_out(1)%profiles_1d%tria_lower(npsi))
410  ALLOCATE (equilibrium_out(1)%profiles_1d%volume(npsi))
411  ALLOCATE (equilibrium_out(1)%profiles_1d%vprime(npsi))
412  ALLOCATE (equilibrium_out(1)%profiles_1d%area(npsi))
413  ALLOCATE (equilibrium_out(1)%profiles_1d%aprime(npsi))
414  ALLOCATE (equilibrium_out(1)%profiles_1d%surface(npsi))
415  ALLOCATE (equilibrium_out(1)%profiles_1d%ftrap(npsi))
416  ALLOCATE (equilibrium_out(1)%profiles_1d%dpsidrho_tor(npsi))
417  ALLOCATE (equilibrium_out(1)%profiles_1d%dvdrho(npsi))
418  ALLOCATE (equilibrium_out(1)%profiles_1d%b_av(npsi))
419 
420  ALLOCATE (equilibrium_out(1)%profiles_1d%gm1(npsi))
421  ALLOCATE (equilibrium_out(1)%profiles_1d%gm2(npsi))
422  ALLOCATE (equilibrium_out(1)%profiles_1d%gm3(npsi))
423  ALLOCATE (equilibrium_out(1)%profiles_1d%gm4(npsi))
424  ALLOCATE (equilibrium_out(1)%profiles_1d%gm5(npsi))
425  ALLOCATE (equilibrium_out(1)%profiles_1d%gm6(npsi))
426  ALLOCATE (equilibrium_out(1)%profiles_1d%gm7(npsi))
427  ALLOCATE (equilibrium_out(1)%profiles_1d%gm8(npsi))
428  ALLOCATE (equilibrium_out(1)%profiles_1d%gm9(npsi))
429 
430 
431  equilibrium_out(1)%profiles_1d%psi = psi * s_psi
432  equilibrium_out(1)%profiles_1d%rho_tor = rho
433  equilibrium_out(1)%profiles_1d%phi = rho**2 * itm_pi * b0 * s_bt
434  equilibrium_out(1)%profiles_1d%q = qsf * s_q
435  equilibrium_out(1)%profiles_1d%pressure = pr
436  equilibrium_out(1)%profiles_1d%jparallel = jpar * s_ip
437  equilibrium_out(1)%profiles_1d%volume = vol
438  equilibrium_out(1)%profiles_1d%vprime = vprime_psi * s_psi
439  equilibrium_out(1)%profiles_1d%area = area
440  equilibrium_out(1)%profiles_1d%aprime = areaprime
441  equilibrium_out(1)%profiles_1d%F_dia = fdia * r0 * b0 * s_bt
442  equilibrium_out(1)%profiles_1d%elongation = elong
443  equilibrium_out(1)%profiles_1d%tria_upper = triang
444  equilibrium_out(1)%profiles_1d%tria_lower = triang
445  equilibrium_out(1)%profiles_1d%r_inboard = rgeo + shafr_shift - amid
446  equilibrium_out(1)%profiles_1d%r_outboard = rgeo + shafr_shift + amid
447  equilibrium_out(1)%profiles_1d%rho_vol = sqrt(vol/vol(npsi))
448 
449  equilibrium_out(1)%profiles_1d%gm1 = g1
450  equilibrium_out(1)%profiles_1d%gm2 = g2
451  equilibrium_out(1)%profiles_1d%gm3 = g3
452  equilibrium_out(1)%profiles_1d%gm4 = g4
453  equilibrium_out(1)%profiles_1d%gm5 = g5
454  equilibrium_out(1)%profiles_1d%gm6 = g6
455  equilibrium_out(1)%profiles_1d%gm7 = g7
456  equilibrium_out(1)%profiles_1d%gm8 = 1._r8/g1**0.5
457  equilibrium_out(1)%profiles_1d%gm9 = g1**0.5
458 
459  CALL l3deriv(pr, rho, npsi, fun, rho, npsi)
460  pprime = (rgeo + shafr_shift) * qsf * s_q / b0 * s_bt / rho * fun !!!!!- was removed
461  IF (rho(1).EQ.0.0_r8) pprime(1) = pprime(2)
462  b2b02 = g5 / b0**2
463  ffprime = fdia * s_bt / b2b02 * (jpar * s_jpar * (rgeo + shafr_shift) / r0 - pprime)
464 
465  equilibrium_out(1)%profiles_1d%pprime = pprime
466  equilibrium_out(1)%profiles_1d%ffprime = ffprime
467  equilibrium_out(1)%profiles_1d%dvdrho = vprime
468  equilibrium_out(1)%profiles_1d%b_av = g5**0.5
469 
470 
471 
472 ! Stabilization of current oscillations should be removed later
473  CALL l3interp(equilibrium_in(1)%profiles_1d%gm2, &
474  equilibrium_in(1)%profiles_1d%rho_tor, &
475  SIZE (equilibrium_in(1)%profiles_1d%rho_tor), &
476  equilibrium_out(1)%profiles_1d%gm2, &
477  equilibrium_out(1)%profiles_1d%rho_tor, &
478  SIZE (equilibrium_out(1)%profiles_1d%rho_tor))
479  equilibrium_out(1)%profiles_1d%gm2 = equilibrium_out(1)%profiles_1d%gm2 * 0.9_r8 + g2 * 0.1_r8
480 
481 ! COORD_SYS
482  ALLOCATE(equilibrium_out(1)%coord_sys%grid_type(4))
483  ALLOCATE(equilibrium_out(1)%coord_sys%grid%dim1(ndim1))
484  ALLOCATE(equilibrium_out(1)%coord_sys%grid%dim2(ndim2))
485  ALLOCATE(equilibrium_out(1)%coord_sys%position%R(ndim1,ndim2))
486  ALLOCATE(equilibrium_out(1)%coord_sys%position%Z(ndim1,ndim2))
487 
488  rho_loop2: DO idim1=1,ndim1
489  amid_2d(idim1) = amin/(ndim1-1)*(idim1-1)
490  END DO rho_loop2
491  CALL l3interp(rho, amid, npsi, rho_2d, amid_2d, ndim1)
492  CALL l3interp(elong, amid, npsi, elong_2d, amid_2d, ndim1)
493  CALL l3interp(triang, amid, npsi, triang_l2d, amid_2d, ndim1)
494  CALL l3interp(triang, amid, npsi, triang_u2d, amid_2d, ndim1)
495  CALL l3interp(shafr_shift, amid, npsi, shafr_shift_2d, amid_2d, ndim1)
496  CALL l3interp(psi, amid, npsi, psi_2d, amid_2d, ndim1)
497  fun = jpar * 1.0e6_r8
498  CALL l3interp(fun, amid, npsi, jpar_2d, amid_2d, ndim1)
499  fun = pr * 1.0e6_r8
500  CALL l3interp(fun, amid, npsi, pr_2d, amid_2d, ndim1)
501 
502 
503  DO itheta=1, ndim2
504  theta = REAL(itheta-1,r8)/REAL(ndim2,r8)*2.0_r8*itm_pi
505  DO ipsi=1, ndim1
506  equilibrium_out(1)%coord_sys%grid_type(1) = '2'
507  equilibrium_out(1)%coord_sys%grid_type(2) = 'inverse'
508  equilibrium_out(1)%coord_sys%grid_type(3) = '3'
509  equilibrium_out(1)%coord_sys%grid_type(4) = 'polar'
510  equilibrium_out(1)%coord_sys%grid%dim1(ipsi) = psi_2d(ipsi) * s_psi
511  equilibrium_out(1)%coord_sys%grid%dim2(itheta) = theta
512  equilibrium_out(1)%coord_sys%position%R(ipsi, itheta) = r0 + shafr_shift_2d(ipsi) + amid_2d(ipsi) &
513  * (cos(theta)-(triang_u2d(ipsi)+triang_l2d(ipsi))/2._r8*(sin(theta))**2)
514  equilibrium_out(1)%coord_sys%position%Z(ipsi, itheta) = zgeo + amid_2d(ipsi) * elong_2d(ipsi) * sin(theta)
515 
516  END DO
517  END DO
518 
519 
520 
521 ! EQ_GEOMETRY
522  IF(ASSOCIATED(equilibrium_out(1)%eqgeometry%boundary)) &
523  CALL deallocate_cpo(equilibrium_out(1)%eqgeometry%boundary)
524  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1))
525  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1)%r(max_npoints))
526  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1)%z(max_npoints))
527 
528  DO itheta=1, max_npoints
529  theta = REAL(itheta-1,r8)/REAL(max_npoints-1,r8)*2.0_r8*itm_pi
530  equilibrium_out(1)%eqgeometry%boundary(1)%r(itheta) = r0 + shafr_shift_2d(npsi) + amin * (cos(theta) - tr*(sin(theta))**2)
531  equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta) = zgeo + amin * el * sin(theta)
532  IF (abs(equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta)).LE.1.e-15_r8) &
533  equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta) = 0.0e0_r8
534  END DO
535  equilibrium_out(1)%eqgeometry%boundarytype = 1
536 
537 
538  equilibrium_out(1)%eqgeometry%geom_axis%r = rgeo
539  equilibrium_out(1)%eqgeometry%geom_axis%z = zgeo
540  equilibrium_out(1)%eqgeometry%a_minor = amin
541  equilibrium_out(1)%eqgeometry%elong_upper = el
542  equilibrium_out(1)%eqgeometry%elong_lower = el
543 
544 
545 
546 
547 
548 ! Allocate 2-D profiles:
549  ALLOCATE (equilibrium_out(1)%profiles_2d(1))
550  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%grid_type(4))
551  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%grid%dim1(ndim1))
552  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%grid%dim2(ndim2))
553  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%R(ndim1, ndim2))
554  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%Z(ndim1, ndim2))
555  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%psi(ndim1, ndim2))
556  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%theta(ndim1, ndim2))
557  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%jphi(ndim1, ndim2))
558  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%jpar(ndim1, ndim2))
559  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%br(ndim1, ndim2))
560  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%bz(ndim1, ndim2))
561  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%bphi(ndim1, ndim2))
562  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%vphi(ndim1, ndim2))
563  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%vtheta(ndim1, ndim2))
564  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%rho_mass(ndim1, ndim2))
565  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%pressure(ndim1, ndim2))
566  ALLOCATE (equilibrium_out(1)%profiles_2d(1)%temperature(ndim1, ndim2))
567 
568 
569  DO itheta=1, ndim2
570  theta = REAL(itheta-1,r8)/REAL(ndim2-1,r8)*2.0_r8*itm_pi
571 
572  DO ipsi=1, ndim1
573  equilibrium_out(1)%profiles_2d(1)%grid_type(1) = '2'
574  equilibrium_out(1)%profiles_2d(1)%grid_type(2) = 'inverse'
575  equilibrium_out(1)%profiles_2d(1)%grid_type(3) = '3'
576  equilibrium_out(1)%profiles_2d(1)%grid_type(4) = 'polar'
577  equilibrium_out(1)%profiles_2d(1)%grid%dim1(ipsi) = psi_2d(ipsi) * s_psi
578  equilibrium_out(1)%profiles_2d(1)%grid%dim2(itheta) = theta
579  IF (theta.LE.itm_pi) tr = triang_u2d(ipsi)
580  IF (theta.GT.itm_pi) tr = triang_l2d(ipsi)
581 
582  equilibrium_out(1)%profiles_2d(1)%R(ipsi, itheta) = r0 + shafr_shift_2d(ipsi) + amid_2d(ipsi) &
583  * (cos(theta)-tr*(sin(theta))**2)
584  equilibrium_out(1)%profiles_2d(1)%Z(ipsi, itheta) = zgeo + amid_2d(ipsi) * elong_2d(ipsi) * sin(theta)
585  equilibrium_out(1)%profiles_2d(1)%theta(ipsi, itheta) = theta
586 
587  equilibrium_out(1)%profiles_2d(1)%psi(ipsi, itheta) = psi_2d(ipsi) * s_psi
588  equilibrium_out(1)%profiles_2d(1)%jpar(ipsi, itheta) = jpar_2d(ipsi) * s_jpar
589  equilibrium_out(1)%profiles_2d(1)%pressure(ipsi, itheta) = pr_2d(ipsi)
590  END DO
591  END DO
592 
593 
594 
595 
596 
597 ! CODE_PARAM
598  ALLOCATE (equilibrium_out(1)%codeparam%codename(1))
599  ALLOCATE (equilibrium_out(1)%codeparam%codeversion(1))
600  ALLOCATE (equilibrium_out(1)%codeparam%parameters(SIZE(code_parameters%parameters)))
601  ALLOCATE (equilibrium_out(1)%codeparam%output_diag(1))
602 
603 
604  equilibrium_out(1)%codeparam%codename(1) = 'emeq'
605  equilibrium_out(1)%codeparam%codeversion(1) = '01.03.2013'
606  equilibrium_out(1)%codeparam%output_diag(1) = 'my_output_diag'
607  equilibrium_out(1)%codeparam%output_flag = 0
608  equilibrium_out(1)%codeparam%parameters = code_parameters%parameters
609 
610 
611 
612 ! DATA_INFO
613  ALLOCATE (equilibrium_out(1)%datainfo%dataprovider(1))
614  equilibrium_out(1)%datainfo%dataprovider(1) = 'Denis Kalupin'
615  equilibrium_out(1)%datainfo%cocos = 13 !!!!12
616 
617 
618 
619 
620  ! +++ Deallocate parameters:
621  DEALLOCATE (rho)
622  DEALLOCATE (jpar)
623  DEALLOCATE (intjpar)
624  DEALLOCATE (qsf)
625  DEALLOCATE (pr)
626  DEALLOCATE (psi)
627  DEALLOCATE (b2b02)
628  DEALLOCATE (g1)
629  DEALLOCATE (g2)
630  DEALLOCATE (g3)
631  DEALLOCATE (g4)
632  DEALLOCATE (g5)
633  DEALLOCATE (g6)
634  DEALLOCATE (g7)
635  DEALLOCATE (fdia)
636  DEALLOCATE (vol)
637  DEALLOCATE (vprime)
638  DEALLOCATE (vprime_psi)
639  DEALLOCATE (area)
640  DEALLOCATE (areaprime)
641  DEALLOCATE (amid)
642  DEALLOCATE (shafr_shift)
643  DEALLOCATE (elong)
644  DEALLOCATE (triang)
645  DEALLOCATE (pprime)
646  DEALLOCATE (ffprime)
647  DEALLOCATE (fun)
648 
649  DEALLOCATE (amid_2d)
650  DEALLOCATE (rho_2d)
651  DEALLOCATE (elong_2d)
652  DEALLOCATE (triang_l2d)
653  DEALLOCATE (triang_u2d)
654  DEALLOCATE (shafr_shift_2d)
655  DEALLOCATE (psi_2d)
656  DEALLOCATE (jpar_2d)
657  DEALLOCATE (pr_2d)
658 
659 
660 
661  RETURN
662 
663 
664 
665 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
666 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
667 
668 
669  CONTAINS
670 
671  SUBROUTINE assign_code_parameters(codeparameters, return_status)
672 
673  !-----------------------------------------------------------------------
674  ! calls the XML parser for the code parameters and assign the
675  ! resulting values to the corresponding variables
676  !TODO: check an alternative and more elegant solution in Perl
677  !-----------------------------------------------------------------------
678 
679  USE mod_f90_kind
680 
681  IMPLICIT NONE
682 
683  TYPE (type_param), INTENT(in) :: codeparameters
684  INTEGER(ikind), INTENT(out) :: return_status
685 
686  TYPE(tree) :: parameter_list
687  TYPE(element), POINTER :: temp_pointer
688  INTEGER(ikind) :: i, nparm, n_values
689  CHARACTER(len = 132) :: cname
690 
691  return_status = 0 ! no error
692 
693  !-- parse xml-string codeparameters
694 
695  WRITE(*,*) 'Calling euitm_xml_parse'
696  CALL euitm_xml_parse(codeparameters, nparm, parameter_list)
697  WRITE(*,*) 'Called euitm_xml_parse'
698 
699  !-- assign variables
700 
701  temp_pointer => parameter_list%first
702 
703  outer: DO
704  cname = char2str(temp_pointer%cname) ! necessary for AIX
705  SELECT CASE (cname)
706  CASE ("parameters")
707  temp_pointer => temp_pointer%child
708  cycle
709 
710 !-- dims parameters
711  CASE ("dims")
712  temp_pointer => temp_pointer%child
713  cycle
714  CASE ("neq")
715  IF (ALLOCATED(temp_pointer%cvalue)) &
716  CALL char2num(temp_pointer%cvalue, neq)
717 
718 !-- solver parameters
719  CASE ("solver")
720  temp_pointer => temp_pointer%child
721  cycle
722  CASE ("convergence")
723  IF (ALLOCATED(temp_pointer%cvalue)) &
724  CALL char2num(temp_pointer%cvalue, convergence)
725  CASE ("aceqlb")
726  IF (ALLOCATED(temp_pointer%cvalue)) &
727  CALL char2num(temp_pointer%cvalue, aceqlb)
728  CASE ("itermax")
729  IF (ALLOCATED(temp_pointer%cvalue)) &
730  CALL char2num(temp_pointer%cvalue, itermax)
731 
732  CASE default
733  WRITE(*, *) 'ERROR: invalid parameter', cname
734  return_status = 1
735  EXIT
736  END SELECT
737  DO
738  IF (ASSOCIATED(temp_pointer%sibling)) THEN
739  temp_pointer => temp_pointer%sibling
740  EXIT
741  END IF
742  IF (ASSOCIATED(temp_pointer%parent, parameter_list%first )) &
743  EXIT outer
744  IF (ASSOCIATED(temp_pointer%parent)) THEN
745  temp_pointer => temp_pointer%parent
746  ELSE
747  WRITE(*, *) 'ERROR: broken list.'
748  RETURN
749  END IF
750  END DO
751  END DO outer
752 
753  !-- destroy tree
754  CALL destroy_xml_tree(parameter_list)
755 
756  RETURN
757 
758  END SUBROUTINE assign_code_parameters
759 
760  SUBROUTINE read_codeparam(in_xml, filename, codeparam)
761 
762  USE euitm_schemas
763  USE ets_version
764  IMPLICIT NONE
765 
766  INTEGER n_lines, in_xml, ios, i
767  CHARACTER (len=*) :: filename
768  TYPE (type_codeparam) :: codeparam
769  CHARACTER(len = 132) :: xml_line
770 
771  OPEN (unit = in_xml, file = filename, status = 'old', &
772  action = 'read', iostat = ios)
773 
774  IF (ios /= 0) THEN
775  WRITE(*,*) 'Could not open ',trim(filename)
776  stop ' ERROR: XML file does not exist '
777  END IF
778 
779  n_lines = 0
780 
781  DO
782  READ (in_xml, '(a)', iostat = ios) xml_line
783  IF (ios == 0) THEN
784  n_lines = n_lines + 1
785  ELSE
786  EXIT
787  END IF
788  END DO
789 
790  rewind in_xml
791 
792  ALLOCATE(codeparam%codename(1))
793  codeparam%codename(1)='EMEQ'
794  ALLOCATE(codeparam%codeversion(1))
795  codeparam%codeversion(1)=version
796  WRITE(*,*) 'Code = ',trim(codeparam%codename(1)),' version = ',trim(codeparam%codeversion(1))
797  ALLOCATE(codeparam%parameters(n_lines))
798  DO i = 1, n_lines
799  READ (in_xml, '(a)', iostat = ios) codeparam%parameters(i)
800  END DO
801 
802  CLOSE(in_xml)
803 
804  RETURN
805  END SUBROUTINE read_codeparam
806 
807  END SUBROUTINE emeq_e3m
808  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
809  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
810 
811 
812 
813 
814 
815  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
816  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
817  SUBROUTINE emeq_interface( &
818  ! +++ EQUILIBRIUM:
819  ! +++ Interface to three moment solver from ASTRA (G.Pereversev) [EMEQ]
820  !
821  ! INPUT FROM THE TRANSPORT CODE TO EQUILIBRIUM:
822  ! values:
823  rgeo, & ! - Geometrical major radius [m]
824  bgeo, & ! - Toroidal magnetic field [T]
825  pc, & ! - Plasma current [A]
826  delta, & ! - Shift of plasma column [m]
827  eln, & ! - elongation, boundary value
828  trn, & ! - triangularity, boundary value
829  amin, & ! - minor radius [m]
830  npsi, & ! - Number of radial points
831  neq, & ! - Number of equilibrium points
832  aceqlb, & !
833  convergence,& !
834  itermax, & !
835  ! arrays(NPSI):
836  rho, & ! - Rho [m]
837  amid, & ! - Normalized minor radius, [m]
838  jpar, & ! - Parallel current density [A/m^2]
839  qsf, & ! - safety factor
840  pr, & ! - pressure profile [Pa]
841  !
842  ! OUTPUT FROM THE EQUILIBRIUM TO THE TRANSPORT CODE:
843  g1, & ! - <1/R^2> [m^-2]
844  g2, & ! - <grad_rho^2/R^2> [m^-2]
845  g3, & ! - <grad_rho^2> [-]
846  g4, & ! - <1/B^2> [1/T^2]
847  g5, & ! - <B^2> [T^2]
848  g6, & ! - <grad_rho^2/B^2> [m^2/T^2]
849  g7, & ! - <grad_rho> [-]
850  b2b02, & ! - <B**2/B0**2> [-]
851  vol, & ! - volume (V) [m^3]
852  vprime, & ! - volume derivative (V') [m^2]
853  area, & ! - cross-sectional area [m^2]
854  areaprime, & ! - c.-s. area derivative [m]
855  fdia, & ! - diamagnetic function [-]
856  shafr_shift,& !
857  elong, & !
858  triang) !
859 
860 
861 
862  USE itm_types
863  USE itm_constants
864  USE emeq_equil
865 
866 
867 
868  IMPLICIT NONE
869 
870  INTEGER :: npsi, neq, i, iter
871 
872  ! +++ input (0-D parameters):
873  REAL (R8) :: rgeo, bgeo, pc
874  REAL (R8) :: delta, eln, trn
875  REAL (R8) :: aceqlb, convergence
876  INTEGER :: itermax
877 
878  ! +++ input (1-D arrays):
879  REAL (R8) :: rho(npsi), rho_new(npsi)
880  REAL (R8) :: amid(npsi), amid_eq(neq)
881  REAL (R8) :: jpar(npsi), qsf(npsi)
882  REAL (R8) :: pr(npsi)
883 
884  ! +++ output (1-D arrays):
885  REAL (R8) :: g1(npsi), g1_eq(neq)
886  REAL (R8) :: g2(npsi), g2_eq(neq)
887  REAL (R8) :: g3(npsi), g3_eq(neq)
888  REAL (R8) :: g4(npsi), g4_eq(neq)
889  REAL (R8) :: g5(npsi), g5_eq(neq)
890  REAL (R8) :: g6(npsi), g6_eq(neq)
891  REAL (R8) :: g7(npsi), g7_eq(neq)
892  REAL (R8) :: area(npsi), area_eq(neq)
893  REAL (R8) :: areaprime(npsi), areaprime_eq(neq)
894  REAL (R8) :: vol(npsi), vol_eq(neq)
895  REAL (R8) :: vprime(npsi), vprime_eq(neq)
896  REAL (R8) :: fdia(npsi), fdia_eq(neq)
897  REAL (R8) :: b2b02(npsi), b2b02_eq(neq)
898  REAL (R8) :: shafr_shift(npsi),gbd(neq)
899  REAL (R8) :: elong(npsi), gl(neq)
900  REAL (R8) :: triang(npsi), gsd(neq)
901 
902 
903  ! +++ internal quantities:
904  REAL (R8) :: ja(npsi), jaeq(neq)
905  REAL (R8) :: jb(npsi), jbeq(neq)
906 
907  REAL (R8) :: time
908  REAL (R8) :: r0, b0, amin
909  REAL (R8) :: dpr(npsi)
910  REAL (R8) :: nabrho(npsi), nabrho_eq(neq)
911  REAL (R8) :: vprime_n(npsi)
912  REAL (R8) :: nabrho_n(npsi)
913  REAL (R8) :: g2_n(npsi)
914  REAL (R8) :: qsf_n(npsi)
915  REAL (R8) :: fun(npsi), intfun(npsi)
916 
917  REAL (R8) :: gr(neq)
918  REAL (R8) :: gra(neq), sqgra(neq)
919  REAL (R8) :: grar(neq), avr2(neq)
920  REAL (R8) :: ai0(neq), dgrda(neq), avsqg(neq)
921  REAL (R8) :: grdaeq(neq), b2b0eq(neq)
922  REAL (R8) :: b0b2eq(neq), b0b0eq(neq), bmaxeq(neq)
923  REAL (R8) :: bmineq(neq), bmodeq(neq), fofbeq(neq)
924 
925  REAL (R8) :: conv, err_vprime, err_nabrho, err_qsf, err_g2
926 
927 
928 
929 
930 ! +++ ALLOCATE WORKING ARRAYS FOR EMEQ:
931  CALL alloc_emeq_equil(neq)
932 
933 
934 ! +++ INITIAL SETTINGS:
935  r0 = (rgeo+delta) ! Major radius
936  b0 = bgeo*rgeo/r0 ! Magnetic field at R0
937 
938  CALL l3deriv(pr, rho, npsi, dpr, rho, npsi)
939 
940 
941 ! +++ this might not be necessary
942  DO i=1,npsi
943  dpr(i) = min(0.0_r8,dpr(i)) ! negative derivative is not acceptable
944  END DO ! by the equilibrium solver
945 
946 
947 
948 ! +++ equidistant equilibrium grid:
949  DO i=1,neq
950  amid_eq(i) = amin /(neq-1)*(i-1)
951  END DO
952 
953 
954 
955  iter = 0
956 
957  vprime = 0.0_r8
958  nabrho = 0.0_r8
959  g2 = 0.0_r8
960 
961 ! +++ CALCULATION OF GEOMETRY COEFFICIENTS:
962 1 CONTINUE
963 
964  iter = iter +1
965 
966  IF(iter.GT.100) THEN
967  WRITE(*,*) '### ERROR! Equilibrium loop in EQUILIBRIUM_INTERFACE not converging'
968  stop 'Error - no convergence in EQUILIBRIUM_INTERFACE'
969  CALL dealloc_emeq_equil
970  RETURN
971 
972  ENDIF
973 
974 
975 
976 ! +++ PARAMETERS FOR ITERATIONS:
977  vprime_n = vprime !copy of equilibrium quantities for iterations
978  nabrho_n = nabrho
979  qsf_n = qsf
980  g2_n = g2
981 
982 
983 
984 ! +++ RHS of Grad-Shafranov equation:
985 2 DO i=2,npsi
986  jb(i) = -r0/b0*qsf(i)/rho(i)*dpr(i)*1.d-6 ![MA/m^2]
987 
988  ja(i) = jb(i)*(1.0_r8-fdia(i)**2/b2b02(i)*rgeo**2/r0**2) &
989  + jpar(i)/1.d6*fdia(i)/b2b02(i)*rgeo/r0 ![MA/m^2]
990  END DO
991 
992 
993  CALL axis_int(npsi, rho, jb)
994  CALL axis_int(npsi, rho, ja)
995 
996 
997 
998 ! +++ Interpolation of input coefficients on equilibrium grid:
999  CALL l3interp(ja, amid ,npsi, jaeq, amid_eq, neq)
1000  CALL l3interp(jb, amid ,npsi, jbeq, amid_eq, neq)
1001 
1002 
1003 
1004 ! +++ CALL TO THE EQUILIBRIUM SOLVER:
1005  CALL e3astr &
1006  !-------------------------------------------------------!
1007  ! The same equidistant grid with respect to the "radial"!
1008  ! variable "A" is assumed to be used both for input and !
1009  ! output. The radial coordinate x=a/a_edge, is !
1010  ! dimensionless, A(1)=0; A(NAEQ1)=1 !
1011  !-------------------------------------------------------!
1012  (jaeq,jbeq & ! j_zeta = BA*(R0/r) + BB*(r/R0-R0/r)[MA/m^2]
1013  ,r0 & ! RGEO+DELTA, Major radius of the edge FS [m]
1014  ,amin & ! minor radius [m]
1015  ,eln & ! elongation, edge value [-]
1016  ,trn*amin & ! triangularity, edge value * minor radius [m]
1017  ,neq & ! radial grid point No.! NEQL
1018  ,aceqlb &
1019  ,itermax & ! max number of iterations
1020  ,b0 & ! vacuum B @ R0 [T]
1021  ,pc & ! total plasma current [MA]
1022  ! Output:
1023  ,gr,gbd,gl,gsd & ! Rho(a), Delta(a), ELONG(a), delta(a)
1024  ,gra & ! <g^{11}>= <[nabla(a)]^2>
1025  ,sqgra & ! <\sqrt{|g^{11}|}>= <|nabla(a)|>
1026  ,grar & ! <g^{11}g^{33}> = <[nabla(a)/r]^2>
1027  ,avr2 & ! <g^{33}>= <1/r^2>=G33/R0**2
1028  ,ai0 & ! I-> IPOL*R0*BTOR
1029  ,dgrda & ! \prti{\rho}{a}
1030  ,avsqg & ! \prti{V}{a}/(4\pi^2)
1031  ,vol_eq & ! V(a)
1032  ,b2b0eq & ! <B**2/B0**2>
1033  ,b0b0eq & ! <B0**2/B**2>
1034  ,bmaxeq & ! BMAXT
1035  ,bmineq & ! BMINT
1036  ,bmodeq & ! <B/BTOR>
1037  ,fofbeq & ! <(BTOR/B)**2*(1.-SQRT(1-B/Bmax)*(1+.5B/Bmax))>
1038  ,grdaeq & ! <grad a>
1039  ,time)
1040 
1041 
1042 
1043 ! +++ OUTPUT PARAMETERS (on equilibrium grid):
1044 
1045 
1046  g1_eq = avr2
1047  g2_eq = grar * dgrda**2
1048  g3_eq = gra * dgrda**2
1049  g4_eq = b0b0eq / b0**2
1050  g5_eq = b2b0eq * b0**2
1051  g6_eq = gra * dgrda**2 / b0**2
1052  g7_eq = sqgra * dgrda
1053 
1054  fdia_eq = ai0 / (b0*r0) !diamagnetic function [-]
1055  b2b02_eq = b2b0eq
1056 
1057  nabrho_eq = gra * dgrda**2
1058  gsd = gsd / amid_eq !triangularity in ETS units
1059  CALL axis_int(neq, amid_eq, gsd)
1060 
1061  CALL l3deriv(vol_eq, gr, neq, vprime_eq, gr, neq) !V' on rho coordinate
1062 
1063  area_eq = vprime_eq * sqgra * dgrda !V'*<grad_rho> [m^2]
1064  CALL l3deriv(area_eq, gr, neq, areaprime_eq, gr, neq)!AREA' on rho coordinate
1065 
1066 
1067 
1068 
1069 ! +++ Interpolation of output quantities on RHO grid:
1070  DO i =1,npsi
1071  rho_new(i) = gr(neq)*(i-1)/(npsi-1)
1072  END DO
1073 
1074  fun = jpar
1075  CALL l3interp(fun, rho, npsi, jpar, rho_new, npsi)
1076  fun = dpr
1077  CALL l3interp(fun, rho, npsi, dpr, rho_new, npsi)
1078 
1079  rho = rho_new
1080 
1081 
1082  CALL l3interp(amid_eq, gr, neq, amid, rho, npsi)
1083  CALL l3interp(nabrho_eq, gr, neq, nabrho, rho, npsi)
1084  CALL l3interp(vol_eq, gr, neq, vol, rho, npsi)
1085  CALL l3interp(vprime_eq, gr, neq, vprime, rho, npsi)
1086  CALL l3interp(g1_eq, gr, neq, g1, rho, npsi)
1087  CALL l3interp(g2_eq, gr, neq, g2, rho, npsi)
1088  CALL l3interp(g3_eq, gr, neq, g3, rho, npsi)
1089  CALL l3interp(g4_eq, gr, neq, g4, rho, npsi)
1090  CALL l3interp(g5_eq, gr, neq, g5, rho, npsi)
1091  CALL l3interp(g6_eq, gr, neq, g6, rho, npsi)
1092  CALL l3interp(g7_eq, gr, neq, g7, rho, npsi)
1093  CALL l3interp(fdia_eq, gr, neq, fdia, rho, npsi)
1094  CALL l3interp(b2b02_eq, gr, neq, b2b02, rho, npsi)
1095  CALL l3interp(area_eq, gr, neq, area, rho, npsi)
1096  CALL l3interp(areaprime_eq, gr, neq, areaprime, rho, npsi)
1097  CALL l3interp(gbd, gr, neq, shafr_shift,rho, npsi)
1098  CALL l3interp(gl, gr, neq, elong , rho, npsi)
1099  CALL l3interp(gsd, gr, neq, triang, rho, npsi)
1100 
1101 
1102 ! +++ recalculation of the safety factor:
1103  intfun = 0.0_r8
1104  DO i=1,npsi
1105  fun(i) = jpar(i)*vprime(i)*itm_mu0/(fdia(i)*b0*r0)**2
1106  IF (rho(i).NE.0.0_r8.AND.i.GT.1) THEN
1107  intfun(i) = intfun(i-1) + (fun(i-1)+fun(i))*(rho(i)-rho(i-1))/2.0_r8
1108  qsf(i) = rho(i)*vprime(i)*g2(i)/(fdia(i)*b0*r0)/intfun(i)
1109  END IF
1110  END DO
1111  CALL axis_int(npsi, rho, qsf)
1112 
1113 
1114 
1115 ! +++ CONVERGENCY CHECK:
1116  conv=0.0_r8
1117 
1118  err_vprime = maxval(abs(1.0_r8-vprime_n(2:)/vprime(2:)))
1119  err_nabrho = maxval(abs(1.0_r8-nabrho_n/nabrho))
1120  err_qsf = maxval(abs(1.0_r8-qsf_n/qsf))
1121  err_g2 = maxval(abs(1.0_r8-g2_n/g2))
1122 
1123  WRITE(*,*) 'err_VPRIME ', err_vprime
1124  WRITE(*,*) 'err_NABRHO ', err_nabrho
1125  WRITE(*,*) ' err_QSF ', err_qsf
1126  WRITE(*,*) ' err_G2 ', err_g2
1127 
1128  conv = max( err_vprime, err_nabrho, err_qsf, err_g2)
1129 
1130  WRITE(*,*) ' CONV ', conv
1131 
1132  IF(conv.GT.convergence) goto 1
1133 
1134 
1135 
1136 ! +++ CHECK FOR NEGATIVE VOLUME:
1137  IF(vol(1).LT.0) THEN
1138  WRITE(*,*) 'VOL(1) = ',vol(1)
1139  vol(1) = 0.0_r8
1140  ENDIF
1141  IF(vprime(1).LT.0.0_r8) THEN
1142  WRITE(*,*) '### ', vprime(1),' changed to 0'
1143  vprime(1) = 0.0_r8
1144  ENDIF
1145 
1146 
1147  CALL dealloc_emeq_equil
1148 
1149 
1150  RETURN
1151 
1152 
1153 
1154  END SUBROUTINE emeq_interface
1155  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1156  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1157 
1158 
1159 
1160 
1161 
1162 
1163  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1164  ! This subroutine finds f(r_1=0) from f(r_2), f(r_3) and f(r_4)
1165  !
1166  ! author: M.Tokar (m.tokar@fz-juelich.de)
1167  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1168  SUBROUTINE axis_int(n, r, f)
1169  !-------------------------------------------------------!
1170  ! !
1171  ! This subroutine finds !
1172  ! f(r_1=0) from f(r_2), f(r_3) and f(r_4) !
1173  ! !
1174  !-------------------------------------------------------!
1175 
1176  USE itm_types
1177  USE itm_constants
1178 
1179  IMPLICIT NONE
1180 
1181  INTEGER :: n, i
1182  REAL (R8) :: h(n), r(n), f(n)
1183 
1184 
1185  DO i=1,3
1186  h(i)=r(i+1)-r(i)
1187  END DO
1188 
1189  f(1) = ((f(2)*r(4)/h(2)+f(4)*r(2)/h(3))*r(3) &
1190  -f(3)*(r(2)/h(2)+r(4)/h(3))*r(2)*r(4)/r(3)) &
1191  /(r(4)-r(2))
1192 
1193  RETURN
1194 
1195 
1196  END SUBROUTINE axis_int
1197  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1198  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1199 
1200 
1201 
1202 END MODULE emeq
emeq::emeq_e3m
subroutine emeq_e3m(EQUILIBRIUM_IN, EQUILIBRIUM_OUT, code_parameters)
Definition: emeq.f90:18
emeq::axis_int
subroutine axis_int(n, r, f)
Definition: emeq.f90:1168
l3deriv
subroutine l3deriv(y_in, x_in, nr_in, dydx_out, x_out, nr_out)
Definition: l3interp.f90:59
assign_code_parameters
subroutine assign_code_parameters(codeparameters, return_status)
Definition: emeq.f90:671
emeq::emeq_interface
subroutine emeq_interface(
Definition: emeq.f90:817
fun
subroutine fun(X, F)
Definition: Ev2.f:10
mod_f90_kind
types
Definition: mod_f90_kind.f90:8
emeq
EMEQ_E3M.
Definition: emeq.f90:11
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
helena21_mod::fdia
real(r8) function fdia(psi_n)
Definition: helena21_mod.f90:139
emeq_equil::dealloc_emeq_equil
subroutine, public dealloc_emeq_equil
Definition: emeq_equil.f:93
emeq_equil::alloc_emeq_equil
subroutine, public alloc_emeq_equil(NP_in)
Definition: emeq_equil.f:53
ets_version
Definition: ets_version.f90:1
emeq_equil::e3astr
subroutine, public e3astr
Definition: emeq_equil.f:178
read_codeparam
subroutine read_codeparam(in_xml, filename, codeparam)
Definition: emeq.f90:760
emeq_equil
replacement for the EMEQ common block
Definition: emeq_equil.f:8