ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
equilibrium_work.f90
Go to the documentation of this file.
1 MODULE equilibrium_work
2 
3 CONTAINS
4 
5 
6 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
8 
9  SUBROUTINE parabolic_prof (COREPROF_IN, EQUILIBRIUM_IN, &
10  coreprof_out, equilibrium_out)
11 
12 !-------------------------------------------------------!
13 ! This routine puts parabolic !
14 ! profiles of psi, q and jparallel !
15 ! !
16 ! information received in: COREPROF_IN !
17 ! !
18 ! information saved in: COREPROF !
19 ! !
20 ! !
21 !-------------------------------------------------------!
22 ! Source: --- !
23 ! Developers: D.Kalupin !
24 ! Kontacts: Denis.Kalupin@efda.org !
25 ! !
26 ! Comments: created for V&V between ETS, JETTO !
27 ! and ASTRA !
28 ! !
29 !-------------------------------------------------------!
30 
31 
32  USE euitm_schemas
33  USE euitm_routines
34  USE copy_structures
35  USE itm_constants
36 
37  IMPLICIT NONE
38 
39  INTEGER :: ifail
40 
41 ! +++ CPO derived types:
42  TYPE (type_coreprof), POINTER :: coreprof_in(:) !input/output CPO with plasma profiles
43  TYPE (type_coreprof), POINTER :: coreprof_out(:) !input/output CPO with plasma profiles
44  TYPE (type_equilibrium), POINTER :: equilibrium_in(:) !input CPO with geometry quantities from previous time
45  TYPE (type_equilibrium), POINTER :: equilibrium_out(:) !input CPO with geometry quantities from previous time
46 
47 
48 ! +++ Dimensions:
49  INTEGER :: nrho, irho !number and index of radial points
50  INTEGER :: nion
51 
52 
53 ! +++ Local variables:
54  REAL (R8), ALLOCATABLE :: rho(:)
55  REAL (R8), ALLOCATABLE :: rhomod(:)
56 
57  REAL (R8), ALLOCATABLE :: psi(:)
58  REAL (R8), ALLOCATABLE :: qsf(:)
59  REAL (R8), ALLOCATABLE :: jpar(:)
60  REAL (R8), ALLOCATABLE :: pr(:)
61 
62  REAL (R8), ALLOCATABLE :: ne(:)
63  REAL (R8), ALLOCATABLE :: te(:)
64  REAL (R8), ALLOCATABLE :: ni(:,:)
65  REAL (R8), ALLOCATABLE :: ti(:,:)
66 
67  REAL (R8), ALLOCATABLE :: fun(:)
68  REAL (R8), ALLOCATABLE :: intjpar(:)
69  REAL (R8), ALLOCATABLE :: intfun(:)
70 
71  REAL (R8), ALLOCATABLE :: fdia(:)
72  REAL (R8), ALLOCATABLE :: gm2(:)
73  REAL (R8), ALLOCATABLE :: vprime(:)
74 
75  REAL (R8) :: curr, curr_total
76  REAL (R8) :: jpar0, jparb
77  REAL (R8) :: pr0, prb
78 
79  REAL (R8) :: r0, b0
80 
81  REAL (R8) :: bnd_value(3)
82  INTEGER :: bnd_type
83 
84  INTEGER :: nx
85  INTEGER :: npsi
86 
87  REAL (R8) :: s_ip, s_bt, s_q, s_psi, s_jpar
88 
89 
90 ! +++ Allocation:
91  ALLOCATE (coreprof_out(1))
92  ALLOCATE (equilibrium_out(1))
93 
94 
95 
96 ! +++ Fill the iterration CPO with parameters obtained on input:
97  CALL copy_cpo(coreprof_in(1), coreprof_out(1))
98  CALL copy_cpo(equilibrium_in(1), equilibrium_out(1))
99 
100 
101 
102 ! +++ Set dimensions and allocate parameters:
103  nrho = SIZE (coreprof_in(1)%rho_tor, dim=1)
104  nion = SIZE (coreprof_in(1)%ni%value, dim=2)
105  npsi = SIZE (equilibrium_in(1)%profiles_1d%psi, dim=1)
106 
107  s_ip = sign(1.0_r8, coreprof_in(1)%globalparam%current_tot)
108  s_bt = sign(1.0_r8, coreprof_in(1)%toroid_field%b0)
109  s_q = sign(1.0_r8, coreprof_in(1)%profiles1d%q%value(1))
110  s_psi = sign(1.0_r8, coreprof_in(1)%psi%value(nrho)-coreprof_in(1)%psi%value(1))
111  s_jpar = sign(1.0_r8, coreprof_in(1)%profiles1d%jtot%value(1))
112 
113  WRITE(*,*) 's_ip, s_bt, s_q, s_psi, s_jpar = ', s_ip, s_bt, s_q, s_psi, s_jpar
114 
115 ! +++ Allocation of local variables:
116  ALLOCATE (rho(nrho) )
117  ALLOCATE (rhomod(nrho) )
118 
119  ALLOCATE (psi(nrho) )
120  ALLOCATE (qsf(nrho) )
121  ALLOCATE (jpar(nrho) )
122  ALLOCATE (pr(nrho) )
123 
124  ALLOCATE (ne(nrho) )
125  ALLOCATE (te(nrho) )
126  ALLOCATE (ni(nrho,nion) )
127  ALLOCATE (ti(nrho,nion) )
128 
129  ALLOCATE (fun(nrho) )
130  ALLOCATE (intjpar(nrho) )
131  ALLOCATE (intfun(nrho))
132 
133  ALLOCATE (fdia(nrho) )
134  ALLOCATE (gm2(nrho) )
135  ALLOCATE (vprime(nrho) )
136 
137 
138 
139 ! +++ Local variables:
140  rho = coreprof_in(1)%rho_tor
141  curr_total = coreprof_in(1)%globalparam%current_tot
142 
143  r0 = coreprof_in(1)%toroid_field%r0
144  b0 = coreprof_in(1)%toroid_field%b0
145 
146  CALL l3deriv(equilibrium_in(1)%profiles_1d%volume, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
147  vprime, rho, nrho)
148  CALL l3interp(equilibrium_in(1)%profiles_1d%F_dia, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
149  fdia, rho, nrho)
150  CALL l3interp(equilibrium_in(1)%profiles_1d%gm2, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
151  gm2, rho, nrho)
152  rhomod = (1.e0_r8-rho**2 / rho(nrho)**2)
153 
154  ne = coreprof_in(1)%ne%value
155  te = coreprof_in(1)%te%value
156  ni = coreprof_in(1)%ni%value
157  ti = coreprof_in(1)%ti%value
158 
159 
160 
161 ! +++ Central and edge values for parabolic profiles:
162  jpar0 = 1.e6_r8 * s_jpar !Jpar [A/m^2]
163  jparb = 0.e0_r8
164 
165 
166 
167 ! +++ Profiles of PSI, Q and CURRENT:
168 7 jpar = (jpar0-jparb) * rhomod + jparb !Jpar [A/m^2]
169 
170 
171 ! parabolic profile for Jpar:
172  jpar = (jpar0-jparb) * rhomod + jparb
173 
174 
175 ! Recalculated profile of QSF:
176  fun = jpar * itm_mu0 * vprime / fdia**2 / rho
177  fun(1) = fun(2)
178  CALL integral(nrho, rho, fun, intfun)
179  qsf = - vprime * gm2 / fdia * rho / intfun
180  qsf(1) = qsf(2)
181 
182 
183 ! Recalculated profile of PSI:
184  fun = 2.e0_r8 * itm_pi * b0 / qsf
185  fun(1) = fun(2)
186  CALL integral(nrho, rho, fun, psi)
187 
188 
189 
190 !+++ Renormalize current:
191  fun = vprime * jpar / 2.0e0_r8 / itm_pi * b0 / fdia**2
192  CALL integral_value(nrho, rho, fun, intjpar)
193  curr = intjpar(nrho) * fdia(nrho)
194 
195  IF (dabs(1.0_r8 - curr/curr_total) .GE. 1.0e-5_r8) THEN
196 
197  jpar0 = jpar0 * curr_total / curr
198  goto 7
199  END IF
200 
201 
202 
203  WRITE (*,*) 'Parabolic current set to', curr
204 
205 
206 ! +++ Parabolic pressure profile:
207  pr = (ne * te + sum(ni * ti, dim=2)) * itm_ev !Pr [Pa]
208 
209  pr0 = pr(1)
210  prb = 0.e0_r8
211 
212  pr = (pr(1)-pr(nrho)) * rhomod + pr(nrho)
213 
214 
215 
216 ! +++ Output profiles in CPO
217  coreprof_out(1)%psi%value = psi
218  coreprof_out(1)%profiles1d%q%value = qsf
219  coreprof_out(1)%profiles1d%jtot%value = jpar
220  coreprof_out(1)%profiles1d%pr_th%value = pr
221  coreprof_out(1)%profiles1d%pr_perp%value = pr
222  coreprof_out(1)%profiles1d%pr_parallel%value = pr
223 
224 
225 ! +++ Deallocation of local variables:
226  DEALLOCATE (rho)
227  DEALLOCATE (rhomod)
228 
229  DEALLOCATE (psi)
230  DEALLOCATE (qsf)
231  DEALLOCATE (jpar)
232  DEALLOCATE (pr)
233 
234  DEALLOCATE (ne)
235  DEALLOCATE (te)
236  DEALLOCATE (ni)
237  DEALLOCATE (ti)
238 
239  DEALLOCATE (fun)
240  DEALLOCATE (intjpar)
241  DEALLOCATE (intfun)
242 
243  DEALLOCATE (fdia)
244  DEALLOCATE (gm2)
245  DEALLOCATE (vprime)
246 
247 
248 
249  RETURN
250 
251  END SUBROUTINE parabolic_prof
252 
253 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
254 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
255 
256 
257 
258 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
259 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
260 
261  SUBROUTINE readjust_profiles &
262  (prof_flag, q0_flag, coreprof_in, equilibrium_in, coreprof_out)
263 
264 !-------------------------------------------------------!
265 ! This routine readjusts consistently !
266 ! profiles of psi, q and jparallel !
267 ! !
268 ! information received in: COREPROF !
269 ! EQUILIBRIUM. !
270 ! !
271 ! information saved in: COREPROF !
272 ! !
273 ! controling parameter: PROF_FLAG !
274 ! !
275 !-------------------------------------------------------!
276 ! Source: --- !
277 ! Developers: D.Kalupin !
278 ! Kontacts: Denis.Kalupin@efda.org !
279 ! !
280 ! Comments: created for V&V between ETS, JETTO !
281 ! and ASTRA !
282 ! !
283 !-------------------------------------------------------!
284 
285 
286  USE euitm_schemas
287  USE itm_constants
288  USE copy_structures
289 
290  IMPLICIT NONE
291 
292  INTEGER :: ifail
293 
294 ! +++ CPO derived types:
295  TYPE (type_equilibrium), POINTER :: equilibrium_in(:) !input CPO with geometry quantities from previous time
296  TYPE (type_coreprof), POINTER :: coreprof_in(:) !input/output CPO with plasma profiles
297  TYPE (type_coreprof), POINTER :: coreprof_out(:) !input/output CPO with plasma profiles
298 
299 
300 ! +++ Dimensions:
301  INTEGER :: nrho, irho !number and index of radial points
302 
303  INTEGER :: prof_flag, flag !determines primary profile
304  INTEGER :: q0_flag !Flag for positive dq/drho: 0-allowed, >0-cut off
305  INTEGER :: npsi
306 
307 
308 
309 
310 ! +++ Local variables:
311  REAL (R8), ALLOCATABLE :: rho(:)
312  REAL (R8), ALLOCATABLE :: rhomod(:)
313  REAL (R8), ALLOCATABLE :: psi(:)
314  REAL (R8), ALLOCATABLE :: dpsi(:)
315  REAL (R8), ALLOCATABLE :: qsf(:)
316  REAL (R8), ALLOCATABLE :: jpar(:)
317 
318  REAL (R8), ALLOCATABLE :: fun(:)
319  REAL (R8), ALLOCATABLE :: dfun(:)
320  REAL (R8), ALLOCATABLE :: intfun(:)
321 
322  REAL (R8), ALLOCATABLE :: fdia(:)
323  REAL (R8), ALLOCATABLE :: gm2(:)
324  REAL (R8), ALLOCATABLE :: vprime(:)
325 
326  REAL (R8) :: curr_total
327  REAL (R8) :: psi0, psib
328  REAL (R8) :: qsf0, qsfb
329  REAL (R8) :: jpar0, jparb
330  REAL (R8) :: pr0, prb
331 
332  REAL (R8) :: r0, b0
333  REAL (R8) :: rho_min, qsf_min
334 
335  REAL (R8) :: s_ip, s_bt, s_q, s_psi, s_jpar
336 
337 
338 ! +++ Set dimensions and parameters:
339  nrho = SIZE (coreprof_in(1)%rho_tor, dim=1)
340  npsi = SIZE (equilibrium_in(1)%profiles_1d%psi, dim=1)
341  flag = prof_flag
342 
343  s_ip = sign(1.0_r8, coreprof_in(1)%globalparam%current_tot)
344  s_bt = sign(1.0_r8, coreprof_in(1)%toroid_field%b0)
345  s_q = sign(1.0_r8, coreprof_in(1)%profiles1d%q%value(1))
346  s_psi = sign(1.0_r8, coreprof_in(1)%psi%value(nrho)-coreprof_in(1)%psi%value(1))
347  s_jpar = sign(1.0_r8, coreprof_in(1)%profiles1d%jtot%value(1))
348 
349  WRITE(*,*) 's_ip, s_bt, s_q, s_psi, s_jpar = ', s_ip, s_bt, s_q, s_psi, s_jpar
350 
351 ! +++ Allocation of local variables:
352  ALLOCATE (rho(nrho))
353  ALLOCATE (rhomod(nrho))
354 
355  ALLOCATE (psi(nrho))
356  ALLOCATE (dpsi(nrho))
357  ALLOCATE (qsf(nrho))
358  ALLOCATE (jpar(nrho))
359 
360  ALLOCATE (fun(nrho))
361  ALLOCATE (dfun(nrho))
362  ALLOCATE (intfun(nrho))
363 
364  ALLOCATE (fdia(nrho))
365  ALLOCATE (gm2(nrho))
366  ALLOCATE (vprime(nrho))
367 
368  ALLOCATE (coreprof_out(1))
369 
370 
371 
372 ! +++ Local variables:
373  rho = coreprof_in(1)%rho_tor
374  psi = coreprof_in(1)%psi%value
375  qsf = coreprof_in(1)%profiles1d%q%value
376  jpar = coreprof_in(1)%profiles1d%jtot%value
377 
378  r0 = coreprof_in(1)%toroid_field%r0
379  b0 = coreprof_in(1)%toroid_field%b0
380 
381  CALL l3deriv(equilibrium_in(1)%profiles_1d%volume, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
382  vprime, rho, nrho)
383  CALL l3interp(equilibrium_in(1)%profiles_1d%F_dia, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
384  fdia, rho, nrho)
385  CALL l3interp(equilibrium_in(1)%profiles_1d%gm2, equilibrium_in(1)%profiles_1d%rho_tor, npsi, &
386  gm2, rho, nrho)
387 
388 
389 
390 ! +++ Recalculate PSI, Q and CURRENT profiles
391  IF (flag.EQ.1) THEN !profile flag = 1: PSI is primary quantity
392 ! Recalculated profile of QSF:
393  CALL derivn(nrho, rho, psi, dpsi)
394  qsf = 2.e0_r8 * itm_pi * b0 * rho / dpsi
395  qsf(1) = qsf(2)
396 
397 
398 
399 ! Recalculated profile of Jpar:
400  fun = dpsi * r0 * b0 * vprime * gm2 / fdia
401  CALL derivn(nrho, rho, fun, dfun)
402  jpar = - dfun / itm_mu0 / r0 / 2.e0_r8 / itm_pi / vprime / b0**2 * fdia**2
403  jpar(1) = jpar(2)
404 
405 
406  IF (qsf(1)*s_q.GT.minval(qsf*s_q).AND.q0_flag.GT.0.0d0) THEN
407  rho_min = rho(1)
408  qsf_min = qsf(1)
409  DO irho = 2,nrho
410  IF (qsf(irho)*s_q.LT.qsf_min*s_q) THEN
411  rho_min = rho(irho)
412  qsf_min = qsf(irho)
413  END IF
414  END DO
415  DO irho = 1,nrho
416  IF (rho(irho).LT.rho_min) THEN
417  qsf(irho) = qsf_min
418  END IF
419  END DO
420  fun = vprime * gm2 / fdia * rho / qsf
421  CALL derivn(nrho, rho, fun, dfun)
422  jpar = dfun / itm_mu0 / vprime * fdia**2
423  jpar(1) = jpar(2)
424  END IF
425 
426 
427 
428  ELSE IF(flag.EQ.2) THEN !profile flag = 2: Q is primary quantity
429 ! Recalculated profile of PSI:
430  fun = 2.e0_r8 * itm_pi * b0 / qsf
431  fun(1) = fun(2)
432  CALL integral(nrho, rho, fun, psi)
433 
434 ! Recalculated profile of Jpar:
435  fun = vprime * gm2 / fdia * rho / qsf
436  CALL derivn(nrho, rho, fun, dfun)
437  jpar = dfun / itm_mu0 / vprime * fdia**2
438  jpar(1) = jpar(2)
439 
440 
441 
442  ELSE IF(flag.EQ.3) THEN !profile flag = 3: current is primary quantity
443 ! Recalculated profile of QSF:
444  fun = jpar * itm_mu0 * vprime / fdia**2 / rho
445  fun(1) = fun(2)
446  CALL integral(nrho, rho, fun, intfun)
447  qsf = - vprime * gm2 / fdia * rho / intfun
448  qsf(1) = qsf(2)
449 
450 ! Recalculated profile of PSI:
451  fun = 2.e0_r8 * itm_pi * b0 / qsf
452  fun(1) = fun(2)
453  CALL integral(nrho, rho, fun, psi)
454 
455 
456  END IF
457 
458 
459 
460 
461 ! +++ Fill the iterration CPO with readjusted profiles:
462  CALL copy_cpo(coreprof_in(1), coreprof_out(1))
463 
464 
465 
466 ! +++ Update profiles
467  coreprof_out(1)%psi%value = psi
468  coreprof_out(1)%profiles1d%q%value = qsf
469  coreprof_out(1)%profiles1d%jtot%value = jpar
470 
471 
472 
473 ! +++ Deallocation of local variables:
474  DEALLOCATE (rho)
475  DEALLOCATE (rhomod)
476 
477  DEALLOCATE (psi)
478  DEALLOCATE (dpsi)
479  DEALLOCATE (qsf)
480  DEALLOCATE (jpar)
481 
482  DEALLOCATE (fun)
483  DEALLOCATE (dfun)
484  DEALLOCATE (intfun)
485 
486  DEALLOCATE (fdia)
487  DEALLOCATE (gm2)
488  DEALLOCATE (vprime)
489 
490 
491 
492  RETURN
493 
494  END SUBROUTINE readjust_profiles
495 
496 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
497 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
498 
499 
500 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
501 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
502 
503  SUBROUTINE correct_current_prof (J0_FLAG, COREPROF_IN, EQUILIBRIUM, COREPROF_OUT)
504 !-------------------------------------------------------!
505 ! This routine corrects profiles !
506 ! of psi, q and jparallel such, that !
507 ! there is no negative current. !
508 ! !
509 ! information received in: COREPROF_IN !
510 ! !
511 ! information saved in: COREPROF_OUT !
512 ! !
513 ! controling parameter: PROF_FLAG !
514 ! !
515 !-------------------------------------------------------!
516 ! Source: --- !
517 ! Developers: D.Kalupin !
518 ! Kontacts: Denis.Kalupin@efda.org !
519 ! !
520 ! Comments: created for V&V between ETS, JETTO !
521 ! and ASTRA !
522 ! !
523 !-------------------------------------------------------!
524 
525 
526  USE euitm_schemas
527  USE euitm_routines
528  USE itm_constants
529  USE copy_structures
530 
531  IMPLICIT NONE
532 
533  INTEGER :: ifail
534 
535 ! +++ CPO derived types:
536  TYPE (type_coreprof), POINTER :: coreprof_in(:) !input/output CPO with plasma profiles
537  TYPE (type_coreprof), POINTER :: coreprof_out(:) !output/output CPO with plasma profiles
538  TYPE (type_equilibrium), POINTER :: equilibrium(:) !input CPO with geometry quantities from previous time
539 
540 
541 
542 ! +++ Dimensions:
543  INTEGER :: nrho, irho !number and index of radial points
544  INTEGER :: npsi, ipsi !number and index of radial points
545  INTEGER :: j0_flag !flag for negative current, if J0_FLAG>0 negative current will be cut off
546 
547 
548 
549 ! +++ Local variables:
550  REAL (R8), ALLOCATABLE :: rho(:)
551  REAL (R8), ALLOCATABLE :: psi(:)
552  REAL (R8), ALLOCATABLE :: qsf(:)
553  REAL (R8), ALLOCATABLE :: jpar(:)
554  REAL (R8), ALLOCATABLE :: intjpar(:)
555 
556  REAL (R8), ALLOCATABLE :: fun(:)
557  REAL (R8), ALLOCATABLE :: intfun(:)
558 
559  REAL (R8), ALLOCATABLE :: vprime(:)
560  REAL (R8), ALLOCATABLE :: fdia(:)
561  REAL (R8), ALLOCATABLE :: gm2(:)
562 
563  REAL (R8) :: curr_total, curr
564  REAL (R8) :: r0, b0
565 
566  REAL (R8), ALLOCATABLE :: jpar_bg(:)
567  REAL (R8) :: curr_bg, bg
568 
569  REAL (R8) :: s_ip, s_bt, s_q, s_psi, s_jpar
570 
571 ! +++ Set dimensions and allocate parameters:
572  nrho = SIZE (coreprof_in(1)%rho_tor, dim=1)
573  npsi = SIZE (equilibrium(1)%profiles_1d%psi, dim=1)
574 
575  s_ip = sign(1.0_r8, coreprof_in(1)%globalparam%current_tot)
576  s_bt = sign(1.0_r8, coreprof_in(1)%toroid_field%b0)
577  s_q = sign(1.0_r8, coreprof_in(1)%profiles1d%q%value(1))
578  s_psi = sign(1.0_r8, coreprof_in(1)%psi%value(nrho)-coreprof_in(1)%psi%value(1))
579  s_jpar = sign(1.0_r8, coreprof_in(1)%profiles1d%jtot%value(1))
580 
581 ! WRITE(*,*) 's_ip, s_bt, s_q, s_psi, s_jpar = ', s_ip, s_bt, s_q, s_psi, s_jpar
582 
583 ! +++ Allocation of local variables:
584  ALLOCATE (rho(nrho))
585 
586  ALLOCATE (psi(nrho))
587  ALLOCATE (qsf(nrho))
588  ALLOCATE (jpar(nrho))
589  ALLOCATE (intjpar(nrho))
590 
591  ALLOCATE (fun(nrho))
592  ALLOCATE (intfun(nrho))
593 
594  ALLOCATE (vprime(nrho))
595  ALLOCATE (fdia(nrho))
596  ALLOCATE (gm2(nrho))
597 
598  ALLOCATE (jpar_bg(nrho))
599 
600  ALLOCATE (coreprof_out(1))
601 
602 
603 
604 ! +++ Local variables:
605  rho = coreprof_in(1)%rho_tor
606  psi = coreprof_in(1)%psi%value
607  qsf = coreprof_in(1)%profiles1d%q%value
608  jpar = coreprof_in(1)%profiles1d%jtot%value
609 
610  curr_total = coreprof_in(1)%globalparam%current_tot
611  r0 = coreprof_in(1)%toroid_field%r0
612  b0 = coreprof_in(1)%toroid_field%b0
613 
614  CALL l3deriv(equilibrium(1)%profiles_1d%volume, equilibrium(1)%profiles_1d%rho_tor, npsi, &
615  vprime, rho, nrho)
616  CALL l3interp(equilibrium(1)%profiles_1d%F_dia, equilibrium(1)%profiles_1d%rho_tor, npsi, &
617  fdia, rho, nrho)
618  CALL l3interp(equilibrium(1)%profiles_1d%gm2, equilibrium(1)%profiles_1d%rho_tor, npsi, &
619  gm2, rho, nrho)
620 
621 
622 ! +++ Determination of negative current area:
623  IF (j0_flag.GT.0) THEN
624  DO irho = 1,nrho
625  IF (jpar(irho)*s_jpar.LE.0.0e0_r8) jpar(irho) = 1.0e-6_r8*s_jpar
626  END DO
627  END IF
628 
629 
630 !+++ Renormalize current:
631  fun = vprime * jpar / 2.0e0_r8 / itm_pi * b0 / fdia**2
632  CALL integral_value(nrho, rho, fun, intjpar)
633  curr = intjpar(nrho) * fdia(nrho)
634 
635 
636 !+++ Difference in total current between required and real:
637  curr_bg = curr_total - curr
638 
639  IF (curr_bg*s_ip.LE.0.0_r8) THEN
640  DO irho = 1,nrho
641  IF (intjpar(irho) * fdia(irho)*s_ip.GE.curr_total*s_ip) jpar(irho) = 1.0e-6_r8*s_jpar
642  END DO
643 
644  ELSE IF (curr_bg*s_ip.GT.0.0_r8) THEN
645  bg = 1.0_r8*s_jpar
646  jpar_bg = bg * (1.0_r8 - 0.5_r8/(1.0_r8+((rho(nrho)-rho)/(rho(nrho)*0.05_r8))**3))
647  7 fun = vprime * max(jpar*s_jpar,jpar_bg*s_jpar)*s_jpar / 2.0e0_r8 / itm_pi * b0 / fdia**2
648  CALL integral_value(nrho, rho, fun, intjpar)
649  curr = intjpar(nrho) * fdia(nrho)
650 
651  IF (abs(1.0_r8-curr/curr_total).GE.1.e-4_r8) THEN
652  bg = bg * curr_total / curr
653  jpar_bg = jpar_bg * curr_total / curr
654  goto 7
655  END IF
656 
657  DO irho = 1,nrho
658  IF (jpar_bg(irho)*s_jpar.GT.jpar(irho)*s_jpar) THEN
659  jpar(irho:nrho) = jpar_bg(irho:nrho)
660  goto 10
661  END IF
662  END DO
663 ! Jpar = MAX(Jpar,Jpar_BG)
664 
665  END IF
666 
667 
668 
669 
670 
671 
672 
673 ! Recalculated profile of QSF:
674  10 fun = jpar * itm_mu0 * vprime / fdia**2 / rho
675  fun(1) = fun(2)
676  CALL integral(nrho, rho, fun, intfun)
677  qsf = - vprime * gm2 / fdia * rho / intfun
678  qsf(1) = qsf(2)
679 
680 
681 ! Recalculated profile of PSI:
682  fun = 2.e0_r8 * itm_pi * b0 / qsf
683  fun(1) = fun(2)
684  CALL integral(nrho, rho, fun, psi)
685 
686 
687 
688 
689 ! +++ Fill the iterration CPO with readjusted profiles:
690  CALL copy_cpo(coreprof_in(1), coreprof_out(1))
691 
692 
693 
694 ! +++ Update profiles
695  coreprof_out(1)%psi%value = psi
696  coreprof_out(1)%profiles1d%q%value = qsf
697  coreprof_out(1)%profiles1d%jtot%value = jpar
698 
699 
700 
701 ! +++ Deallocation of local variables:
702  DEALLOCATE (rho)
703 
704  DEALLOCATE (psi)
705  DEALLOCATE (qsf)
706  DEALLOCATE (jpar)
707  DEALLOCATE (intjpar)
708 
709  DEALLOCATE (fun)
710  DEALLOCATE (intfun)
711 
712  DEALLOCATE (vprime)
713  DEALLOCATE (fdia)
714  DEALLOCATE (gm2)
715 
716  DEALLOCATE (jpar_bg)
717 
718 
719 
720  RETURN
721 
722 
723  END SUBROUTINE correct_current_prof
724 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
725 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
726 
727 
728 
729 
730 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
731 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
732 
733  SUBROUTINE negative_current (J0_FLAG, COREPROF_IN, EQUILIBRIUM, COREPROF_OUT)
734 !-------------------------------------------------------!
735 ! This routine corrects profiles !
736 ! of psi, q and jparallel such, that !
737 ! there is no negative current. !
738 ! !
739 ! information received in: COREPROF_IN !
740 ! !
741 ! information saved in: COREPROF_OUT !
742 ! !
743 ! controling parameter: PROF_FLAG !
744 ! !
745 !-------------------------------------------------------!
746 ! Source: --- !
747 ! Developers: D.Kalupin !
748 ! Kontacts: Denis.Kalupin@efda.org !
749 ! !
750 ! Comments: created for V&V between ETS, JETTO !
751 ! and ASTRA !
752 ! !
753 !-------------------------------------------------------!
754 
755 
756  USE euitm_schemas
757  USE euitm_routines
758  USE itm_constants
759  USE copy_structures
760 
761  IMPLICIT NONE
762 
763  INTEGER :: ifail
764 
765 ! +++ CPO derived types:
766  TYPE (type_coreprof), POINTER :: coreprof_in(:) !input/output CPO with plasma profiles
767  TYPE (type_coreprof), POINTER :: coreprof_out(:) !output/output CPO with plasma profiles
768  TYPE (type_equilibrium), POINTER :: equilibrium(:) !input CPO with geometry quantities from previous time
769 
770 
771 
772 ! +++ Dimensions:
773  INTEGER :: nrho, irho !number and index of radial points
774  INTEGER :: npsi
775  INTEGER :: j0_flag !flag for negative current, if J0_FLAG>0 negative current will be cut off
776  INTEGER :: izero !inner position of negative current region
777 
778 
779 
780 ! +++ Local variables:
781  REAL (R8), ALLOCATABLE :: rho(:)
782  REAL (R8), ALLOCATABLE :: psi(:)
783  REAL (R8), ALLOCATABLE :: qsf(:)
784  REAL (R8), ALLOCATABLE :: jpar(:)
785  REAL (R8), ALLOCATABLE :: intjpar(:)
786 
787  REAL (R8), ALLOCATABLE :: fun(:)
788  REAL (R8), ALLOCATABLE :: intfun(:)
789 
790  REAL (R8), ALLOCATABLE :: vprime(:)
791  REAL (R8), ALLOCATABLE :: fdia(:)
792  REAL (R8), ALLOCATABLE :: gm2(:)
793 
794  REAL (R8) :: curr_total, curr
795  REAL (R8) :: r0, b0
796 
797  REAL (R8) :: s_ip, s_bt, s_q, s_psi, s_jpar
798 
799 
800 
801 ! +++ Set dimensions and allocate parameters:
802  nrho = SIZE (coreprof_in(1)%rho_tor, dim=1)
803  npsi = SIZE (equilibrium(1)%profiles_1d%psi, dim=1)
804 
805  s_ip = sign(1.0_r8, coreprof_in(1)%globalparam%current_tot)
806  s_bt = sign(1.0_r8, coreprof_in(1)%toroid_field%b0)
807  s_q = sign(1.0_r8, coreprof_in(1)%profiles1d%q%value(1))
808  s_psi = sign(1.0_r8, coreprof_in(1)%psi%value(nrho)-coreprof_in(1)%psi%value(1))
809  s_jpar = sign(1.0_r8, coreprof_in(1)%profiles1d%jtot%value(1))
810 
811  WRITE(*,*) 's_ip, s_bt, s_q, s_psi, s_jpar = ', s_ip, s_bt, s_q, s_psi, s_jpar
812 
813 ! +++ Allocation of local variables:
814  ALLOCATE (rho(nrho))
815 
816  ALLOCATE (psi(nrho))
817  ALLOCATE (qsf(nrho))
818  ALLOCATE (jpar(nrho))
819  ALLOCATE (intjpar(nrho))
820 
821  ALLOCATE (fun(nrho))
822  ALLOCATE (intfun(nrho))
823 
824  ALLOCATE (vprime(nrho))
825  ALLOCATE (fdia(nrho))
826  ALLOCATE (gm2(nrho))
827 
828 
829  ALLOCATE (coreprof_out(1))
830 
831 
832 
833 ! +++ Local variables:
834  rho = coreprof_in(1)%rho_tor
835  psi = coreprof_in(1)%psi%value
836  qsf = coreprof_in(1)%profiles1d%q%value
837  jpar = coreprof_in(1)%profiles1d%jtot%value
838  curr_total = coreprof_in(1)%globalparam%current_tot
839 
840  r0 = coreprof_in(1)%toroid_field%r0
841  b0 = coreprof_in(1)%toroid_field%b0
842 
843 
844  CALL l3deriv(equilibrium(1)%profiles_1d%volume, equilibrium(1)%profiles_1d%rho_tor, npsi, &
845  vprime, rho, nrho)
846  CALL l3interp(equilibrium(1)%profiles_1d%F_dia, equilibrium(1)%profiles_1d%rho_tor, npsi, &
847  fdia, rho, nrho)
848  CALL l3interp(equilibrium(1)%profiles_1d%gm2, equilibrium(1)%profiles_1d%rho_tor, npsi, &
849  gm2, rho, nrho)
850 
851 
852 ! +++ Determination of negative current area:
853  DO irho = 1,nrho
854  IF (jpar(irho)*s_jpar.LE.0.0e0_r8) jpar(irho) = 1.0e-6_r8*s_jpar
855  END DO
856 
857 
858 !+++ Renormalize current:
859  fun = vprime * jpar / 2.0e0_r8 / itm_pi * b0 / fdia**2
860  CALL integral_value(nrho, rho, fun, intjpar)
861  curr = intjpar(nrho) * fdia(nrho)
862 
863  DO irho = 1,nrho
864  IF (intjpar(irho) * fdia(irho)*s_jpar.GT.curr_total*s_jpar) jpar(irho) = 1.0e-6_r8*s_jpar
865  END DO
866 
867 
868 ! +++ Fill the iterration CPO with readjusted profiles:
869  CALL copy_cpo(coreprof_in(1), coreprof_out(1))
870 
871 
872 
873 ! +++ Update profiles
874  IF (j0_flag.GT.0) coreprof_out(1)%profiles1d%jtot%value = jpar
875 
876 
877 
878 
879 ! +++ Deallocation of local variables:
880  DEALLOCATE (rho)
881 
882  DEALLOCATE (psi)
883  DEALLOCATE (qsf)
884  DEALLOCATE (jpar)
885  DEALLOCATE (intjpar)
886 
887  DEALLOCATE (fun)
888  DEALLOCATE (intfun)
889 
890  DEALLOCATE (vprime)
891  DEALLOCATE (fdia)
892  DEALLOCATE (gm2)
893 
894 
895 
896  RETURN
897 
898 
899  END SUBROUTINE negative_current
900 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
901 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
902 
903 
904 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
913 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
914 SUBROUTINE integral(N,X,Y,INTY)
915 !-------------------------------------------------------!
916 ! This subroutine calculates integral of function !
917 ! Y(X)*X from X=0 until X=X(N) !
918 !-------------------------------------------------------!
919 
920  USE itm_types
921 
922  IMPLICIT NONE
923 
924  INTEGER :: n ! number of radial points (input)
925  INTEGER :: i
926 
927  REAL (R8) :: x(n), & ! argument array (input)
928  y(n), & ! function array (input)
929  inty(n) ! function integral array (output)
930 
931  inty(1)=y(1)*x(1)**2/2.e0_r8
932  DO i=2,n
933  inty(i)=inty(i-1)+(y(i-1)*x(i-1)+y(i)*x(i))*(x(i)-x(i-1))/2.e0_r8
934  END DO
935 
936  RETURN
937 END SUBROUTINE integral
938 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
939 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
940 
941 
942 
943 
944 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
945 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
946 ! This subroutine calculates integral of a function y(x)
947 ! from X=0
948 SUBROUTINE integral_value(N,X,Y,INTY)
949  USE itm_types
950  IMPLICIT NONE
951 
952  INTEGER :: n ! number of radial points (input)
953  INTEGER :: i
954 
955  REAL (R8) :: x(n), & ! argument array (input)
956  y(n), & ! function array (input)
957  inty(n) ! function integral array (output)
958 
959 
960  inty(1)=y(1)*x(1)/2.0_r8
961  DO i=2,n
962  inty(i)=inty(i-1)+(y(i)+y(i-1))*(x(i)-x(i-1))/2.0_r8
963  ENDDO
964  RETURN
965  END SUBROUTINE integral_value
966  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
967  ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
968 
969 
970 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
979 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
980 SUBROUTINE derivn(N,X,Y,DY1)
981 
982 !-------------------------------------------------------!
983 ! These subroutines calculate first and second !
984 ! derivatives, DY1 and DY2, of function Y respect !
985 ! to argument X !
986 !-------------------------------------------------------!
987 
988  USE itm_types
989  IMPLICIT NONE
990 
991  INTEGER :: n ! number of radial points (input)
992  INTEGER :: i
993 
994  REAL (R8) :: x(n), & ! argument array (input)
995  y(n), & ! function array (input)
996  dy1(n) ! function derivative array (output)
997  REAL (R8) :: h(n),dy2(n)
998 
999  DO i=1,n-1
1000  h(i)=x(i+1)-x(i)
1001  END DO
1002 
1003  DO i=2,n-1
1004  dy1(i)=((y(i+1)-y(i))*h(i-1)/h(i)+(y(i)-y(i-1))*h(i)/h(i-1)) &
1005  /(h(i)+h(i-1))
1006  dy2(i)=2.e0_r8*((y(i-1)-y(i))/h(i-1)+(y(i+1)-y(i))/h(i)) &
1007  /(h(i)+h(i-1))
1008  END DO
1009 
1010  dy1(1)=dy1(2)-dy2(2)*h(1)
1011  dy1(n)=dy1(n-1)+dy2(n-1)*h(n-1)
1012 
1013  RETURN
1014 END SUBROUTINE derivn
1015 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1016 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1017 
1018 
1019 END MODULE equilibrium_work
l3deriv
subroutine l3deriv(y_in, x_in, nr_in, dydx_out, x_out, nr_out)
Definition: l3interp.f90:59
fun
subroutine fun(X, F)
Definition: Ev2.f:10
derivn
subroutine derivn(N, X, Y, DY1)
These subroutines calculate first and second derivatives, DY1 and DY2, of function Y respect to argum...
Definition: main_plasma.f90:3856
equilibrium_work
Definition: equilibrium_work.f90:1
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
equilibrium_work::parabolic_prof
subroutine parabolic_prof(COREPROF_IN, EQUILIBRIUM_IN, COREPROF_OUT, EQUILIBRIUM_OUT)
Definition: equilibrium_work.f90:9
helena21_mod::fdia
real(r8) function fdia(psi_n)
Definition: helena21_mod.f90:139
equilibrium_work::negative_current
subroutine negative_current(J0_FLAG, COREPROF_IN, EQUILIBRIUM, COREPROF_OUT)
Definition: equilibrium_work.f90:733
equilibrium_work::readjust_profiles
subroutine readjust_profiles(PROF_FLAG, Q0_FLAG, COREPROF_IN, EQUILIBRIUM_IN, COREPROF_OUT)
Definition: equilibrium_work.f90:261
integral_value
subroutine integral_value(N, X, Y, INTY)
Definition: prepare_input_cpos.F90:1897
integral
subroutine integral(n, h, r, f, int)
Definition: solution2.f90:527
equilibrium_work::correct_current_prof
subroutine correct_current_prof(J0_FLAG, COREPROF_IN, EQUILIBRIUM, COREPROF_OUT)
Definition: equilibrium_work.f90:503