istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

Modified Jcd outputs to be consistent with JINTRAC

Browse Source
This commit is contained in:
Daniela Farina 2012-10-10 14:43:17 +00:00
parent 8348f52783
commit b4d877f6a3
1 changed files with 67 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

116
src/gray.f
View File

@ -508,7 +508,7 @@ c
write(7,*)'#Icd Pa Jphimx dPdVmx '//
.'rhotj rhotjava rhotp rhotpav '//
.'drhotjava drhotpav stmx polpsi polchi index_rt'
write(48,*) '#psi rhot Jphi Jcda dPdV Icdins Pins P% index_rt'
write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
else
@ -2128,8 +2128,9 @@ c
dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1)
dimension cbmx(nnintp,4),cbmn(nnintp,4),crbav(nnintp,4)
dimension cvol(nnintp,4),crri(nnintp,4),carea(nnintp,4)
dimension cratj1(nnintp,4),cratj2(nnintp,4),cfc(nnintp,4)
dimension vratj1(nnintp),vratj2(nnintp),crhotq(nnintp,4)
dimension cratjpl(nnintp,4),cratja(nnintp,4),cratjb(nnintp,4)
dimension cfc(nnintp,4),crhotq(nnintp,4)
dimension vratjpl(nnintp),vratja(nnintp),vratjb(nnintp)
dimension alam(nlam),fhlam(nnintp,nlam)
dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam)
dimension tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4))
@ -2144,7 +2145,7 @@ c
common/pstab/rpstab
common/flav/vvol,rri,rbav,bmxpsi,bmnpsi
common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc
common/cratj/cratj1,cratj2
common/cratj/cratja,cratjb,cratjpl
common/crhotq/crhotq
common/cnt/rup,zup,rlw,zlw
common/bound/zbmin,zbmax
@ -2183,6 +2184,7 @@ c computation of flux surface averaged quantities
dvdpsi=2.0d0*pi*anorm
dadpsi=2.0d0*pi/abs(btaxis)
ratio_cdator=abs(btaxis/btrcen)
ratio_cdbtor=1.0d0
ratio_pltor=1.0d0
qq=1.0d0
fc=1.0d0
@ -2199,8 +2201,9 @@ c computation of flux surface averaged quantities
rri(1)=rmaxis
varea(1)=0.0d0
vvol(1)=0.0d0
vratj1(1)=ratio_pltor
vratj2(1)=ratio_cdator
vratjpl(1)=ratio_pltor
vratja(1)=ratio_cdator
vratjb(1)=ratio_cdbtor
ffc(1)=fc
rhot2q=0.0d0
@ -2245,23 +2248,23 @@ c
dlp=sqrt((rcon(inc1)-rcon(inc))**2+
. (zcon(inc1)-zcon(inc))**2)
drc=(rcon(inc1)-rcon(inc))
c
c compute length, area and volume defined by psi=height^2
c
ph=0.5d0*(dla+dlb+dlp)
area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp)
area=area+sqrt(area2)
rzp=rcon(inc1)*zcon(inc1)
rz=rcon(inc)*zcon(inc)
volume=pi*(rzp+rz)*drc+volume
c
c compute line integral on the contour psi=height^2
c
rpsim=rcon(inc1)
zpsim=zcon(inc1)
call bfield(rpsim,zpsim,bphi,brr,bzz)
call tor_curr(rpsim,zpsim,ajphi)
c
btot=sqrt(bphi**2+brr**2+bzz**2)
bpoloid=sqrt(brr**2+bzz**2)
dlpv(inc)=dlp
@ -2284,20 +2287,20 @@ c
rpsim0=rpsim
bpoloid0=bpoloid
btot0=btot
c
c computation maximum/minimum B values on given flux surface
c
if(btot.le.bmmn) bmmn=btot
if(btot.ge.bmmx) bmmx=btot
end do
c
c bav=<B> [T] , b2av=<B^2> [T^2] , rbav=<B>/b_min
c anorm = int d l_p/B_p = dV/dpsi/(2pi)
c r2iav=<1/R^2> [m^-2] ,
c riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)),
c rri = <B>/(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1]
c currp = plasma current within psi=const
c
bbav=bbav/anorm
r2iav=r2iav/anorm
dvdpsi=2.0d0*pi*anorm
@ -2305,10 +2308,13 @@ c
b2av=b2av/anorm
vcurrp(jp)=ccj*currp
vajphiav(jp)=ajphiav/dadpsi
c
c area == varea, volume == vvol
c flux surface minor radius == (area/pi)^1/2
c
c area == varea, volume == vvol
c flux surface minor radius == (area/pi)^1/2
c ratio_cdator = Jcd_astra/J_phi Jcd_astra = <Jcd.B>/B0
c ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = <Jcd.B>/<B>
c ratio_pltor = Jcd_||/J_phi Jcd_|| = <Jcd.b>
pstab(jp)=height2
rpstab(jp)=height
vvol(jp)=abs(volume)
@ -2319,14 +2325,16 @@ c
bmnpsi(jp)=bmmn
rri(jp)=bav(jp)/abs(fpolv*r2iav)
ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen))
ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav))
ratio_pltor=abs(bbav*riav/(fpolv*r2iav))
vratj1(jp)=ratio_pltor
vratj2(jp)=ratio_cdator
vratjpl(jp)=ratio_pltor
vratja(jp)=ratio_cdator
vratjb(jp)=ratio_cdbtor
qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi))
qqv(jp)=qq
c
c computation of rhot from calculated q profile
c
rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1))
. /dble(nintp-1)
rhotqv(jp)=sqrt(rhot2q)
@ -2334,7 +2342,7 @@ c
c computation of fraction of circulating/trapped fraction fc, ft
c and of function H(lambda,rhop)
c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
c
fc=0.0d0
shlam=0.0d0
do l=nlam,1,-1
@ -2376,7 +2384,7 @@ c
end do
write(56,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
.' Area Vol |I_pl| <J_phi> qq fc'
.' Area Vol |I_pl| <J_phi> qq fc ratio_cdbtor'
qqv(1)=qqv(2)
vajphiav(1)=vajphiav(2)
@ -2392,13 +2400,13 @@ c
. bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp),
. vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)
end do
c
rarea=sqrt(varea(nintp)/pi)
c
c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi
c used for computations of dP/dV and J_cd
c spline coefficients of rhot
c
iopt=0
call difcs(rpstab,vvol,nintp,iopt,cvol,ier)
iopt=0
@ -2410,9 +2418,11 @@ c
iopt=0
call difcs(rpstab,bmnpsi,nintp,iopt,cbmn,ier)
iopt=0
call difcs(rpstab,vratj1,nintp,iopt,cratj1,ier)
call difcs(rpstab,vratja,nintp,iopt,cratja,ier)
iopt=0
call difcs(rpstab,vratj2,nintp,iopt,cratj2,ier)
call difcs(rpstab,vratjb,nintp,iopt,cratjb,ier)
iopt=0
call difcs(rpstab,vratjpl,nintp,iopt,cratjpl,ier)
iopt=0
call difcs(rpstab,varea,nintp,iopt,carea,ier)
iopt=0
@ -4265,19 +4275,21 @@ c
c
c
c
subroutine ratioj(rpsi,ratj1i,ratj2i)
subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli)
implicit real*8 (a-h,o-z)
parameter(nnintp=41)
dimension rpstab(nnintp),cratj1(nnintp,4),cratj2(nnintp,4)
dimension rpstab(nnintp)
dimension cratja(nnintp,4),cratjb(nnintp,4),cratjpl(nnintp,4)
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
common/cratj/cratj1,cratj2
common/cratj/cratja,cratjb,cratjpl
ip=int((nintp-1)*rpsi+1)
if(ip.eq.0) ip=1
if(ip.eq.nintp) ip=nintp-1
dps=rpsi-rpstab(ip)
ratj1i=spli(cratj1,nintp,ip,dps)
ratj2i=spli(cratj2,nintp,ip,dps)
ratjai=spli(cratja,nintp,ip,dps)
ratjbi=spli(cratjb,nintp,ip,dps)
ratjpli=spli(cratjpl,nintp,ip,dps)
return
end
c
@ -5705,8 +5717,8 @@ c
dimension rtab1(0:nndmx)
dimension ajphiv(nndmx),dpdv(nndmx)
dimension dvolt(nndmx),darea(nndmx)
dimension ratj1v(nndmx),ratj2v(nndmx)
dimension ajplv(nndmx),ajcdv(nndmx)
dimension ratjav(nndmx),ratjbv(nndmx),ratjplv(nndmx)
dimension ajplv(nndmx),ajcdav(nndmx),ajcdbv(nndmx)
dimension pins(nndmx),currins(nndmx),fi(nndmx)
parameter(llmx=21)
dimension isev(llmx)
@ -5771,9 +5783,10 @@ c
darea(it)=abs(areai1-areai0)
voli0=voli1
areai0=areai1
call ratioj(sqrt(psit),ratj1i,ratj2i)
ratj1v(it)=ratj1i
ratj2v(it)=ratj2i
call ratioj(sqrt(psit),ratjai,ratjbi,ratjpli)
ratjav(it)=ratjai
ratjbv(it)=ratjbi
ratjplv(it)=ratjpli
else
rhotv(it)=sqrt(psit)
area1=pi*psit1*rpam**2
@ -5926,8 +5939,9 @@ c of gaussian profile
do i=1,nd
dpdv(i)=facpds*dpdv(i)
ajphiv(i)=facjs*ajphiv(i)
ajcdv(i)=ajphiv(i)*ratj2v(i)
ajplv(i)=ajphiv(i)*ratj1v(i)
ajcdav(i)=ajphiv(i)*ratjav(i)
ajcdbv(i)=ajphiv(i)*ratjbv(i)
ajplv(i)=ajphiv(i)*ratjplv(i)
end do
call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx,
@ -5936,6 +5950,7 @@ c of gaussian profile
call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi,
. drhotjfi,drhopfi)
xps=rhopfi
call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx)
else
rhotjfi=0.0d0
rhopfi=0.0d0
@ -5943,6 +5958,9 @@ c of gaussian profile
drhotjfi=0.0d0
drhopfi=0.0d0
xps=rhopp
ratjamx=1.0d0
ratjbmx=1.0d0
ratjplmx=1.0d0
end if
iif1=iiv(1,1)
@ -5975,10 +5993,10 @@ c of gaussian profile
ajmxfi=0.0d0
dpdvmx=0.0d0
rhotjfi=1.0d0
rhotjav=1.0d0
rhotjava=1.0d0
rhotp=1.0d0
rhotpav=1.0d0
rhotjav=1.0d0
rhotjava=1.0d0
rhotp=1.0d0
rhotpav=1.0d0
drhotjfi=0.0d0
drhotjava=0.0d0
drhotp=0.0d0
@ -6000,7 +6018,7 @@ c dPdV [MW/m^3], Jcd [MA/m^2]
write(6,*)' '
write(6,*)'#beta alpha Icd Pa Jphimx dPdVmx '//
.'rhotj rhotjava rhotp rhotpav '//
.'drhotjava drhotpav stmx polpsi polchi index_rt'
.'drhotjava drhotpav ratjbmx stmx polpsi polchi index_rt'
write(6,99) betac,alfac,currtka,pabstot,ajmxfi,dpdvmx,
. rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav,
@ -6008,7 +6026,7 @@ c dPdV [MW/m^3], Jcd [MA/m^2]
write(7,99) currtka,pabstot,ajmxfi,dpdvmx,
. rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav,
. drhotjava,drhotpav,ratjbmx,
. stmx,psipol,chipol,real(index_rt)
do i=1,nd
@ -6021,7 +6039,7 @@ c dPdV [MW/m^3], Jcd [MA/m^2]
end if
pinsr=0.0d0
if(pabstot.gt.0) pinsr=pins(i)/pabstot
write(48,99) psin,rhotv(i),ajphiv(i),ajcdv(i),dpdv(i),
write(48,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i),
. currins(i),pins(i),pinsr,real(index_rt)
end do