ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
equilibrium_tools.f90
Go to the documentation of this file.
1 MODULE equilibrium_tools
2 !--------------------------------------------------------
3 ! The module contains subroutines !
4 ! used during initialization of ETS !
5 ! for defining magnetic configuration !
6 ! !
7 ! author: Denis Kalupin !
8 ! e-mail: Denis.Kalupin@euro-fusion.org !
9 ! revisions: Pär Strand, par.strand@chalmers.se !
10 !--------------------------------------------------------
11 ! revisions: !
12 ! 2016-08-17 (Pär Strand, par.strand@chalmers.se !
13 ! !
14 ! - added read protection for invalid input in !
15 ! geom_from_ cpo !
16 ! - corrected append line output_diag to maintain !
17 ! history of messages !
18 ! - added calculations from R,Z boundary !
19 ! - Modifed diagnostic error handling to avoid !
20 ! self healing of errors and to comply with !
21 ! external usage !
22 ! diag%ierr = 0 - no error !
23 ! diag%ierr = 1 - warnings !
24 ! diag%ierr = 2 - Blocking errors !
25 ! diag%ierr = -1 - Blocking errors !
26 !--------------------------------------------------------
27 CONTAINS
28 !--------------------------------------------------------
29 !--------------------------------------------------------
30 SUBROUTINE geometry_from_wf_parameters(EQUILIBRIUM, EQUILIBRIUM_OUT, &
31 !AXIS:
32  geo_ax, mag_ax, plasma_ax, &
33 !GEOMETRICAL QUANTITIES:
34  ip, amin, &
35  elong_up, elong_low, &
36  tria_up, tria_low, &
37 !SPACE_RESOLUTION:
38  npsi, ndim1, ndim2, npoints)
39 
40 !--------------------------------------------------------
41 ! This subroutine prefiles equilibrium !
42 ! variables from parameters provided in the !
43 ! workflow. !
44 ! It culates the boundary !
45 ! and RHO_TOR (at the boundary) for the ETS grid. !
46 !--------------------------------------------------------
47 
48 ! +++ Declaration of variables:
49  USE euitm_schemas
50  USE itm_types
51  USE itm_constants
52  USE ets_plasma
53 
54  IMPLICIT NONE
55 
56 ! +++ CPO derived types:
57  TYPE (type_equilibrium), POINTER :: equilibrium(:) !time independent CPO slice
58  TYPE (type_equilibrium), POINTER :: equilibrium_out(:) !time independent CPO slice
59 
60 ! +++ Local derived types:
61  TYPE (diagnostic) :: diag !diagnostic output
62 
63 ! +++ Internal derived types:
64  REAL (R8) :: geo_ax(3) !{R, Z, B} at device geometrical centre
65  REAL (R8) :: mag_ax(3) !{R, Z, B} at magnetic axis
66  REAL (R8) :: plasma_ax(3) !{R, Z, B} at LCMS centre
67 
68  REAL (R8) :: ip !Total plasma current [A]
69  REAL (R8) :: amin !minor radius [m]
70  REAL (R8) :: elong !averaged elongation oof LCMS
71  REAL (R8) :: elong_up, elong_low !upper and lower elongation oof LCMS
72  REAL (R8) :: tria !averaged triangularity oof LCMS
73  REAL (R8) :: tria_up, tria_low !upper and lower triangularity oof LCMS
74 
75  REAL (R8) :: theta !poloidal angle
76  REAL (R8) :: rho_bnd !toroidal minor radius, boundary value[m]
77  REAL (R8), ALLOCATABLE :: rho(:) !toroidal minor radius [m]
78  REAL (R8), ALLOCATABLE :: rho_2d(:) !toroidal minor radius for 2D signals [m]
79 
80  INTEGER :: npsi, ipsi !total number and index of rho grid on grid and output CPOs
81  INTEGER :: ndim1, idim1 !total number and index of first dimension of EQ 2D signals
82  INTEGER :: ndim2, idim2 !total number and index of second dimension of EQ 2D signals
83  INTEGER :: npoints, itheta !total number and index of boundary points
84  INTEGER :: nline, i
85  CHARACTER*1000 :: output_diag
86 
87 
88 
89 ! +++ Nullify diagnostic output:
90  diag%ERROR_MESSAGE = " "
91  diag%IERR = 0
92 
93 
94 ! +++ Allocate output CPO:
95  IF(ASSOCIATED(equilibrium_out)) DEALLOCATE (equilibrium_out)
96  ALLOCATE (equilibrium_out(1))
97 
98 ! +++ GLOBAL_PARAM:
99  equilibrium_out(1)%global_param%i_plasma = ip
100  equilibrium_out(1)%global_param%toroid_field%r0 = geo_ax(1)
101  equilibrium_out(1)%global_param%toroid_field%b0 = geo_ax(3)
102 
103  equilibrium_out(1)%global_param%mag_axis%position%r = mag_ax(1)
104  equilibrium_out(1)%global_param%mag_axis%position%z = mag_ax(2)
105  equilibrium_out(1)%global_param%mag_axis%bphi = mag_ax(3)
106 
107 
108 ! +++ RHO BOUNDARY:
109  CALL calculate_rho_tor(plasma_ax, mag_ax, amin, elong_up, elong_low, tria_up, tria_low, &
110  rho_bnd, diag)
111  IF (diag%IERR == -1 .OR. diag%IERR == 2 ) goto 90
112 
113  ALLOCATE (rho(npsi))
114  ALLOCATE (equilibrium_out(1)%profiles_1d%rho_tor(npsi))
115 
116  DO ipsi=1,npsi
117  rho(ipsi) = rho_bnd/(npsi-1)*(ipsi-1)
118  END DO
119  equilibrium_out(1)%profiles_1d%rho_tor = rho
120 
121 
122 
123 ! +++ PROFILES_1D:
124  ALLOCATE (equilibrium_out(1)%profiles_1d%elongation(npsi))
125  ALLOCATE (equilibrium_out(1)%profiles_1d%tria_upper(npsi))
126  ALLOCATE (equilibrium_out(1)%profiles_1d%tria_lower(npsi))
127 
128  ALLOCATE (equilibrium_out(1)%profiles_1d%r_inboard(npsi))
129  ALLOCATE (equilibrium_out(1)%profiles_1d%r_outboard(npsi))
130 
131  ALLOCATE (equilibrium_out(1)%profiles_1d%r_inboard(npsi))
132  ALLOCATE (equilibrium_out(1)%profiles_1d%r_outboard(npsi))
133  ALLOCATE (equilibrium_out(1)%profiles_1d%gm1(npsi))
134  ALLOCATE (equilibrium_out(1)%profiles_1d%gm2(npsi))
135  ALLOCATE (equilibrium_out(1)%profiles_1d%gm3(npsi))
136  ALLOCATE (equilibrium_out(1)%profiles_1d%gm4(npsi))
137  ALLOCATE (equilibrium_out(1)%profiles_1d%gm5(npsi))
138  ALLOCATE (equilibrium_out(1)%profiles_1d%gm6(npsi))
139  ALLOCATE (equilibrium_out(1)%profiles_1d%gm7(npsi))
140  ALLOCATE (equilibrium_out(1)%profiles_1d%gm8(npsi))
141  ALLOCATE (equilibrium_out(1)%profiles_1d%gm9(npsi))
142  ALLOCATE (equilibrium_out(1)%profiles_1d%volume(npsi))
143  ALLOCATE (equilibrium_out(1)%profiles_1d%vprime(npsi))
144  ALLOCATE (equilibrium_out(1)%profiles_1d%area(npsi))
145  ALLOCATE (equilibrium_out(1)%profiles_1d%aprime(npsi))
146  ALLOCATE (equilibrium_out(1)%profiles_1d%F_dia(npsi))
147  ALLOCATE (equilibrium_out(1)%profiles_1d%rho_vol(npsi))
148 
149  equilibrium_out(1)%profiles_1d%elongation(npsi) = (elong_up + elong_low) / 2._r8
150  equilibrium_out(1)%profiles_1d%tria_upper(npsi) = tria_up
151  equilibrium_out(1)%profiles_1d%tria_lower(npsi) = tria_low
152 
153  equilibrium_out(1)%profiles_1d%r_inboard(npsi) = plasma_ax(1) - amin
154  equilibrium_out(1)%profiles_1d%r_outboard(npsi) = plasma_ax(1) + amin
155  equilibrium_out(1)%profiles_1d%gm1 = 1.e0_r8 / plasma_ax(1)**2
156  equilibrium_out(1)%profiles_1d%gm2 = 1.e0_r8 / plasma_ax(1)**2
157  equilibrium_out(1)%profiles_1d%gm3 = 1.e0_r8
158  equilibrium_out(1)%profiles_1d%gm4 = 1.e0_r8 / plasma_ax(3)**2
159  equilibrium_out(1)%profiles_1d%gm5 = plasma_ax(3)**2
160  equilibrium_out(1)%profiles_1d%gm6 = 1.e0_r8 / plasma_ax(3)**2
161  equilibrium_out(1)%profiles_1d%gm7 = 1.e0_r8
162  equilibrium_out(1)%profiles_1d%gm8 = plasma_ax(1)
163  equilibrium_out(1)%profiles_1d%gm9 = 1.e0_r8 / plasma_ax(1)
164  equilibrium_out(1)%profiles_1d%volume = 2.e0_r8 * itm_pi**2 * rho**2 * plasma_ax(1)
165  equilibrium_out(1)%profiles_1d%vprime = 4.e0_r8 * itm_pi**2 * rho * plasma_ax(1)
166  equilibrium_out(1)%profiles_1d%area = itm_pi * rho**2
167  equilibrium_out(1)%profiles_1d%aprime = 4.e0_r8 * itm_pi**2 * plasma_ax(1)
168  equilibrium_out(1)%profiles_1d%F_dia = plasma_ax(3) * plasma_ax(1)
169  equilibrium_out(1)%profiles_1d%rho_vol = sqrt(equilibrium_out(1)%profiles_1d%volume/equilibrium_out(1)%profiles_1d%volume(npsi))
170 
171 
172 ! +++ EQ_GEOMETRY:
173  equilibrium_out(1)%eqgeometry%geom_axis%r = plasma_ax(1)
174  equilibrium_out(1)%eqgeometry%geom_axis%z = plasma_ax(2)
175 
176  equilibrium_out(1)%eqgeometry%a_minor = amin
177  equilibrium_out(1)%eqgeometry%elongation = (elong_up + elong_low) / 2._r8
178  equilibrium_out(1)%eqgeometry%elong_upper = elong_up
179  equilibrium_out(1)%eqgeometry%elong_lower = elong_low
180  equilibrium_out(1)%eqgeometry%tria_upper = tria_up
181  equilibrium_out(1)%eqgeometry%tria_lower = tria_low
182 
183 
184 ! +++ BOUNDARY:
185  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1))
186  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1)%r(npoints))
187  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1)%z(npoints))
188 
189  DO itheta=1, npoints
190  theta = REAL(itheta-1,r8)/REAL(npoints-1,r8)*2.0_r8*itm_pi
191  IF (theta.LE.itm_pi) tria = tria_up
192  IF (theta.LE.itm_pi) elong = elong_up
193  IF (theta.GT.itm_pi) tria = tria_low
194  IF (theta.GT.itm_pi) elong = elong_low
195  equilibrium_out(1)%eqgeometry%boundary(1)%r(itheta) = plasma_ax(1) + amin * (cos(theta) - tria*(sin(theta))**2)
196  equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta) = plasma_ax(2) + amin * elong * sin(theta)
197 
198  IF (abs(equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta)).LE.1.e-15_r8) &
199  equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta) = 0.0e0_r8
200  END DO
201 
202 
203 ! +++ COORD_SYS:
204  ALLOCATE (rho_2d(ndim1))
205  rho_loop2: DO idim1=1,ndim1
206  rho_2d(idim1) = rho_bnd/(ndim1-1)*(idim1-1)
207  END DO rho_loop2
208 
209  ALLOCATE(equilibrium_out(1)%coord_sys%grid_type(1))
210  ALLOCATE(equilibrium_out(1)%coord_sys%grid%dim1(ndim1))
211  ALLOCATE(equilibrium_out(1)%coord_sys%grid%dim2(ndim2))
212  ALLOCATE(equilibrium_out(1)%coord_sys%position%R(ndim1,ndim2))
213  ALLOCATE(equilibrium_out(1)%coord_sys%position%Z(ndim1,ndim2))
214 
215 
216 
217  equilibrium_out(1)%coord_sys%grid_type = '(rho,theta)'
218  equilibrium_out(1)%coord_sys%grid%dim1 = rho_2d
219 
220  DO idim2=1, ndim2
221  theta = REAL(idim2-1,r8)/REAL(ndim2-1,r8)*2.0_r8*itm_pi
222  equilibrium_out(1)%coord_sys%grid%dim2(idim2) = theta
223  DO idim1=1, ndim1
224 
225  IF (theta.LE.itm_pi) tria = tria_up
226  IF (theta.GT.itm_pi) tria = tria_low
227  IF (theta.LE.itm_pi) elong = elong_up
228  IF (theta.GT.itm_pi) elong = elong_low
229 
230  equilibrium_out(1)%coord_sys%position%R(idim1,idim2) = plasma_ax(1) + amin * (cos(theta) - tria*(sin(theta))**2) * rho_2d(idim1)/rho_2d(ndim1)
231  equilibrium_out(1)%coord_sys%position%Z(idim1,idim2) = plasma_ax(2) + amin * rho_2d(idim1)/rho_2d(ndim1) * elong * sin(theta)
232  IF (abs(equilibrium_out(1)%eqgeometry%boundary(1)%z(idim2)).LE.1.e-15_r8) &
233  equilibrium_out(1)%coord_sys%position%Z(idim1,idim2) = 0.0e0_r8
234  END DO
235  END DO
236 
237  DEALLOCATE (rho)
238  DEALLOCATE (rho_2d)
239 
240 ! +++ CODEPARAM:
241  90 output_diag = "GEOMETRY_FROM_WF_PARAMETERS: "//trim(adjustl(diag%ERROR_MESSAGE))
242  nline = floor(len_trim(adjustl(output_diag))/132.001)+1
243  ALLOCATE (equilibrium_out(1)%codeparam%codename(1))
244  ALLOCATE (equilibrium_out(1)%codeparam%codeversion(1))
245  ALLOCATE (equilibrium_out(1)%codeparam%output_diag(nline))
246 
247  equilibrium_out(1)%codeparam%codename = 'GEOM_FROM_WF_PARAMETERS'
248  equilibrium_out(1)%codeparam%codeversion = 'GEOM_FROM_WF_PARAMETERS_4.10b.10'
249  equilibrium_out(1)%codeparam%output_flag = diag%IERR
250  DO i = 1,nline
251  equilibrium_out(1)%codeparam%output_diag(i) = output_diag(((i-1)*132+1):(i*132))
252  END DO
253 
254 
255 RETURN
256 END SUBROUTINE geometry_from_wf_parameters
257 !--------------------------------------------------------
258 !--------------------------------------------------------
259 
260 
261 
262 !--------------------------------------------------------
263 !--------------------------------------------------------
264 SUBROUTINE geometry_from_cpo(EQUILIBRIUM, EQUILIBRIUM_OUT, &
265  !coord_sys resolution:
266  ndim1, ndim2,ierr_out,error_mes)
267 
268 !--------------------------------------------------------
269 ! This subroutine checks key !
270 ! equilibrium variables. !
271 ! If necessary it recalculates the boundary !
272 ! and RHO_TOR (at the boundary) for the ETS grid. !
273 !--------------------------------------------------------
274 
275 
276 ! +++ Declaration of variables:
277  USE euitm_schemas
278  USE itm_types
279  USE itm_constants
280  USE ets_plasma
281  USE copy_structures
282  USE deallocate_structures
283 
284  IMPLICIT NONE
285 
286 ! +++ CPO derived types:
287  TYPE (type_equilibrium), POINTER :: equilibrium(:) !time independent CPO slice
288  TYPE (type_equilibrium), POINTER :: equilibrium_out(:) !time independent CPO slice
289 
290 
291 ! +++ Local derived types:
292  TYPE (diagnostic) :: diag !diagnostic output
293 
294 
295 ! +++ Internal derived types:
296  REAL (R8) :: amin !minor radius [m]
297  REAL (R8) :: elong !averaged elongation oof LCMS
298  REAL (R8) :: elong_up, elong_low !upper and lower elongation oof LCMS
299  REAL (R8) :: tria !averaged triangularity oof LCMS
300  REAL (R8) :: tria_up, tria_low !upper and lower triangularity oof LCMS
301 
302  REAL (R8) :: geo_ax(3) !{R, Z, B} at device geometrical centre
303  REAL (R8) :: mag_ax(3) !{R, Z, B} at magnetic axis
304  REAL (R8) :: plasma_ax(3) !{R, Z, B} at LCMS centre
305 
306  REAL (R8) :: theta !poloidal angle
307  REAL (R8) :: rho_bnd !toroidal minor radius, boundary value[m]
308 
309  INTEGER :: npsi=200, ipsi !total number and index of rho grid on grid and output CPOs
310  INTEGER :: npoints=300, itheta !total number and index of boundary points
311  INTEGER :: ndim1, idim1 !total number and index of first dimension of EQ 2D signals
312  INTEGER :: ndim2, idim2 !total number and index of second dimension of EQ 2D signals
313  INTEGER :: nline, i, iwarning
314  CHARACTER*1000 :: output_diag
315  character*255 :: error_mes
316  integer :: ierr_out
317 
318 ierr_out=0
319 error_mes=''
320 
321 if (equilibrium(1)%time.eq.itm_r8_invalid) then
322  ierr_out=-1
323  error_mes='equilibrium is not present on input, you might check your basic paramters like username, shot_in, run_in, etc.'
324  allocate(equilibrium_out(1))
325  return
326 end if
327 
328 
329 ! +++ Nullify diagnostic output:
330  diag%ERROR_MESSAGE = " "
331  diag%IERR = 0
332  iwarning = 0
333 
334 ! +++ Copy input Equlibrium to output:
335  ALLOCATE (equilibrium_out(1))
336  CALL copy_cpo(equilibrium(1), equilibrium_out(1))
337 
338 
339 
340 ! +++ Check Plasma Current:
341  IF (equilibrium(1)%global_param%i_plasma.NE.equilibrium(1)%global_param%i_plasma &
342  .OR. .NOT. (itm_is_valid(equilibrium(1)%global_param%i_plasma))) THEN
343  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"PLASMA CURRENT is not present in input EQUILIBRIUM. "
344  diag%IERR = -1
345  goto 100
346  END IF
347 
348  CALL eq_param_from_cpo(equilibrium, geo_ax, mag_ax, plasma_ax,amin, &
349  elong_up, elong_low, tria_up, tria_low, diag)
350  IF (diag%IERR == -1 .OR. diag%IERR == 2) goto 100
351 
352 ! +++ Check RHO_TOR:
353  IF (.NOT.ASSOCIATED (equilibrium(1)%profiles_1d%rho_tor)) THEN
354  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"RHO_TOR is not present in input EQUILIBRIUM. "
355  iwarning = iwarning + 1
356 
357  CALL calculate_rho_tor(plasma_ax, mag_ax, amin, elong_up, elong_low, tria_up, tria_low, &
358  rho_bnd, diag)
359  IF (diag%IERR == -1 .OR. diag%IERR == 2 ) goto 100
360 
361  IF (ASSOCIATED (equilibrium_out(1)%profiles_1d%rho_tor)) THEN
362  npsi = SIZE(equilibrium_out(1)%profiles_1d%rho_tor)
363  ELSE
364  ALLOCATE (equilibrium_out(1)%profiles_1d%rho_tor(npsi))
365  END IF
366 
367  DO ipsi=1,npsi
368  equilibrium_out(1)%profiles_1d%rho_tor(ipsi) = rho_bnd/(npsi-1)*(ipsi-1)
369  END DO
370 
371  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"RHO_TOR was recalculated. "
372 
373  END IF
374 
375 
376 ! jofe> fill some of the missing fields with computed/derived quantities,
377 ! which may be needed by some of the equilibrium codes (eg. EMEQ)
378 
379 ! +++ GLOBAL_PARAM:
380  IF (equilibrium_out(1)%global_param%toroid_field%r0.NE.equilibrium_out(1)%global_param%toroid_field%r0 &
381  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%global_param%toroid_field%r0))) THEN
382  equilibrium_out(1)%global_param%toroid_field%r0 = geo_ax(1)
383  ENDIF
384  IF (equilibrium_out(1)%global_param%toroid_field%b0.NE.equilibrium_out(1)%global_param%toroid_field%b0 &
385  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%global_param%toroid_field%b0))) THEN
386  equilibrium_out(1)%global_param%toroid_field%b0 = geo_ax(3)
387  ENDIF
388  IF (equilibrium_out(1)%global_param%mag_axis%position%r.NE.equilibrium_out(1)%global_param%mag_axis%position%r &
389  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%global_param%mag_axis%position%r))) THEN
390  equilibrium_out(1)%global_param%mag_axis%position%r = mag_ax(1)
391  ENDIF
392  IF (equilibrium_out(1)%global_param%mag_axis%position%z.NE.equilibrium_out(1)%global_param%mag_axis%position%z &
393  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%global_param%mag_axis%position%z))) THEN
394  equilibrium_out(1)%global_param%mag_axis%position%z = mag_ax(2)
395  ENDIF
396  IF (equilibrium_out(1)%global_param%mag_axis%bphi.NE.equilibrium_out(1)%global_param%mag_axis%bphi &
397  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%global_param%mag_axis%bphi))) THEN
398  equilibrium_out(1)%global_param%mag_axis%bphi = mag_ax(3)
399  ENDIF
400 
401 ! +++ EQ_GEOMETRY:
402  IF (equilibrium_out(1)%eqgeometry%geom_axis%r.NE.equilibrium_out(1)%eqgeometry%geom_axis%r &
403  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%eqgeometry%geom_axis%r))) THEN
404  equilibrium_out(1)%eqgeometry%geom_axis%r = plasma_ax(1)
405  ENDIF
406  IF (equilibrium_out(1)%eqgeometry%geom_axis%z.NE.equilibrium_out(1)%eqgeometry%geom_axis%z &
407  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%eqgeometry%geom_axis%z))) THEN
408  equilibrium_out(1)%eqgeometry%geom_axis%z = plasma_ax(2)
409  ENDIF
410  IF (equilibrium_out(1)%eqgeometry%a_minor.NE.equilibrium_out(1)%eqgeometry%a_minor &
411  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%eqgeometry%a_minor))) THEN
412  equilibrium_out(1)%eqgeometry%a_minor = amin
413  ENDIF
414  IF (equilibrium_out(1)%eqgeometry%elongation.NE.equilibrium_out(1)%eqgeometry%elongation &
415  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%eqgeometry%elongation))) THEN
416  equilibrium_out(1)%eqgeometry%elongation = (elong_up + elong_low) / 2._r8
417  ENDIF
418  IF (equilibrium_out(1)%eqgeometry%elong_upper.NE.equilibrium_out(1)%eqgeometry%elong_upper &
419  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%eqgeometry%elong_upper))) THEN
420  equilibrium_out(1)%eqgeometry%elong_upper = elong_up
421  ENDIF
422  IF (equilibrium_out(1)%eqgeometry%elong_lower.NE.equilibrium_out(1)%eqgeometry%elong_lower &
423  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%eqgeometry%elong_lower))) THEN
424  equilibrium_out(1)%eqgeometry%elong_lower = elong_low
425  ENDIF
426  IF (equilibrium_out(1)%eqgeometry%tria_upper.NE.equilibrium_out(1)%eqgeometry%tria_upper &
427  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%eqgeometry%tria_upper))) THEN
428  equilibrium_out(1)%eqgeometry%tria_upper = tria_up
429  ENDIF
430  IF (equilibrium_out(1)%eqgeometry%tria_lower.NE.equilibrium_out(1)%eqgeometry%tria_lower &
431  .OR. .NOT. (itm_is_valid(equilibrium_out(1)%eqgeometry%tria_lower))) THEN
432  equilibrium_out(1)%eqgeometry%tria_lower = tria_low
433  ENDIF
434 ! jofe>
435 
436 ! +++ Check Boundary:
437 
438  IF(.NOT.ASSOCIATED (equilibrium_out(1)%eqgeometry%boundary)) &
439  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1))
440 
441  IF (.NOT.ASSOCIATED (equilibrium_out(1)%eqgeometry%boundary(1)%r) .OR. &
442  .NOT.ASSOCIATED (equilibrium_out(1)%eqgeometry%boundary(1)%z)) THEN
443 
444  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"Boundary is not defined in input EQUILIBRIUM. "
445  iwarning = iwarning + 1
446 
447 
448  IF(ASSOCIATED (equilibrium_out(1)%eqgeometry%boundary(1)%r)) THEN
449  npoints = SIZE(equilibrium_out(1)%eqgeometry%boundary(1)%r)
450  IF(.NOT.ASSOCIATED(equilibrium_out(1)%eqgeometry%boundary(1)%z)) &
451  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1)%z(npoints))
452 
453  ELSE IF (ASSOCIATED(equilibrium_out(1)%eqgeometry%boundary(1)%z) ) THEN
454  npoints = SIZE(equilibrium_out(1)%eqgeometry%boundary(1)%z)
455  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1)%r(npoints))
456 
457  ELSE
458  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1)%r(npoints))
459  ALLOCATE (equilibrium_out(1)%eqgeometry%boundary(1)%z(npoints))
460 
461  END IF
462 
463  DO itheta=1, npoints
464  theta = REAL(itheta-1,r8)/REAL(npoints-1,r8)*2.0_r8*itm_pi
465  IF (theta.LE.itm_pi) tria = tria_up
466  IF (theta.LE.itm_pi) elong = elong_up
467  IF (theta.GT.itm_pi) tria = tria_low
468  IF (theta.GT.itm_pi) elong = elong_low
469  equilibrium_out(1)%eqgeometry%boundary(1)%r(itheta) = plasma_ax(1) + amin * (cos(theta) - tria*(sin(theta))**2)
470  equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta) = plasma_ax(2) + amin * elong * sin(theta)
471 
472  IF (abs(equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta)).LE.1.e-15_r8) &
473  equilibrium_out(1)%eqgeometry%boundary(1)%z(itheta) = 0.0e0_r8
474  END DO
475 
476  END IF
477 
478 
479  ! +++ Allocate coord_sys if empty:
480 
481  IF (.not.ASSOCIATED(equilibrium_out(1)%coord_sys%position%R)) THEN
482  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%grid_type)) &
483  DEALLOCATE(equilibrium_out(1)%coord_sys%grid_type)
484  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%grid%dim1)) &
485  DEALLOCATE(equilibrium_out(1)%coord_sys%grid%dim1)
486  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%grid%dim2)) &
487  DEALLOCATE(equilibrium_out(1)%coord_sys%grid%dim2)
488  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%position%R)) &
489  DEALLOCATE(equilibrium_out(1)%coord_sys%position%R)
490  IF(ASSOCIATED(equilibrium_out(1)%coord_sys%position%Z)) &
491  DEALLOCATE(equilibrium_out(1)%coord_sys%position%Z)
492 
493  ALLOCATE(equilibrium_out(1)%coord_sys%grid_type(1))
494  ALLOCATE(equilibrium_out(1)%coord_sys%grid%dim1(ndim1))
495  ALLOCATE(equilibrium_out(1)%coord_sys%grid%dim2(ndim2))
496  ALLOCATE(equilibrium_out(1)%coord_sys%position%R(ndim1,ndim2))
497  ALLOCATE(equilibrium_out(1)%coord_sys%position%Z(ndim1,ndim2))
498 
499  equilibrium_out(1)%coord_sys%grid_type = 'empty'
500  equilibrium_out(1)%coord_sys%grid%dim1 = 0.0_r8
501  equilibrium_out(1)%coord_sys%grid%dim2 = 0.0_r8
502  equilibrium_out(1)%coord_sys%position%R = 0.0_r8
503  equilibrium_out(1)%coord_sys%position%Z = 0.0_r8
504  END IF
505 
506 
507 ! +++ Output Error messages
508 
509 ! CODEPARAM:
510 100 output_diag = "GEOMETRY_FROM_CPO: "//trim(adjustl(diag%ERROR_MESSAGE))
511  nline = floor(len_trim(adjustl(output_diag))/132.001)+1
512  IF(ASSOCIATED(equilibrium_out(1)%codeparam%codename)) &
513  DEALLOCATE(equilibrium_out(1)%codeparam%codename)
514  IF(ASSOCIATED(equilibrium_out(1)%codeparam%codeversion)) &
515  DEALLOCATE(equilibrium_out(1)%codeparam%codeversion)
516  IF(ASSOCIATED(equilibrium_out(1)%codeparam%output_diag)) &
517  DEALLOCATE(equilibrium_out(1)%codeparam%output_diag)
518  ALLOCATE (equilibrium_out(1)%codeparam%codename(1))
519  ALLOCATE (equilibrium_out(1)%codeparam%codeversion(1))
520  ALLOCATE (equilibrium_out(1)%codeparam%output_diag(nline))
521 
522  equilibrium_out(1)%codeparam%codename = 'GEOM_FROM_CPO'
523  equilibrium_out(1)%codeparam%codeversion = 'GEOM_FROM_CPO_4.10b.10'
524 
525  IF (diag%IERR == 0 .AND. iwarning /= 0 ) diag%IERR = iwarning
526 
527  equilibrium_out(1)%codeparam%output_flag = diag%IERR
528  DO i = 1,nline
529  equilibrium_out(1)%codeparam%output_diag(i) = output_diag(((i-1)*132+1):(i*132))
530  END DO
531 
532 
533 RETURN
534 END SUBROUTINE geometry_from_cpo
535 !--------------------------------------------------------
536 !--------------------------------------------------------
537 
538 
539 
540 !--------------------------------------------------------
541 !--------------------------------------------------------
542 SUBROUTINE calculate_rho_tor(PLASMA_AX, MAG_AX, &
543  amin, &
544  elong_up, elong_low, &
545  tria_up, tria_low, &
546 ! OUTPUT:
547  rho, diag)
548 !-------------------------------------------------------!
549 ! This routine calculates RHO_TOR based !
550 ! on three moment description of the boundary !
551 ! !
552 !-------------------------------------------------------!
553 ! Source: ASTRA !
554 ! Developers: D.Kalupin !
555 ! Contacts: Denis.Kalupin@efda.org !
556 ! !
557 ! Comments: created for the ETS workflow !
558 ! !
559 !-------------------------------------------------------!
560 
561  USE itm_constants
562  USE ets_plasma
563  USE itm_types
564 
565  IMPLICIT NONE
566 
567 ! +++ Local derived types:
568  TYPE (diagnostic) :: diag !diagnostic output
569 
570 ! +++ Input/Output variables:
571  REAL (R8) :: plasma_ax(3)! PLASMA AXIS: {R[m],Z[m],B[T]}
572  REAL (R8) :: mag_ax(3) ! MAGNETIC AXIS: {R[m],Z[m],B[T]}
573  REAL (R8) :: amin ! minor radius a[m]
574  REAL (R8) :: elong_up ! elongation upper
575  REAL (R8) :: elong_low ! elongation lower
576  REAL (R8) :: tria_up ! triangularity upper
577  REAL (R8) :: tria_low ! triangularity lower
578 
579  REAL (R8) :: rho ! toroidal minor radius [m]
580 
581 ! +++ Internal (ASTRA) variables:
582  REAL (R8) :: yr ! Major radius [m]
583  REAL (R8) :: yd ! Shafranov shift [m]
584  REAL (R8) :: ya ! Minor radius [m]
585  REAL (R8) :: ye ! Elongation [d/l]
586  REAL (R8) :: yt ! Triangularity [d/l]
587 
588 
589  REAL (R8) :: yd1,yt2,yge
590  REAL (R8) :: y1,y2,y3,y4,y5,y6
591 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
592 
593 
594 
595 ! +++ Copy values to ASTRA variables:
596  yr = plasma_ax(1)
597  yd = mag_ax(1) - plasma_ax(1)
598  ya = amin
599  ye = (elong_up + elong_low) / 2._r8
600  yt = (tria_up + tria_low) / 2._r8
601 
602 
603 
604 
605 ! +++ ASTRA routine:
606  yt2 = 2._r8*yt
607  yge = ya/(yr+yd)
608  IF (abs(yge) .GT. 1._r8) goto 8
609  IF (yge) 9,5,1
610 
611 
612  1 CONTINUE
613  yd1 = 1._r8+yt2*(yt2-2._r8/yge)
614  IF (yd1) 3,2,2
615 
616 
617  2 CONTINUE
618  yd1 = sqrt(yd1)
619  y1 = (2._r8-yt2*yge)*yd1/(1._r8+yd1)
620  y2 = yge-yt2
621  y3 = sqrt(1._r8-yt2*yge+yge*yd1)
622  y4 = sqrt(1._r8-yt2*yge-yge*yd1)
623  y5 = sqrt(1._r8+yge)
624  y6 = sqrt(1._r8-yge)
625  y1 = (y1+y2)/(y3+(1._r8-yt2)*y5)-(y1-y2)/(y4+(1._r8+yt2)*y6)
626  y1 = y1/sqrt(2._r8*(1._r8-yt2*yge+y5*y6))
627  goto 4
628 
629 
630  3 CONTINUE
631  y5 = sqrt(1._r8+yge)
632  y6 = sqrt(1._r8-yge)
633  y1 = (1._r8-yt2)*y5+(1._r8+yt2)*y6
634  y1 = (1._r8-y1/sqrt(2._r8*(1._r8-yt2*yge+y5*y6)))/yt2
635 
636 
637 ! +++ OUTPUT:
638  4 CONTINUE
639  rho = 2._r8*sqrt(yr*ya*ye*y1)
640  goto 100
641 
642 
643  5 CONTINUE
644  rho = 0._r8
645  goto 100
646 
647 
648  8 diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"CALCULATE_RHO_TOR >>> Illegal input: R+Delta < a. "
649  diag%IERR = -1
650  goto 100
651 
652 
653  9 diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"CALCULATE_RHO_TOR >>> Illegal input. "
654  diag%IERR = -1
655 
656 
657  100 CONTINUE
658 
659 
660 RETURN
661 END SUBROUTINE calculate_rho_tor
662 !--------------------------------------------------------
663 !--------------------------------------------------------
664 
665 
666 
667 !--------------------------------------------------------
668 !--------------------------------------------------------
669 SUBROUTINE eq_param_from_cpo(EQUILIBRIUM, &
670 !AXIS:
671  geo_ax, mag_ax, plasma_ax, &
672 !GEOMETRICAL QUANTITIES:
673  amin, &
674  elong_up, elong_low, &
675  tria_up, tria_low, &
676  diag)
677 
678 !--------------------------------------------------------
679 ! This subroutine checks key !
680 ! equilibrium variables. !
681 !--------------------------------------------------------
682 
683 
684 ! +++ Declaration of variables:
685  USE euitm_schemas
686  USE itm_types
687  USE ets_plasma
688  USE copy_structures
689 
690  IMPLICIT NONE
691 
692 ! +++ CPO derived types:
693  TYPE (type_equilibrium), POINTER :: equilibrium(:) !time independent CPO slice
694 
695 
696 ! +++ Local derived types:
697  TYPE (diagnostic) :: diag !diagnostic output
698 
699 
700 ! +++ Internal derived types:
701  REAL (R8) :: r_0, z_0, b_0
702  REAL (R8) :: r_mag, z_mag, b_mag
703  REAL (R8) :: r_pla, z_pla, b_pla
704 
705  REAL (R8) :: amin
706  REAL (R8) :: elong
707  REAL (R8) :: elong_up, elong_low
708  REAL (R8) :: tria_up, tria_low
709  REAL (R8) :: rho_bnd
710 
711  REAL (R8) :: geo_ax(3)
712  REAL (R8) :: mag_ax(3)
713  REAL (R8) :: plasma_ax(3)
714 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
715  REAL (R8) :: amin_tmp
716  REAL (R8) :: rgeo_tmp
717  REAL (R8) :: zgeo_tmp
718  REAL (R8) :: tria_upper_tmp
719  REAL (R8) :: tria_lower_tmp
720  REAL (R8) :: tria_tmp
721  REAL (R8) :: elong_upper_tmp
722  REAL (R8) :: elong_lower_tmp
723  REAL (R8) :: elong_tmp
724  LOGICAL :: l_have_boundary
725  INTEGER :: iwarning, ierror
726 
727 ! +++ Global parameters:
728  l_have_boundary = .false.
729  ierror = 0
730  iwarning = 0
731 
732 
733  IF(.NOT. ASSOCIATED (equilibrium(1)%eqgeometry%boundary)) &
734  ALLOCATE (equilibrium(1)%eqgeometry%boundary(1))
735 
736  IF (ASSOCIATED (equilibrium(1)%eqgeometry%boundary(1)%r) &
737  .AND. &
738  ASSOCIATED (equilibrium(1)%eqgeometry%boundary(1)%z) ) THEN
739 
740  amin_tmp = (maxval(equilibrium(1)%eqgeometry%boundary(1)%r) - &
741  minval(equilibrium(1)%eqgeometry%boundary(1)%r)) / 2._r8
742 
743  rgeo_tmp = (maxval(equilibrium(1)%eqgeometry%boundary(1)%r) + &
744  minval(equilibrium(1)%eqgeometry%boundary(1)%r)) / 2._r8
745 
746  zgeo_tmp = (maxval(equilibrium(1)%eqgeometry%boundary(1)%z) + &
747  minval(equilibrium(1)%eqgeometry%boundary(1)%z)) / 2._r8
748 
749  tria_upper_tmp = (equilibrium(1)%eqgeometry%boundary(1)%r(maxloc(equilibrium(1)%eqgeometry%boundary(1)%z,dim=1)) &
750  - rgeo_tmp)/2.0_r8
751  tria_lower_tmp = (equilibrium(1)%eqgeometry%boundary(1)%r(minloc(equilibrium(1)%eqgeometry%boundary(1)%z,dim=1)) &
752  - rgeo_tmp)/2.0_r8
753  tria_tmp =(tria_upper_tmp + tria_lower_tmp)/2.0_r8
754 
755  elong_upper_tmp = (maxval(equilibrium(1)%eqgeometry%boundary(1)%z) - zgeo_tmp)/amin_tmp
756  elong_lower_tmp = (zgeo_tmp-minval(equilibrium(1)%eqgeometry%boundary(1)%z))/amin_tmp
757  elong_tmp = (elong_upper_tmp+elong_lower_tmp)/2.0_r8
758  l_have_boundary =.true.
759  END IF
760 
761 
762 
763 ! Minor radius
764  IF (equilibrium(1)%eqgeometry%a_minor.EQ.equilibrium(1)%eqgeometry%a_minor &
765  .AND. itm_is_valid(equilibrium(1)%eqgeometry%a_minor)) THEN
766  amin = equilibrium(1)%eqgeometry%a_minor
767  ELSE
768  IF (ASSOCIATED (equilibrium(1)%profiles_1d%r_inboard).AND. &
769  ASSOCIATED (equilibrium(1)%profiles_1d%r_outboard)) THEN
770  amin = (equilibrium(1)%profiles_1d%r_outboard(SIZE(equilibrium(1)%profiles_1d%r_outboard)) - &
771  equilibrium(1)%profiles_1d%r_inboard(SIZE(equilibrium(1)%profiles_1d%r_inboard))) / 2._r8
772  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"MINOR RADIUS is approximated from other signals. "
773  iwarning = iwarning + 1
774  ELSE
775  IF (ASSOCIATED (equilibrium(1)%eqgeometry%boundary(1)%r)) THEN
776  amin = (maxval(equilibrium(1)%eqgeometry%boundary(1)%r) - &
777  minval(equilibrium(1)%eqgeometry%boundary(1)%r)) / 2._r8
778  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"MINOR RADIUS is approximated from boundary. "
779  iwarning = iwarning + 1
780  ELSE
781  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"MINOR RADIUS is not present in input EQUILIBRIUM. "
782  ierror = ierror + 1
783  END IF
784  END IF
785  END IF
786 
787 
788 ! Elongation
789  IF (equilibrium(1)%eqgeometry%elongation.EQ.equilibrium(1)%eqgeometry%elongation &
790  .AND. itm_is_valid(equilibrium(1)%eqgeometry%elongation)) THEN
791  elong = equilibrium(1)%eqgeometry%elongation
792  ELSE
793  IF (ASSOCIATED (equilibrium(1)%profiles_1d%elongation)) THEN
794  elong = equilibrium(1)%profiles_1d%elongation(SIZE(equilibrium(1)%profiles_1d%elongation))
795  ELSE
796  IF (l_have_boundary) THEN
797  elong = elong_tmp
798  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"ELONGATION is approximated from bounday. "
799  iwarning = iwarning + 1
800  ELSE
801  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"ELONGATION is not present in input EQUILIBRIUM. "
802  ierror = ierror + 1
803  END IF
804  END IF
805  END IF
806 
807 
808 ! Elongation Upper
809  IF (equilibrium(1)%eqgeometry%elong_upper.EQ.equilibrium(1)%eqgeometry%elong_upper &
810  .AND. itm_is_valid(equilibrium(1)%eqgeometry%elong_upper)) THEN
811  elong_up = equilibrium(1)%eqgeometry%elong_upper
812  ELSE
813  IF (ASSOCIATED (equilibrium(1)%profiles_1d%elongation)) THEN
814  elong_up = equilibrium(1)%profiles_1d%elongation(SIZE(equilibrium(1)%profiles_1d%elongation))
815  ELSE
816  IF (l_have_boundary) THEN
817  elong_up = elong_upper_tmp
818  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"ELONG_UPPER is approximated from boundary. "
819  iwarning = iwarning + 1
820  ELSE
821  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"ELONG_UPPER is not present in input EQUILIBRIUM. "
822  ierror = ierror + 1
823  END IF
824  END IF
825  END IF
826 
827 
828 ! Elongation Lower
829  IF (equilibrium(1)%eqgeometry%elong_lower.EQ.equilibrium(1)%eqgeometry%elong_lower &
830  .AND. itm_is_valid(equilibrium(1)%eqgeometry%elong_lower) ) THEN
831  elong_low = equilibrium(1)%eqgeometry%elong_lower
832  ELSE
833  IF (ASSOCIATED (equilibrium(1)%profiles_1d%elongation)) THEN
834  elong_low = equilibrium(1)%profiles_1d%elongation(SIZE(equilibrium(1)%profiles_1d%elongation))
835  ELSE
836  IF (l_have_boundary) THEN
837  elong_low = elong_lower_tmp
838  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"ELONG_LOWER is approximated from boundary. "
839  iwarning = iwarning +1
840  ELSE
841  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"ELONG_LOWER is not present in input EQUILIBRIUM. "
842  ierror = ierror + 1
843  ENDIF
844  END IF
845  END IF
846 
847 
848 ! Triangularity Upper
849  IF (equilibrium(1)%eqgeometry%tria_upper.EQ.equilibrium(1)%eqgeometry%tria_upper &
850  .AND. itm_is_valid(equilibrium(1)%eqgeometry%tria_upper)) THEN
851  tria_up = equilibrium(1)%eqgeometry%tria_upper
852  ELSE
853  IF (ASSOCIATED (equilibrium(1)%profiles_1d%tria_upper)) THEN
854  tria_up= equilibrium(1)%profiles_1d%tria_upper(SIZE(equilibrium(1)%profiles_1d%tria_upper))
855  ELSE
856  IF (l_have_boundary) THEN
857  tria_up = tria_upper_tmp
858  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"TRIANG_UPPER is approximated from bounday. "
859  iwarning = iwarning + 1
860  ELSE
861  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"TRIANG_UPPER is not present in input EQUILIBRIUM. "
862  ierror = ierror + 1
863  END IF
864  END IF
865  END IF
866 
867 
868 ! Triangularity Lower
869  IF (equilibrium(1)%eqgeometry%tria_lower.EQ.equilibrium(1)%eqgeometry%tria_lower &
870  .AND. itm_is_valid(equilibrium(1)%eqgeometry%tria_lower)) THEN
871  tria_low = equilibrium(1)%eqgeometry%tria_lower
872  ELSE
873  IF (ASSOCIATED (equilibrium(1)%profiles_1d%tria_lower)) THEN
874  tria_low= equilibrium(1)%profiles_1d%tria_lower(SIZE(equilibrium(1)%profiles_1d%tria_lower))
875  ELSE
876  IF (l_have_boundary) THEN
877  tria_low = tria_lower_tmp
878  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"TRIANG_LOWER is approximated from boundary. "
879  iwarning = iwarning + 1
880  ELSE
881  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"TRIANG_LOWER is not present in input EQUILIBRIUM. "
882  ierror = ierror + 1
883  END IF
884  END IF
885  END IF
886 
887 
888 ! +++ Axis:
889 
890 ! GEOMETRICAL CENTRE OF THE DEVICE
891  IF (equilibrium(1)%global_param%toroid_field%r0.EQ.equilibrium(1)%global_param%toroid_field%r0 &
892  .AND. itm_is_valid(equilibrium(1)%global_param%toroid_field%r0)) THEN
893  r_0 = equilibrium(1)%global_param%toroid_field%r0
894  ELSE
895  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"R_0 is not present in input EQUILIBRIUM. "
896  iwarning = iwarning + 1
897  IF (ASSOCIATED (equilibrium(1)%profiles_1d%r_inboard).AND. &
898  ASSOCIATED (equilibrium(1)%profiles_1d%r_outboard)) THEN
899  r_0 = (equilibrium(1)%profiles_1d%r_inboard(SIZE(equilibrium(1)%profiles_1d%r_inboard)) + &
900  equilibrium(1)%profiles_1d%r_outboard(SIZE(equilibrium(1)%profiles_1d%r_outboard))) / 2._r8
901  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"R_0 is calculated from other sources. "
902  iwarning = iwarning + 1
903  ELSE
904  IF (ASSOCIATED (equilibrium(1)%eqgeometry%boundary(1)%r)) THEN
905  r_0 = (maxval(equilibrium(1)%eqgeometry%boundary(1)%r) + &
906  minval(equilibrium(1)%eqgeometry%boundary(1)%r)) / 2._r8
907  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"R_0 is calculated from other sources. "
908  iwarning = iwarning + 1
909  ELSE
910  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"R_0 can not be defined from other signals. "
911  ierror = ierror + 1
912  END IF
913  END IF
914  END IF
915 
916 
917  IF (equilibrium(1)%global_param%toroid_field%b0.EQ.equilibrium(1)%global_param%toroid_field%b0 &
918  .AND. itm_is_valid(equilibrium(1)%global_param%toroid_field%b0)) THEN
919  b_0 = equilibrium(1)%global_param%toroid_field%b0
920  ELSE
921  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"B_0 is not present in input EQUILIBRIUM. "
922  ierror = ierror + 1
923  END IF
924 
925  z_0 = 0._r8
926 
927 ! PLASMA AXIS
928  IF (equilibrium(1)%eqgeometry%geom_axis%r.EQ.equilibrium(1)%eqgeometry%geom_axis%r &
929  .AND. itm_is_valid(equilibrium(1)%eqgeometry%geom_axis%r)) THEN
930  r_pla = equilibrium(1)%eqgeometry%geom_axis%r
931  ELSE
932 !jofe R_PLA = R_0
933  r_pla = rgeo_tmp
934  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"R_plasma is approximated by Rgeo. "
935  iwarning = iwarning + 1
936  END IF
937 
938  IF (equilibrium(1)%eqgeometry%geom_axis%z.EQ.equilibrium(1)%eqgeometry%geom_axis%z &
939  .AND. itm_is_valid(equilibrium(1)%eqgeometry%geom_axis%z)) THEN
940  z_pla = equilibrium(1)%eqgeometry%geom_axis%z
941  ELSE
942 !jofe Z_PLA = Z_0
943  z_pla = zgeo_tmp
944  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"Z_plasma is set to Zgeo "
945  iwarning = iwarning + 1
946  END IF
947 
948  b_pla = b_0 * r_0 / r_pla
949 
950 ! MAGNETIC AXIS
951  IF (equilibrium(1)%global_param%mag_axis%position%r.EQ.equilibrium(1)%global_param%mag_axis%position%r &
952  .AND. itm_is_valid(equilibrium(1)%global_param%mag_axis%position%r)) THEN
953  r_mag = equilibrium(1)%global_param%mag_axis%position%r
954  ELSE
955  r_mag = r_pla
956  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"R_mag is set to R_plasma. "
957  iwarning = iwarning + 1
958  END IF
959 
960  IF (equilibrium(1)%global_param%mag_axis%position%z.EQ.equilibrium(1)%global_param%mag_axis%position%z &
961  .AND. itm_is_valid(equilibrium(1)%global_param%mag_axis%position%z)) THEN
962  z_mag = equilibrium(1)%global_param%mag_axis%position%z
963  ELSE
964  z_mag = z_pla
965  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"Z_mag is set to Z_plasma. "
966  iwarning = iwarning + 1
967  END IF
968 
969  IF (equilibrium(1)%global_param%mag_axis%bphi.EQ.equilibrium(1)%global_param%mag_axis%bphi &
970  .AND. itm_is_valid(equilibrium(1)%global_param%mag_axis%bphi)) THEN
971  b_mag = equilibrium(1)%global_param%mag_axis%bphi
972  ELSE
973  b_mag = b_pla * r_pla / r_mag
974  diag%ERROR_MESSAGE = trim(adjustl(diag%ERROR_MESSAGE))//"B_mag is recalculated from B_plasma* R_plasma/R_magM. "
975  iwarning = iwarning + 1
976  END IF
977 
978  geo_ax(1) = r_0
979  geo_ax(2) = z_0
980  geo_ax(3) = b_0
981 
982  mag_ax(1) = r_mag
983  mag_ax(2) = z_mag
984  mag_ax(3) = b_mag
985 
986  plasma_ax(1) = r_pla
987  plasma_ax(2) = z_pla
988  plasma_ax(3) = b_pla
989 
990  IF (iwarning == 0 .and. ierror == 0 ) THEN
991  diag%IERR = 0
992  ELSE
993  IF (iwarning > 0 ) diag%IERR = 1
994  IF (ierror > 0 ) diag%IERR = 2
995  END IF
996 
997 RETURN
998 END SUBROUTINE eq_param_from_cpo
999 !--------------------------------------------------------
1000 !--------------------------------------------------------
1001 
1002 
1003 
1004 
1005 !--------------------------------------------------------
1006 !--------------------------------------------------------
1007 END MODULE equilibrium_tools
equilibrium_tools::geometry_from_wf_parameters
subroutine geometry_from_wf_parameters(EQUILIBRIUM,EQUILIBRIUM_OUT,
Definition: equilibrium_tools.f90:30
equilibrium_tools
Definition: equilibrium_tools.f90:1
ets_plasma::diagnostic
Definition: ets_plasma.f90:463
equilibrium_tools::calculate_rho_tor
subroutine calculate_rho_tor(PLASMA_AX, MAG_AX,AMIN,ELONG_UP,ELONG_LOW,TRIA_UP,TRIA_LOW,
Definition: equilibrium_tools.f90:542
equilibrium_tools::geometry_from_cpo
subroutine geometry_from_cpo(EQUILIBRIUM,EQUILIBRIUM_OUT,
Definition: equilibrium_tools.f90:264
ets_plasma
The module declares types of variables used in ETS (transport code)
Definition: ets_plasma.f90:8
equilibrium_tools::eq_param_from_cpo
subroutine eq_param_from_cpo(EQUILIBRIUM,
Definition: equilibrium_tools.f90:669