From b4d877f6a3999d84e973f5209d930b89ad0c6e63 Mon Sep 17 00:00:00 2001 From: Daniela Farina Date: Wed, 10 Oct 2012 14:43:17 +0000 Subject: [PATCH] Modified Jcd outputs to be consistent with JINTRAC --- src/gray.f | 116 +++++++++++++++++++++++++++++++---------------------- 1 file changed, 67 insertions(+), 49 deletions(-) diff --git a/src/gray.f b/src/gray.f index 94acbe8..f71ee97 100644 --- a/src/gray.f +++ b/src/gray.f @@ -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= [T] , b2av= [T^2] , rbav=/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 = /(|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 = /B0 +c ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / +c ratio_pltor = Jcd_||/J_phi Jcd_|| = + 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/ -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 || |Bmx| |Bmn|'// - .' Area Vol |I_pl| qq fc' + .' Area Vol |I_pl| 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