ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
Urs.f
Go to the documentation of this file.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2  real*8 function funppp(psi)
3 
4  implicit real*8(a-h,o-z)
5  include 'urs.inc'
6  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
7 
8  if(psi.lt.0.d0) then
9  funppp =0.d0
10  return
11  endif
12  if(psi.gt.1.d0) then
13  funppp =ppp(nurs)
14  return
15  endif
16 
17  i=1+psi*(nurs-1)
18 
19  if(i.eq.nurs) i=nurs-1
20  if(i.lt.1) i=1
21 
22  dpsip=psit(i+1)-psi
23  dpsim=psi-psit(i)
24 
25  zppp=(ppp(i)*dpsip+ppp(i+1)*dpsim)/(dpsip+dpsim)
26 
27  funppp = zppp
28 
29  return
30  end
31 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32  real*8 function funfff(psi)
33 
34  implicit real*8(a-h,o-z)
35  include 'urs.inc'
36  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
37 
38  if(psi.lt.0.d0) then
39  funfff =0.d0
40  return
41  endif
42  if(psi.gt.1.d0) then
43  funfff =fff(nurs)
44  return
45  endif
46 
47  i=1+psi*(nurs-1)
48 
49  if(i.eq.nurs) i=nurs-1
50  if(i.lt.1) i=1
51 
52  dpsip=psit(i+1)-psi
53  dpsim=psi-psit(i)
54 
55  zfff=(fff(i)*dpsip+fff(i+1)*dpsim)/(dpsip+dpsim)
56 
57  funfff = zfff
58 
59  return
60  end
61 
62 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
63  real*8 function funcur(ri,psi)
64 
65  implicit real*8(a-h,o-z)
66  include 'urs.inc'
67 
68  if(psi.lt.0.d0) then
69  funcur =0.d0
70  return
71  endif
72 
73  if(psi.gt.1.d0) then
74  funcur =purs(nurs)*ri+furs(nurs)/ri+wurs(nurs)*ri**3
75  return
76  endif
77 
78  i=1+psi*(nurs-1)
79 
80  if(i.ge.nurs) i=nurs-1
81  if(i.lt.1) i=1
82 
83  dpsip=psit(i+1)-psi
84  dpsim=psi-psit(i)
85 
86  pp=(purs(i)*dpsip+purs(i+1)*dpsim)/(dpsip+dpsim)
87  fp=(furs(i)*dpsip+furs(i+1)*dpsim)/(dpsip+dpsim)
88  wp=(wurs(i)*dpsip+wurs(i+1)*dpsim)/(dpsip+dpsim)
89 
90  funcur = ri*pp + fp/ri + wp*ri**3
91 
92  return
93  end
94 
95 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
96 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
97  real*8 function funcur_p(ri,psi)
98 
99  implicit real*8(a-h,o-z)
100  include 'urs.inc'
101 
102  if(psi.lt.0.d0) then
103  funcur_p =0.d0
104  return
105  endif
106 
107  if(psi.gt.1.d0) then
108  funcur_p =purs(nurs)*ri
109  return
110  endif
111 
112  i=1+psi*(nurs-1)
113 
114  if(i.ge.nurs) i=nurs-1
115  if(i.lt.1) i=1
116 
117  dpsip=psit(i+1)-psi
118  dpsim=psi-psit(i)
119 
120  pp=(purs(i)*dpsip+purs(i+1)*dpsim)/(dpsip+dpsim)
121  !!!fp=(furs(i)*dpsip+furs(i+1)*dpsim)/(dpsip+dpsim)
122 
123  funcur_p = ri*pp
124 
125  return
126  end
127 
128 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
129 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
130  real*8 function funcur_f(ri,psi)
131 
132  implicit real*8(a-h,o-z)
133  include 'urs.inc'
134 
135  if(psi.lt.0.) then
136  funcur_f =0.d0
137  return
138  endif
139 
140  if(psi.gt.1.d0) then
141  funcur_f =furs(nurs)/ri
142  return
143  endif
144 
145  i=1+psi*(nurs-1)
146 
147  if(i.ge.nurs) i=nurs-1
148  if(i.lt.1) i=1
149 
150  dpsip=psit(i+1)-psi
151  dpsim=psi-psit(i)
152 
153  !!!pp=(purs(i)*dpsip+purs(i+1)*dpsim)/(dpsip+dpsim)
154  fp=(furs(i)*dpsip+furs(i+1)*dpsim)/(dpsip+dpsim)
155 
156  funcur_f = fp/ri
157 
158  return
159  end
160 
161 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
162 
163 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
164 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
165 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
166 
167  subroutine taburs(ien,coin,nursb)
168 
169  use ppf_modul
170 
171  include 'double.inc'
172  include 'urs.inc'
173  parameter(nursp4=nursp+4,nursp6=nursp4*6)
174  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
175  common/comsav/ alf0n
176  common/com_flag/kastr
177  !dimension pstab(nursp),pptab(nursp),fptab(nursp),wptab(nursp)
178  dimension wptab(nursp)
179  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
180  real*8 cwk(4)
181 
182 ! if(.not.allocated(pstab))
183 ! %allocate( pstab(nutab), pptab(nutab), fptab(nutab) )
184 
185 c!!!!!!table for fast calculation P' and FF' as funct. of PSI_nor
186 c!!!!!! P' - pptab
187 c!!!!!! FF'- fptab
188 c!!!!!! PSI_nor - pstab
189 
190  if(kastr.eq.1) then
191  alw0p=0.d0
192  alw1p=0.d0
193  alw2p=0.d0
194  endif
195 
196  if(ien.eq.0) then !!!!!!!!!!!!!!!!!!!!!!
197 
198  nurs=3999
199 
200 c if nursb<0 p'and ff'are assumed to be defined as tab.(file 'tabppf.dat'
201 
202 
203 c write(6,*) nurs
204 c
205 c write(6,*) alf0p
206 c write(6,*) alf1p
207 c write(6,*) alf2p
208 c
209 c write(6,*) bet0f
210 c write(6,*) bet1f
211 c write(6,*) bet2f
212 
213  psit(1)=0.d0
214  dpsi=1.d0/(nurs-1.d0)
215 
216  do i=2,nurs-1
217 
218  zpsi=psit(i-1)+dpsi
219  psit(i)=zpsi
220  enddo
221 
222  psit(nurs)=1.d0
223 
224  if(nursb.gt.0) then
225 
226  !write(6,*) 'taburs:initial calculation P and FF'
227 
228  do 10 i=1,nurs
229 
230  zpsi=psit(i)
231  purs(i)=funpp(zpsi)
232  furs(i)=funfp(zpsi)
233  wurs(i)=funwp(zpsi)
234 
235  !write(6,*) 'f/p,ps:',purs(i)/(furs(i)+1.d-15),zpsi,i
236  !pause ' '
237 
238  10 continue
239 
240  else
241 
242 ! write(fname,'(a,a)') path(1:kname),'tabppf.dat'
243 ! open(1,file=fname,form='formatted')
244 ! !open(1,file='tabppf.dat')
245 
246 ! read(1,*) nutab
247 
248 ! do i=1,nutab
249 
250 ! read(1,*) pstab(i),pptab(i),fptab(i)
251 
252 ! enddo
253 
254  do i=1,nutab
255 
256  pstab(i)=pstab(i)/pstab(nutab)
257 
258  enddo
259 
260 ! close(1)
261 
262  key_int=0
263  if(key_int.eq.0) then
264 
265  do i=2,nurs-1
266  zpsi=1.d0-psit(i)
267  do j=1,nutab-1
268  if(zpsi.ge.pstab(j) .AnD. zpsi.lt.pstab(j+1)) then
269  dpsip=pstab(j+1)-zpsi
270  dpsim=zpsi-pstab(j)
271  purs(i)=( pptab(j)*dpsip+pptab(j+1)*dpsim )/(dpsip+dpsim)
272  furs(i)=( fptab(j)*dpsip+fptab(j+1)*dpsim )/(dpsip+dpsim)
273  endif
274  enddo
275  enddo
276 
277  else
278 
279  CALL e01baf(nutab,pstab,pptab,rrk,cck,
280  * nutab+4,wrk,6*nutab+16,ifail)
281 
282  if(ifail.ne.0) write(*,*) 'ifail=',ifail
283  do i=2,nurs-1
284 
285  zpsi=1.d0-psit(i)
286 
287  CALL e02bcf(nutab+4,rrk,cck,zpsi,0,cwk,ifail)
288  if(ifail.ne.0) write(*,*) 'ifail=',ifail
289 
290  purs(i)=cwk(1)
291 
292  enddo
293 
294  CALL e01baf(nutab,pstab,fptab,rrk,cck,
295  * nutab+4,wrk,6*nutab+16,ifail)
296  if(ifail.ne.0) write(*,*) 'ifail=',ifail
297 
298  do i=2,nurs-1
299 
300  zpsi=1.d0-psit(i)
301 
302  CALL e02bcf(nutab+4,rrk,cck,zpsi,0,cwk,ifail)
303  if(ifail.ne.0) write(*,*) 'ifail=',ifail
304 
305  furs(i)=cwk(1)
306 
307  enddo
308 
309  endif
310 
311  purs(1)=pptab(nutab)
312  purs(nurs)=pptab(1)
313 
314  furs(1)=fptab(nutab)
315  furs(nurs)=fptab(1)
316 
317  !!>>>>>>>>>>>>>
318  do i=1,nurs
319  zpsi=psit(i)
320  wurs(i)=funwp(zpsi)
321  !wurs(i)=0.d0
322  enddo
323  !!>>>>>>>>>>>>>
324 
325  endif
326  ! open(1,file='tabus.wr')
327  ! write(1,*) nurs
328  ! write(1,*)(psit(i),i=1,nurs)
329  ! write(1,*)(purs(i),i=1,nurs)
330  ! write(1,*)(furs(i),i=1,nurs)
331  ! close(1)
332 
333  elseif(ien.eq.1) then !!!!!!!!!!
334 
335 
336  do i=1,nurs
337  purs(i)=purs(i)*coin
338  enddo
339 
340  !write(6,*) 'taburs:purs recalculation'
341  !write(6,*) 'coin=',coin
342 
343 
344  endif !!!!!!!!!!!!!!!!!!!!!!!!!
345 
346 
347  ppp(1)=0.d0
348  fff(1)=0.d0
349  www(1)=0.d0
350  dpsi=1.d0/(nurs-1.d0)
351 
352  do 20 i=2,nurs
353 
354  ppp(i)=ppp(i-1) +(purs(i-1)+purs(i))*dpsi*0.5d0
355  fff(i)=fff(i-1) +(furs(i-1)+furs(i))*dpsi*0.5d0
356  www(i)=www(i-1) +(wurs(i-1)+wurs(i))*dpsi*0.5d0
357 
358 ccc write(6,*) 'i,ppp(i)',i,ppp(i)
359  20 continue
360 
361  !deallocate( pstab, pptab, fptab )
362 
363  return
364  end
365 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
366  subroutine tabper
367  include 'double.inc'
368  include 'urs.inc'
369  parameter(nursp4=nursp+4,nursp6=nursp4*6)
370 
371  do i=1,nurs
372  purs(i)=purs(i)*1.d3
373  furs(i)=furs(i)
374  enddo
375 
376  return
377  end
378 
379 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
380 
381 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
382  subroutine tabnor(cnor)
383  include 'double.inc'
384  include 'urs.inc'
385  parameter(nursp4=nursp+4,nursp6=nursp4*6)
386  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
387 
388  do i=1,nurs
389  purs(i)=cnor*purs(i)
390  furs(i)=cnor*furs(i)
391  ppp(i)=cnor*ppp(i)
392  fff(i)=cnor*fff(i)
393  enddo
394 
395  return
396  end
397 
398 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
399 
400  real*8 function funpp(psi)
401 c
402  implicit real*8(a-h,o-z)
403  include 'urs.inc'
404  epss=1.d-10
405 
406  if(psi.lt.0.d0) then
407  funpp=0.d0
408  return
409  endif
410 
411  psm = 1.d0-psi + epss
412 
413  zpp = dabs( 1.d0 -(dabs(psm))**alf1p ) + epss
414  pp = alf0p*( zpp**alf2p )
415 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
416 !! Profile for TCV, shot 10396, time = 0.24
417 !
418 ! c0 = 0.0d0
419 ! c1 = 0.730930d0
420 ! c2 = 0.d0
421 ! c3 = 0.d0
422 !
423 ! pp = alf0p*( c0 + c1*psi + c2*psi*2 + c3*psi**3 )
424 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
425 c-------
426 c. cbi=0.86d0
427 c. cr0=8.025d0
428 c. cnp=-1.4d0
429 c. cnf=-1.4d0
430 c. cal=2.d0
431 c
432 c. pp=(cbi/cr0)*( dexp( cnp*(1-psm**cal) ) - 1. )/
433 c. / ( dexp( cnp ) - 1. )
434 c-------
435 !------------------------------------------------
436 ! Profile for TCV, shot 20333, time = 0.500
437 !
438 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
439  del=0.0d0
440  x=(psi-del)/(1.d0-del)
441 !------------------------------------------------
442 ! TCV, shot 20333, time=0.500
443  c0 = 0.0d0
444  c1 = 1.7791d0
445  c2 = 0.d0
446  c3 = 0.d0
447 !---------------------------------------------------------------------
448 !
449 ! if(x.gt.0.d0) then
450 ! pp = alf0p*( c0 + c1*x + c2*x**2 + c3*x**3 )
451 ! !pp = alf0p*( c0 + c1*psi + c2*psi**2 + c3*psi**3 )
452 ! else
453 ! pp = 0.d0
454 ! endif
455 !
456 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
457  funpp=pp
458 c-------
459  return
460  end
461 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
462 
463  real*8 function funfp(psi)
464 c
465  implicit real*8(a-h,o-z)
466  include 'urs.inc'
467 
468  epss=1.d-10
469 
470  if(psi.lt.0.d0) then
471  funfp=0.d0
472  return
473  endif
474 
475  psm = 1.d0- psi + epss
476 
477  zzfp = dabs( 1.d0 - psm**bet1f ) + epss
478  fp = bet0f*( zzfp**bet2f )
479 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
480 ! Profile for TCV, shot 10396, time = 0.24
481 !
482 ! c0 = 0.0d0
483 ! c1 =-0.147740d0
484 ! c2 = 2.542300d0
485 ! c3 = 0.d0
486 !
487 ! fp = bet0f*( c0 + c1*psi + c2*psi**2 + c3*psi**3 )
488 c----------------------------------------------------------
489 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
490 c-------
491 c. fp= 1.d0
492 c-------
493 c. cbi=0.86d0
494 c. cr0=8.025d0
495 c. cnp=-1.4d0
496 c. cnf=-1.4d0
497 c. cal=2.d0
498 c
499 c. fp=(1.-cbi)*cr0*( dexp( cnf*(1-psm**cal) ) -1. )/
500 c. / ( dexp( cnf ) - 1. )
501 c-------
502 !------------------------------------------------
503  del=0.0d0
504  x=(psi-del)/(1.d0-del)
505 ! Profile for TCV, shot 20333, time = 0.500
506  c0 = 0.0d0
507  c1 = 2.467100d0
508  c2 =-1.473200d0
509  c3 = 0.d0
510 !----------------------------------------------------------------------
511 !
512 ! if(x.gt.0.d0) then
513 ! !fp = bet0f*( c0 + c1*psi + c2*psi**2 + c3*psi**3 )
514 ! fp = bet0f*( c0 + c1*x + c2*x**2 + c3*x**3 )
515 ! else
516 ! fp = 0.d0
517 ! endif
518 !
519 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
520  funfp=fp
521 c-------
522  return
523  end
524 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
525 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
526 
527  real*8 function tabw(psi)
528 
529  implicit real*8(a-h,o-z)
530  include 'urs.inc'
531 
532  if(psi.lt.0.d0) then
533  wp =0.d0
534  return
535  endif
536 
537  if(psi.gt.1.d0) then
538  wp =wurs(nurs)
539  return
540  endif
541 
542  i=1+psi*(nurs-1)
543 
544  if(i.eq.nurs) i=nurs-1
545 
546  if(i.lt.1) i=1
547 
548  dpsip=psit(i+1)-psi
549  dpsim=psi-psit(i)
550 
551  wp=(wurs(i)*dpsip+wurs(i+1)*dpsim)/(dpsip+dpsim)
552 
553  tabw = wp
554 
555  return
556  end
557 
558 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
559 
560  real*8 function tabp(psi)
561 
562  implicit real*8(a-h,o-z)
563  include 'urs.inc'
564 
565  if(psi.lt.0.d0) then
566  pp =0.d0
567  return
568  endif
569 
570  if(psi.gt.1.d0) then
571  pp =purs(nurs)
572  return
573  endif
574 
575  i=1+psi*(nurs-1)
576 
577  if(i.eq.nurs) i=nurs-1
578 
579  if(i.lt.1) i=1
580 
581  dpsip=psit(i+1)-psi
582  dpsim=psi-psit(i)
583 
584  pp=(purs(i)*dpsip+purs(i+1)*dpsim)/(dpsip+dpsim)
585 
586  tabp = pp
587 
588  return
589  end
590 
591 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
592 
593 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
594 
595  real*8 function tabf(psi)
596 
597  implicit real*8(a-h,o-z)
598  include 'urs.inc'
599 
600  if(psi.lt.0.d0) then
601 
602  fp =0.d0
603 
604  return
605 
606  endif
607 
608  if(psi.gt.1.d0) then
609 
610  fp =furs(nurs)
611 
612  return
613 
614  endif
615 
616  i=1+psi*(nurs-1)
617 
618  if(i.eq.nurs) i=nurs-1
619 
620  if(i.lt.1) i=1
621 
622  dpsip=psit(i+1)-psi
623  dpsim=psi-psit(i)
624 
625  fp=(furs(i)*dpsip+furs(i+1)*dpsim)/(dpsip+dpsim)
626 
627  tabf = fp
628 
629  return
630  end
631 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
632 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
633 
634  real*8 function funwp(psi)
635 
636  implicit real*8(a-h,o-z)
637  include 'urs.inc'
638 
639  epss=1.d-10
640  pi=3.14159265359d0
641 
642  if(psi.lt.0.d0) then
643  funwp=0.d0
644  return
645  endif
646 
647  psm = 1.d0 - psi + epss
648 
649  zwp = dabs( 1.d0 -(dabs(psm))**alw1p ) + epss
650  wp = alw0p*( zwp**alw2p )
651  !wp = -alw0p*dsin( 2.d0*pi*psm )
652 
653  funwp = wp
654  !funwp = 0.d0
655 
656  return
657  end
658 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
e02bcf
subroutine e02bcf(NCAP7, K, C, X, LEFT, S, IFAIL)
Definition: NAG.f:2761
e01baf
subroutine e01baf(M, X, Y, K, C, LCK, WRK, LWRK, IFAIL)
Definition: NAG.f:2511
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
ppp
subroutine ppp
Definition: ppplib.f:1