ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
phys_functions.f90
Go to the documentation of this file.
1 module phys_functions
2 
3  use itm_types
4  use mod_parameter
5 
6  implicit none
7 
8  integer(itm_i4), parameter, private :: iu6 = 6
9 
10  contains
11 
12 !---------------------------------------------------------------------------
13  function identity(a, x, xr, xs, yr, ys, psi, psir, F_dia) result(f_identity)
14 !---------------------------------------------------------------------------
15 ! simply returns 1
16 !---------------------------------------------------------------------------
17 
18  implicit none
19 
20  real(r8), intent(in) :: a
21  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
22  real(r8), intent(in) :: f_dia
23  real(r8) :: f_identity
24 
25  f_identity = 1._r8
26 
27  end function identity
28 
29 !---------------------------------------------------------------------------
30  function rh_side(a, x, xr, xs, yr, ys, psi, psir, F_dia) result(f_rh_side)
31 !---------------------------------------------------------------------------
32 ! calculates the right hand side of the Grad Shafranov equation rh_side
33 ! depending on hbt and isol
34 !---------------------------------------------------------------------------
35 
36  use mod_dat
37 
38  implicit none
39 
40  real(r8), intent(in) :: a
41  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
42  real(r8), intent(in) :: f_dia
43  real(r8) :: f_rh_side
44  real(r8) :: dgamma_dpsi, dp_dpsi
45 
46  if (hbt) then
47  f_rh_side = dgamma_dpsi(psi) + b * x * (1._r8 + eps * x / 2._r8) &
48  * dp_dpsi(psi)
49  else
50  f_rh_side = dgamma_dpsi(psi) + b * (1._r8 + eps * x)**2 &
51  * dp_dpsi(psi)
52  end if
53  f_rh_side = a * f_rh_side / (1._r8 + eps * x)
54 !-- Soloviev rhs
55  if (isol == 1) f_rh_side = a * (1._r8 + b * x * (1._r8 + eps * x &
56  / 2._r8)) / (1._r8 + eps * x)
57 
58  end function rh_side
59 
60 !---------------------------------------------------------------------------
61  function one_over_r(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
62  result(f_one_over_r)
63 !---------------------------------------------------------------------------
64 ! returns 1 / (1 + eps * x)
65 !---------------------------------------------------------------------------
66 
67  use mod_dat
68 
69  implicit none
70 
71  real(r8), intent(in) :: a
72  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
73  real(r8), intent(in) :: f_dia
74  real(r8) :: f_one_over_r
75 
76  f_one_over_r = 1._r8 / (1._r8 + eps * x)
77 
78  end function one_over_r
79 
80 !---------------------------------------------------------------------------
81  function rquad_minus_1_over_2_r(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
82  result(f_rquad_minus_1_over_2_r)
83 !---------------------------------------------------------------------------
84 ! returns x * (1 + eps * x / 2) / (1 + eps * x)
85 !---------------------------------------------------------------------------
86 
87  use mod_dat
88 
89  implicit none
90 
91  real(r8), intent(in) :: a
92  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
93  real(r8), intent(in) :: f_dia
94  real(r8) :: f_rquad_minus_1_over_2_r
95 
96  f_rquad_minus_1_over_2_r = x * (1._r8 + eps * x / 2._r8) &
97  / (1._r8 + eps * x)
98 
99  end function rquad_minus_1_over_2_r
100 
101 !---------------------------------------------------------------------------
102  function r_major(a, x, xr, xs, yr, ys, psi, psir, F_dia) result(f_R_major)
103 !---------------------------------------------------------------------------
104 ! returns 1 + eps * x
105 !---------------------------------------------------------------------------
106 
107  use mod_dat
108 
109  implicit none
110 
111  real(r8), intent(in) :: a
112  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
113  real(r8), intent(in) :: f_dia
114  real(r8) :: f_r_major
115 
116  f_r_major = 1._r8 + eps * x
117 
118  end function r_major
119 
120 !---------------------------------------------------------------------------
121  function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia) result(f_R2)
122 !---------------------------------------------------------------------------
123 ! returns (1 + eps * x)**2
124 !---------------------------------------------------------------------------
125 
126  use mod_dat
127 
128  implicit none
129 
130  real(r8), intent(in) :: a
131  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
132  real(r8), intent(in) :: f_dia
133  real(r8) :: f_r2
134 
135  f_r2 = r_major(a, x, xr, xs, yr, ys, psi, psir, f_dia)**2
136 
137  end function r2
138 
139 !---------------------------------------------------------------------------
140  function one_over_r2(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
141  result(f_one_over_r2)
142 !---------------------------------------------------------------------------
143 ! returns 1 / (1 + eps * x)**2
144 !---------------------------------------------------------------------------
145 
146  use mod_dat
147 
148  implicit none
149 
150  real(r8), intent(in) :: a
151  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
152  real(r8), intent(in) :: f_dia
153  real(r8) :: f_one_over_r2
154 
155  f_one_over_r2 = 1._r8 / r2(a, x, xr, xs, yr, ys, psi, psir, f_dia)
156 
157  end function one_over_r2
158 
159 !---------------------------------------------------------------------------
160  function gradpsi2(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
161  result(f_gradpsi2)
162 !---------------------------------------------------------------------------
163 ! returns g11_hel / (2*pi)**2 normalized to (cpsurfin / a_minor)**2
164 ! with g11_hel = grad(psi) and psi the total poloidal flux
165 !---------------------------------------------------------------------------
166 
167  use mod_dat
168 
169  implicit none
170 
171  real(r8), intent(in) :: a
172  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
173  real(r8), intent(in) :: f_dia
174  real(r8) :: jacobian
175 
176  real(r8) :: f_gradpsi2
177 
178  jacobian = xr * ys - xs * yr
179 
180  if (jacobian /= 0._r8) then
181  f_gradpsi2 = psir**2 * (xs**2 + ys**2) / jacobian**2
182  else
183  if (psir /= 0._r8 .and. (xs /= 0._r8 .or. ys /= 0._r8)) then
184  write(iu6, *) 'ERROR: Division by zero in function gradpsi2! '
185  stop
186  else
187  f_gradpsi2 = 0._r8
188  end if
189  end if
190 
191  f_gradpsi2 = f_gradpsi2 / twopi**2
192 
193  end function gradpsi2
194 
195 !---------------------------------------------------------------------------
196  function gradpsi(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
197  result(f_gradpsi)
198 !---------------------------------------------------------------------------
199 ! returns sqrt(g11_hel) / (2 * pi) normalized to (cpsurfin / a_minor)
200 ! with g11_hel = grad(psi) and psi the total poloidal flux
201 !---------------------------------------------------------------------------
202 
203  use mod_dat
204 
205  implicit none
206 
207  real(r8), intent(in) :: a
208  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
209  real(r8), intent(in) :: f_dia
210 
211  real(r8) :: f_gradpsi
212 
213  f_gradpsi = sqrt(gradpsi2(a, x, xr, xs, yr, ys, psi, psir, f_dia))
214 
215  end function gradpsi
216 
217 !---------------------------------------------------------------------------
218  function gradpsi2_over_r2(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
219  result(f_gradpsi2_over_r2)
220 !---------------------------------------------------------------------------
221 ! returns g11_hel / (2* pi * R)**2 normalized to
222 ! (cpsurfin / (a_minor * Rvac))**2
223 ! with g11_hel = grad(psi) and psi the total poloidal flux
224 !---------------------------------------------------------------------------
225 
226  use mod_dat
227 
228  implicit none
229 
230  real(r8), intent(in) :: a
231  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
232  real(r8), intent(in) :: f_dia
233 
234  real(r8) :: f_gradpsi2_over_r2
235 
236  f_gradpsi2_over_r2 = gradpsi2(a, x, xr, xs, yr, ys, psi, psir, f_dia) &
237  / r2(a, x, xr, xs, yr, ys, psi, psir, f_dia)
238 
239  end function gradpsi2_over_r2
240 
241 !---------------------------------------------------------------------------
242  function f2_plus_gradpsi2(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
243  result(f_f2_plus_gradpsi2)
244 !---------------------------------------------------------------------------
245 ! returns F**2 + g11_hel / (2 * pi)**2 normalized to (cpsurfin / a_minor)**2
246 ! with g11_hel = grad(psi) and psi the total poloidal flux
247 !---------------------------------------------------------------------------
248 
249  use mod_dat
250 
251  implicit none
252 
253  real(r8), intent(in) :: a
254  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
255  real(r8), intent(in) :: f_dia
256 
257  real(r8) :: f_f2_plus_gradpsi2
258 
259  f_f2_plus_gradpsi2 = f_dia**2 + gradpsi2(a, x, xr, xs, yr, ys, psi, &
260  psir, f_dia)
261 
262  end function f2_plus_gradpsi2
263 
264 !---------------------------------------------------------------------------
265  function b_fun(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
266  result(f_b2)
267 !---------------------------------------------------------------------------
268 ! returns sqrt(F**2 + g11_hel / (2 * pi)**2) / R**2 normalized to
269 ! (cpsurfin / (a_minor * Rvac))
270 ! with g11_hel = grad(psi) and psi the total poloidal flux
271 !---------------------------------------------------------------------------
272 
273  use mod_dat
274 
275  implicit none
276 
277  real(r8), intent(in) :: a
278  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
279  real(r8), intent(in) :: f_dia
280 
281  real(r8) :: f_b2
282 
283  f_b2 = sqrt((f_dia**2 + gradpsi2(a, x, xr, xs, yr, ys, psi, psir, f_dia)) &
284  / r2(a, x, xr, xs, yr, ys, psi, psir, f_dia))
285 
286  end function b_fun
287 
288 !---------------------------------------------------------------------------
289  function b2(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
290  result(f_b2)
291 !---------------------------------------------------------------------------
292 ! returns (F**2 + g11_hel / (2 * pi)**2) / R**2 normalized to
293 ! (cpsurfin / (a_minor * Rvac))**2
294 ! with g11_hel = grad(psi) and psi the total poloidal flux
295 !---------------------------------------------------------------------------
296 
297  use mod_dat
298 
299  implicit none
300 
301  real(r8), intent(in) :: a
302  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
303  real(r8), intent(in) :: f_dia
304 
305  real(r8) :: f_b2
306 
307  f_b2 = (f_dia**2 + gradpsi2(a, x, xr, xs, yr, ys, psi, psir, f_dia)) &
308  / r2(a, x, xr, xs, yr, ys, psi, psir, f_dia)
309 
310  end function b2
311 
312 !---------------------------------------------------------------------------
313  function one_over_b2(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
314  result(f_one_over_b2)
315 !---------------------------------------------------------------------------
316 ! returns R**2 / (F**2 + g11_hel / (2 * pi)**2) normalized to
317 ! (a_minor * Rvac / cpsurfin)**2
318 ! with g11_hel = grad(psi) and psi the total poloidal flux
319 !---------------------------------------------------------------------------
320 
321  use mod_dat
322 
323  implicit none
324 
325  real(r8), intent(in) :: a
326  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
327  real(r8), intent(in) :: f_dia
328 
329  real(r8) :: f_one_over_b2
330 
331  f_one_over_b2 = 1._r8 / b2(a, x, xr, xs, yr, ys, psi, psir, f_dia)
332 
333  end function one_over_b2
334 
335 !---------------------------------------------------------------------------
336  function gradrho2_over_b2(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
337  result(f_gradrho2_over_b2)
338 !---------------------------------------------------------------------------
339 ! returns R**2 * g11_hel / (2 *pi)**2 / (F**2 + g11_hel / (2 * pi)**2)
340 ! normalized to Rvac**2
341 ! with g11_hel = grad(psi) and psi the total poloidal flux
342 !---------------------------------------------------------------------------
343 
344  use mod_dat
345 
346  implicit none
347 
348  real(r8), intent(in) :: a
349  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
350  real(r8), intent(in) :: f_dia
351 
352  real(r8) :: f_gradrho2_over_b2
353 
354  f_gradrho2_over_b2 = gradpsi2(a, x, xr, xs, yr, ys, psi, psir, f_dia) &
355  / b2(a, x, xr, xs, yr, ys, psi, psir, f_dia)
356 
357  end function gradrho2_over_b2
358 
359 !---------------------------------------------------------------------------
360  function bp2(a, x, xr, xs, yr, ys, psi, psir, F_dia) &
361  result(f_bp2)
362 !---------------------------------------------------------------------------
363 ! returns g11_hel / (2 * pi * R)**2 normalized to
364 ! (cpsurfin / (a_minor * Rvac))**2
365 ! with g11_hel = grad(psi) and psi the total poloidal flux
366 !---------------------------------------------------------------------------
367 
368  use mod_dat
369 
370  implicit none
371 
372  real(r8), intent(in) :: a
373  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
374  real(r8), intent(in) :: f_dia
375 
376  real(r8) :: f_bp2
377 
378  f_bp2 = gradpsi2_over_r2(a, x, xr, xs, yr, ys, psi, psir, f_dia)
379 
380  end function bp2
381 
382 !---------------------------------------------------------------------------
383  function p(a, x, xr, xs, yr, ys, psi, psir, F_dia) result(f_p)
384 !---------------------------------------------------------------------------
385 ! returns p (in SI)
386 !---------------------------------------------------------------------------
387 
388  use mod_dat, only : b, alfa, bvac, p_sep, hbt, eps, pscale, bmag
389  use mod_parameter, only : mu0
390 
391  implicit none
392 
393  real(r8), intent(in) :: a
394  real(r8), intent(in) :: x, xr, xs, yr, ys, psi, psir
395  real(r8), intent(in) :: f_dia
396  real(r8) :: f_p
397 
398  if (hbt) then
399  f_p = 0.5_r8 * a * b * pressure(psi) * sign(1.0_r8, -alfa) &
400  * sign(1.0_r8, bvac)
401  else
402  f_p = eps * a * b * pressure(psi) * sign(1.0_r8, -alfa) &
403  * sign(1.0_r8, bvac)
404  end if
405 !-- denormalize
406  f_p = f_p * pscale
407  f_p = f_p * bmag**2 / mu0
408 !-- add separatrix pressure
409  f_p = f_p + p_sep
410 
411  end function p
412 
413 !---------------------------------------------------------------------------
414  function pressure(flux) result(f_pressure)
415 !------------------------------------------------------------------
416 ! the normalized profile of the pressure versus flux
417 ! this routine must be initialized by a call to pressure_profile
418 !------------------------------------------------------------------
419 
420  use itm_types
421  use mod_profiles, only : psivec, npts
422  use mod_spline, only : p_spline
423  use helena_spline
424 
425  implicit none
426 
427  real(r8), intent(in) :: flux
428  real(r8) :: f_pressure
429 
430  real(r8), dimension(3) :: ablt
431  real(r8) :: flux_tmp
432 
433  flux_tmp = flux
434 
435  if (flux > 1._r8) then
436  flux_tmp = 1._r8
437  end if
438 
439  f_pressure = spwert(npts, flux_tmp, p_spline, psivec, ablt, 0)
440 
441 end function pressure
442 
443 end module phys_functions
mod_dat
Definition: mod_dat.f90:1
mod_spline
Definition: mod_spline.f90:1
phys_functions::bp2
real(r8) function bp2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:360
phys_functions::b_fun
real(r8) function b_fun(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:265
phys_functions::gradpsi2
real(r8) function gradpsi2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:160
phys_functions::one_over_r2
real(r8) function one_over_r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:140
phys_functions::gradpsi
real(r8) function gradpsi(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:196
flux
subroutine flux(psitok, rk, zk, nk)
Definition: EQ1_m.f:786
spwert
REAL *8 function spwert(N, XWERT, A, B, C, D, X, ABLTG)
Definition: solution6.f90:155
phys_functions::rh_side
real(r8) function rh_side(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:30
phys_functions::f2_plus_gradpsi2
real(r8) function f2_plus_gradpsi2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:242
phys_functions::one_over_b2
real(r8) function one_over_b2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:313
mod_profiles
Definition: mod_profiles.f90:1
dgamma_dpsi
real(r8) function dgamma_dpsi(flux)
Definition: dgamma_dpsi.f90:1
phys_functions::b2
real(r8) function b2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:289
mod_parameter
Definition: mod_parameter.f90:1
phys_functions::gradrho2_over_b2
real(r8) function gradrho2_over_b2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:336
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
phys_functions::identity
real(r8) function identity(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:13
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
phys_functions::rquad_minus_1_over_2_r
real(r8) function rquad_minus_1_over_2_r(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:81
phys_functions::one_over_r
real(r8) function one_over_r(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:61
phys_functions::r_major
real(r8) function r_major(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:102
dp_dpsi
real(r8) function dp_dpsi(flux)
Definition: dp_dpsi.f90:1
phys_functions
Definition: phys_functions.f90:1
phys_functions::pressure
real(r8) function pressure(flux)
Definition: phys_functions.f90:414
phys_functions::gradpsi2_over_r2
real(r8) function gradpsi2_over_r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:218
helena_spline
Definition: helena_spline.f90:1