ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
mod_spider.f90
Go to the documentation of this file.
1 module mod_spider
2 
3 contains
4 
5 ! program main_astr_fixed
6 ! Create the library ~/MAstra/.lbr/Intel/esp.a as (see makefile)
7 ! $> cd ~/MAstra/esp/ ; make ALL
8 
9  SUBROUTINE int_spider(&
10 ! Control:
11  key_equil &
12 ! Input:
13  ,neql & ! radial grid in SPIDER (150)
14  ,nteta & ! poloidal grid in SPIDER (100)
15  ,nbnd & ! # boundary points (101)
16  ,rzbnd & ! array(2*nbnd) of boundary points
17  ,na1 & ! radial grid for input (dimension of all arrays below)
18  ,eqpf & !
19  ,eqff & ! A
20  ,rho & ! S
21  ,ipl & ! T
22  ,rtor & ! R
23  ,ztor & ! Z
24  ,btor & ! A
25  ,roc & !
26 ! Output:
27  ,rocnew & ! n
28 ! Input:
29  ,nstep & ! o
30 ! Input (and output if key_dmf = -2, -3)
31  ,fp & ! t
32 ! Output:
33  ,g11 & ! a
34  ,g22 & ! t
35  ,g33 & ! i
36  ,vr & ! o
37  ,vrs & ! n
38  ,slat & ! s
39  ,gradro & !
40  ,mu & ! &
41  ,ipol & !
42  ,bmaxt &
43  ,bmint &
44  ,bdb02 &
45  ,bdb0 &
46  ,b0db2 &
47  ,droda &
48  ,rout &
49  ,zout &
50  ,shafr_shift &!
51  ,elong &!
52  ,elong_u &!
53  ,elong_l &!
54  ,triang_l &!
55  ,triang_u &!
56  ,amid &!
57  ,ftrap &!
58 ! Input:
59  ,cc & ! u
60  ,te & ! n
61  ,cubs & ! i
62  ,cd & ! t
63 ! Input:
64  ,pres & ! s
65  ,cu & !
66  )
67 
68  USE itm_constants
69  USE parameters
70 
71  IMPLICIT NONE
72 
73 
74 ! +++ I/O:
75  INTEGER :: na1
76  INTEGER :: nbnd, nstep, key_equil
77  INTEGER :: neql, nteta
78 
79  REAL (R8) :: cc(na1), cubs(na1), cd(na1),cu(na1)
80  REAL (R8) :: te(na1), pres(na1)
81  REAL (R8) :: cu_out(na1)
82  REAL (R8) :: pres_s(na1)
83  REAL (R8) :: ymu(na1),yfp(na1)
84  REAL (R8) :: rzbnd(2*nbnd),eqpf(na1),eqff(na1),fp(na1),ipl,ybetpl
85  REAL (R8) :: yli3,rtor,ztor,btor,rho(na1),roc, ftrap(na1)
86  REAL (R8) :: g11(na1),g22(na1),g33(na1)
87  REAL (R8) :: vr(na1),vrs(na1),slat(na1),gradro(na1),rocnew
88  REAL (R8) :: mu(na1),ipol(na1),bmaxt(na1),bmint(na1),bdb02(na1)
89  REAL (R8) :: bdb0(na1),b0db2(na1),droda(na1)
90  REAL (R8) :: shafr_shift(na1),elong(na1),elong_u(na1),elong_l(na1),triang_l(na1),triang_u(na1),amid(na1)
91  REAL (R8) :: roc_input
92  REAL (R8) :: rout(neql,nteta), zout(neql,nteta)
93 
94 
95 ! +++ locals:
96  REAL (R8) :: hro, yreler, platok,time, dt
97 
98 
99 ! +++ SPIDER I/O:
100  TYPE(type_parameters) :: parameters_spider
101 
102 
103 
104 
105  CALL input2spider2(neql,nteta,nbnd,rzbnd,na1,eqpf,eqff,fp,ipl, &
106  rtor,btor,rho,pres,cu, &
107  nstep, key_equil, &
108  parameters_spider)
109 
110 
111  CALL spider_run2(parameters_spider)
112 
113 
114  CALL spider2output2(parameters_spider, &
115  neql,nteta,rtor,ztor,btor,rho,roc,na1, &
116  g11,g22,g33,vr,vrs,slat,gradro,rocnew, &
117  ymu,ipol,bmaxt,bmint,bdb02,b0db2,bdb0,droda, &
118  yreler,yli3,platok,cu_out,yfp, &
119  pres_s,rout,zout,shafr_shift,elong,elong_u,elong_l, &
120  triang_l,triang_u,amid,ftrap)
121 
122 
123  RETURN
124 
125  END SUBROUTINE int_spider
126 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
127 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
128 
129 
130 
131 
132 
133 
134 
135 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
136 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
137  SUBROUTINE input2spider2(neql,nteta,nbnd,rzbnd,na1,eqpf,eqff,fp,ipl, &
138  rtor,btor,rho,pres,cu, &
139  nstep, key_equil, &
140  parameters_spider)
141 
142  USE itm_constants
143  USE parameters
144 
145  IMPLICIT NONE
146 
147 
148 ! +++ I/O:
149  INTEGER :: nbnd, key_equil, na1, nstep
150  INTEGER :: neql, nteta
151  REAL (R8) :: dt,time,dpsdt,yreler
152  REAL (R8) :: ipl, rtor, btor,roc
153  REAL (R8) :: rzbnd(1:2*nbnd)
154  REAL (R8) :: eqpf(na1), eqff(na1), fp(na1)
155  REAL (R8) :: rho(na1), pres(na1), cu(na1)
156 
157 
158 ! +++ locals:
159  INTEGER :: kpr,k_con,kname
160  INTEGER :: k_fixfree,key_ini,key_start,key_0stp,key_pres,key_dmf
161  INTEGER :: k_grid,k_auto
162  REAL (R8), ALLOCATABLE :: mu(:), cc(:),te(:),cubc(:),cd(:)
163 
164  CHARACTER*40 :: prename
165  CHARACTER*40 :: eqdfn
166 
167  TYPE(type_parameters) :: parameters_spider
168 
169 
170 
171  ALLOCATE ( mu(na1), cc(na1), te(na1), cubc(na1), cd(na1) )
172 
173  mu = 0.e0_r8
174  cc = 0.e0_r8
175  te = 0.e0_r8
176  cubc = 0.e0_r8
177  cd = 0.e0_r8
178 
179 
180 
181 ! +++ setup parameters for equilibrium problem to solve
182  parameters_spider%nstep = nstep
183 
184 
185  IF (nstep .NE. 0) THEN
186  parameters_spider%key_dmf = key_equil
187  ELSE
188  IF(key_equil .eq. 0) THEN
189  parameters_spider%key_0stp = key_equil
190  ELSE IF (key_equil .eq. -10) THEN
191  parameters_spider%key_0stp = 1
192  ELSE
193  parameters_spider%key_0stp = key_equil
194  END IF
195  END IF
196 
197 
198 
199 ! +++ p' input for key_equil=0 (otherwise the default key_pres is used)
200  IF (key_equil .eq. 0) THEN
201  parameters_spider%key_pres = 0
202  END IF
203 
204  kpr = parameters_spider%kpr !print in spider
205  k_con = 0
206  prename = parameters_spider%prename
207  kname = parameters_spider%kname
208 
209 
210  CALL kpr_calc(kpr)
211  CALL put_key_con(k_con)
212  CALL put_name(prename,kname)
213 
214  nstep = parameters_spider%nstep
215  time = parameters_spider%time
216  dt = parameters_spider%dt
217  key_dmf = parameters_spider%key_dmf
218  k_grid = parameters_spider%k_grid ! k_grid= 0 rect. grid; k_grid= 1 adap. grid
219  k_auto = parameters_spider%k_auto ! k_auto= 1-> full initialization; k_auto= 0-> preinitialization is assumed to be done
220  k_fixfree = parameters_spider%k_fixfree !=0->only fixed boundary spider
221  dpsdt = parameters_spider%dpsdt
222  key_ini = parameters_spider%key_ini ! =1 astra profiles, =0 start from EQDSK and SPIDER profiles
223  key_start = parameters_spider%key_start
224  eqdfn = parameters_spider%eqdfn
225  key_0stp = parameters_spider%key_0stp
226  key_pres = parameters_spider%key_pres
227 ! neql = PARAMETERS_SPIDER%neql
228 ! nteta = PARAMETERS_SPIDER%nteta
229 
230 
231 ! PARAMETERS_SPIDER%nstep = nstep
232  parameters_spider%neql = neql
233  parameters_spider%nteta = nteta
234 
235  CALL aspid_flag(1)
236 ! CALL put_key_fix(k_fixfree)
237 
238  roc = rho(na1)
239 
240 
241 
242  CALL astra2spider(neql,nteta,nbnd,rzbnd,key_dmf, &
243  na1,eqpf,eqff,fp,ipl, &
244  rtor,btor,rho,roc,nstep,yreler,mu, &
245  cc,te,cubc,cd,key_ini,eqdfn,pres,cu, &
246  key_0stp,key_pres )
247 
248 
249 
250 
251  DEALLOCATE ( mu, cc, te, cubc, cd )
252 
253 
254 
255  RETURN
256 
257  END SUBROUTINE input2spider2
258 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
259 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
260 
261 
262 
263 
264 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
265 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
266  SUBROUTINE spider_run2(PARAMETERS_SPIDER)
267 
268  USE itm_constants
269  USE parameters
270 
271  IMPLICIT NONE
272 
273  !locals
274  INTEGER :: nstep, key_dmf, k_grid, k_auto, k_fixfree, key_ini
275  INTEGER :: kpr, kname
276  REAL (R8) :: time, dt, dpsdt
277  INTEGER :: i,neql
278 
279  REAL (R8), ALLOCATABLE :: voltpf(:), phi(:)
280 
281  TYPE(type_parameters) :: parameters_spider
282 
283 
284 
285  ALLOCATE ( voltpf(100) )
286 
287  voltpf = 0.e0_r8
288 
289  nstep = parameters_spider%nstep
290  time = parameters_spider%time
291  dt = parameters_spider%dt
292  key_dmf = parameters_spider%key_dmf
293  k_grid = parameters_spider%k_grid
294  k_auto = parameters_spider%k_auto
295  k_fixfree = parameters_spider%k_fixfree
296  dpsdt = parameters_spider%dpsdt
297  key_ini = parameters_spider%key_ini
298  kpr = parameters_spider%kpr
299  kname = parameters_spider%kname
300 
301 
302 
303  CALL kpr_calc(kpr)
304  CALL put_name(parameters_spider%prename,kname)
305 
306  CALL put_tim(dt,time)
307  CALL savepsi
308 
309  CALL spider(nstep,time,dt,key_dmf,k_grid,k_auto,k_fixfree, &
310  dpsdt,key_ini,voltpf)
311 
312 ! CALL PLA_VOLT(dt)
313  CALL cur_avg
314  CALL wrb
315  CALL wr_spik
316 
317 
318  DEALLOCATE ( voltpf )
319 
320 
321  RETURN
322 
323  END SUBROUTINE spider_run2
324 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
325 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
326 
327 
328 
329 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
330 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
331  SUBROUTINE spider2output2(PARAMETERS_SPIDER, &
332  neql,nteta,rtor,ztor,btor,rho,roc,na1, &
333  g11,g22,g33,vr,vrs,slat,gradro,rocnew, &
334  mu,ipol,bmaxt,bmint,bdb02,b0db2,bdb0,droda, &
335  yreler,yli3,platok,cu_out,fp,pres_s, &
336  rout,zout,shafr_shift,elong,elong_u,elong_l,&
337  triang_l, triang_u,amid,ftrap)
338 
339  USE itm_constants
340  USE parameters
341 
342  IMPLICIT NONE
343 
344 
345  INTEGER :: na1,ni_p,nj_p
346  INTEGER :: i, j, neql, nteta, ndim2
347  INTEGER :: k1, k2
348 
349  REAL (R8) :: rtor,ztor,btor, platok,yli3,roc,rocnew,yreler
350  REAL (R8) :: g11(na1), g22(na1), g33(na1)
351  REAL (R8) :: rho(na1), pres_s(na1), cu_out(na1),fp(na1),ftrap(na1)
352  REAL (R8) :: vr(na1),vrs(na1),slat(na1)
353  REAL (R8) :: gradro(na1),droda(na1),mu(na1)
354  REAL (R8) :: ipol(na1),bmaxt(na1),bmint(na1)
355  REAL (R8) :: bdb02(na1),b0db2(na1),bdb0(na1)
356  REAL (R8) :: rout(neql,nteta), zout(neql,nteta)
357  REAL (R8) :: shafr_shift(na1),elong(na1),elong_u(na1),elong_l(na1),triang_l(na1),triang_u(na1),amid(na1)
358  REAL (R8) :: r_in, r_out, z_min, z_max, rz_min, rz_max
359  REAL (R8) :: r1,r2,r3,z1,z2,z3,a,b
360 
361  !locals
362  REAL (R8), ALLOCATABLE :: w_dj(:), yfofb(:)
363 
364  TYPE(type_parameters) :: parameters_spider
365 
366 
367  ALLOCATE ( w_dj(na1), yfofb(na1) )
368 
369 
370  w_dj = 0.e0_r8
371  yfofb = 0.e0_r8
372 
373 
374  ni_p = parameters_spider%neql
375  nj_p = parameters_spider%nteta
376 
377 
378 
379  rout = 0.e0_r8
380  zout = 0.e0_r8
381 
382 
383 
384  CALL spider2astra(rout,zout,rtor,btor,rho,roc,na1, &
385  g11,g22,g33,vr,vrs,slat,gradro,rocnew, &
386  mu,ipol,bmaxt,bmint,bdb02,b0db2,bdb0,droda, &
387  yreler,yli3,ni_p,nj_p,platok,cu_out,fp, &
388  pres_s,w_dj,yfofb)
389 
390  ftrap = yfofb
391 
392 
393  CALL output2d(neql,nteta,na1,rtor,ztor,rout,zout,shafr_shift,elong,elong_u,elong_l,triang_l,triang_u,amid)
394 
395 
396  DEALLOCATE ( w_dj, yfofb )
397 
398 
399 
400  RETURN
401 
402  END SUBROUTINE spider2output2
403 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
404 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
405 
406 
407 
408 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
409 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
410  SUBROUTINE output2d (neql,nteta,na1,rtor,ztor,rout,zout,shafr_shift,elong,elong_u,elong_l,triang_l,triang_u,amid)
411 
412  USE itm_constants
413 
414  IMPLICIT NONE
415 
416 
417  INTEGER :: na1, neql, nteta
418  INTEGER :: i, j, ndim2
419  INTEGER :: k1, k2
420 
421  REAL (R8) :: rout(neql,nteta), zout(neql,nteta), rtor,ztor
422  REAL (R8) :: shafr_shift(na1),elong(na1),elong_u(na1),elong_l(na1),triang_l(na1),triang_u(na1),amid(na1)
423  REAL (R8) :: r_in, r_out, z_in, z_out, z_min, z_max, rz_min, rz_max
424  REAL (R8) :: r1,r2,r3,z1,z2,z3,a,b
425 
426 
427  ndim2 = SIZE(rout, dim=2)
428 
429  DO i = 2, na1
430  r_in = minval(rout(i,1:ndim2))
431  r_out = maxval(rout(i,1:ndim2))
432  z_min = minval(zout(i,1:ndim2))
433  z_max = maxval(zout(i,1:ndim2))
434 
435  rz_min = 0.5d0*(r_in+r_out)
436  rz_max = 0.5d0*(r_in+r_out)
437 
438  DO j = 1, ndim2
439  IF (zout(i,j).EQ.z_min) THEN
440  k1 = 2
441  k2 = 2
442  IF (ndim2-j.LT.2) k2 = 2 - ndim2
443  IF (j.LT.2) k1 = 2 - ndim2
444  r1 = rout(i,j-k1)
445  r2 = rout(i,j)
446  r3 = rout(i,j+k2)
447  z1 = zout(i,j-k1)
448  z2 = zout(i,j)
449  z3 = zout(i,j+k2)
450 
451  a = ((z2-z3)/(r2-r3)-(z1-z2)/(r1-r2))/(r3-r1)
452  b = (z1-z2)/(r1-r2) - a*(r1+r2)
453 
454  rz_min = - b/2.d0/a
455  z_min = a*(rz_min**2-r2**2) + b*(rz_min - r2)+ z2
456 
457  IF (abs(1.d0-rz_min/rout(i,j)).ge.abs(1.d0-rout(i,1)/rout(i,2))) rz_min = rout(i,j)
458  END IF
459  IF (zout(i,j).EQ.z_max) THEN
460  k1 = 2
461  k2 = 2
462  IF (ndim2-j.LT.2) k2 = 2 - ndim2
463  IF (j.LT.2) k1 = 2 - ndim2
464  r1 = rout(i,j-k1)
465  r2 = rout(i,j)
466  r3 = rout(i,j+k2)
467  z1 = zout(i,j-k1)
468  z2 = zout(i,j)
469  z3 = zout(i,j+k2)
470 
471  a = ((z2-z3)/(r2-r3)-(z1-z2)/(r1-r2))/(r3-r1)
472  b = (z1-z2)/(r1-r2) - a*(r1+r2)
473 
474  rz_max = - b/2.d0/a
475  z_max = a*(rz_max**2-r2**2) + b*(rz_max - r2)+ z2
476  IF (abs(1.d0-rz_max/rout(i,j)).ge.abs(1.d0-rout(i,1)/rout(i,2))) rz_max = rout(i,j)
477  END IF
478  IF (rout(i,j).EQ.r_out) THEN
479  k1 = 2
480  k2 = 2
481  IF (ndim2-j.LT.2) k2 = 2 - ndim2
482  IF (j.LT.2) k1 = 2 - ndim2
483  r1 = rout(i,j-k1)
484  r2 = rout(i,j)
485  r3 = rout(i,j+k2)
486  z1 = zout(i,j-k1)
487  z2 = zout(i,j)
488  z3 = zout(i,j+k2)
489 
490  a = ((r2-r3)/(z2-z3)-(r1-r2)/(z1-z2))/(z3-z1)
491  b = (r1-r2)/(z1-z2) - a*(z1+z2)
492 
493  z_out = - b/2.d0/a
494  r_out = a*(z_out**2-z2**2) + b*(z_out-z2)+ r2
495  IF (abs(1.d0-z_out/zout(i,j)).ge.abs(1.d0-zout(i,1)/zout(i,2))) z_out = zout(i,j)
496  END IF
497  IF (rout(i,j).EQ.r_in) THEN
498  k1 = 2
499  k2 = 2
500  IF (ndim2-j.LT.2) k2 = 2 - ndim2
501  IF (j.LT.2) k1 = 2 - ndim2
502  r1 = rout(i,j-k1)
503  r2 = rout(i,j)
504  r3 = rout(i,j+k2)
505  z1 = zout(i,j-k1)
506  z2 = zout(i,j)
507  z3 = zout(i,j+k2)
508 
509  a = ((r2-r3)/(z2-z3)-(r1-r2)/(z1-z2))/(z3-z1)
510  b = (r1-r2)/(z1-z2) - a*(z1+z2)
511 
512  z_in = - b/2.d0/a
513  r_in = a*(z_in**2-z2**2) + b*(z_in-z2)+ r2
514  IF (abs(1.d0-z_in/zout(i,j)).ge.abs(1.d0-zout(i,1)/zout(i,2))) z_in = zout(i,j)
515  END IF
516  END DO
517 
518  amid(i) = 0.5d0*(r_out-r_in)
519  shafr_shift(i) = 0.5d0*(r_in+r_out)-rtor
520  elong(i) = (z_max-z_min)/(r_out-r_in)
521  elong_u(i) = 2.d0*abs(z_max-ztor)/(r_out-r_in)
522  elong_l(i) = 2.d0*abs(ztor-z_min)/(r_out-r_in)
523  triang_l(i) = (0.5d0*(r_in+r_out)-rz_min)/amid(i)
524  triang_u(i) = (0.5d0*(r_in+r_out)-rz_max)/amid(i)
525  END DO
526  shafr_shift(1) = rout(1,1)-rtor
527  elong(1) = elong(2)
528  triang_l(1) = 0.d0
529  triang_u(1) = 0.d0
530  amid(1) = 0.d0
531  shafr_shift(na1) = 0.d0
532 
533  RETURN
534 
535  END SUBROUTINE output2d
536 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
537 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
538 
539 end module mod_spider
mod_spider::input2spider2
subroutine input2spider2(neql, nteta, nbnd, rzbnd, na1, eqpf, eqff, fp, ipl, rtor, btor, rho, pres, cu, nstep, key_equil, PARAMETERS_SPIDER)
Definition: mod_spider.f90:137
savepsi
subroutine savepsi
Definition: interfac.f:2970
mod_spider
Definition: mod_spider.f90:1
put_tim
subroutine put_tim(dt, time)
Definition: B_eqb.f:944
mod_spider::output2d
subroutine output2d(neql, nteta, na1, rtor, ztor, rout, zout, shafr_shift, elong, elong_u, elong_l, triang_l, triang_u, amid)
Definition: mod_spider.f90:410
kpr_calc
subroutine kpr_calc(kpr_xx)
Definition: Spider_call.f:215
mod_spider::spider2output2
subroutine spider2output2(PARAMETERS_SPIDER, neql, nteta, rtor, ztor, btor, rho, roc, na1, g11, g22, g33, vr, vrs, slat, gradro, rocnew, mu, ipol, bmaxt, bmint, bdb02, b0db2, bdb0, droda, yreler, yli3, platok, cu_out, fp, pres_s, rout, zout, shafr_shift, elong, elong_u, elong_l, triang_l, triang_u, amid, ftrap)
Definition: mod_spider.f90:331
mod_spider::spider_run2
subroutine spider_run2(PARAMETERS_SPIDER)
Definition: mod_spider.f90:266
cur_avg
subroutine cur_avg
Definition: com_sub.f:813
wrb
subroutine wrb
Definition: B_wrd.f:1
wr_spik
subroutine wr_spik
Definition: B_wrd.f:646
parameters::type_parameters
Definition: parameters.f90:7
astra2spider
subroutine astra2spider(neql, nteta, nbnd, rzbnd, key_dmf,
Definition: interfac.f:1
put_key_con
subroutine put_key_con(k)
Definition: controller.f:1061
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
spider
subroutine spider(nstep, time, dt, key_dmf, k_grid, k_auto, k_fixfree,
Definition: Spider_call.f:1
put_name
subroutine put_name(name, ksym)
Definition: Spider_call.f:234
mod_spider::int_spider
subroutine int_spider(
Definition: mod_spider.f90:9
spider2astra
subroutine spider2astra(rout, zout, rtor, btor, rho, roc, na1,
Definition: interfac.f:369