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

gray-jintrac gray-externals.f: added dadrhotv(1) and dvdrhotv(1) initialization to 0 in fluxaverage, spline degree in read_beams set to 1

Browse Source
This commit is contained in:
Daniele Micheletti 2015-02-13 15:18:32 +00:00
parent 7a14671b2a
commit a78099b8b8
1 changed files with 41 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
src/gray-externals.f
View File

@ -319,6 +319,8 @@ c
common/nplr/anpl,anpr common/nplr/anpl,anpr
common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11
common/nprw/anprr,anpri common/nprw/anprr,anpri
c
common/angles/alpha0,beta0
c c
x=ywrk(1,j,k) x=ywrk(1,j,k)
y=ywrk(2,j,k) y=ywrk(2,j,k)
@ -382,7 +384,7 @@ c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1
tau11=taujk tau11=taujk
pt11=exp(-taujk) pt11=exp(-taujk)
write(604,99) stm,rrm,zzm,phideg, write(604,99) alpha0,beta0,stm,rrm,zzm,phideg,
. psi,rhot,dens11,tekev, . psi,rhot,dens11,tekev,
. bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100, . bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100,
. alpha11,tau11,pt11,ajcd11,dids11, . alpha11,tau11,pt11,ajcd11,dids11,
@ -398,7 +400,7 @@ c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
if(istpl.eq.istpl0) then if(istpl.eq.istpl0) then
if(j.eq.nrayr) then if(j.eq.nrayr) then
c if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then c if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
write(633,111) i,j,k,stm,xxm,yym,rrm,zzm, write(633,112) alpha0,beta0,i,j,k,stm,xxm,yym,rrm,zzm,
. psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt) . psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt)
end if end if
c if(k.eq.nrayth) write(633,*) ' ' c if(k.eq.nrayth) write(633,*) ' '
@ -407,6 +409,7 @@ c
return return
99 format(30(1x,e16.8e3)) 99 format(30(1x,e16.8e3))
111 format(3i5,16(1x,e16.8e3)) 111 format(3i5,16(1x,e16.8e3))
112 format(2(1x,e16.8e3),3i5,16(1x,e16.8e3))
end end
c c
c c
@ -418,17 +421,19 @@ c
If(index_rt.eq.1) then If(index_rt.eq.1) then
write(604,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '// write(604,*)' #alpha0 beta0 sst R z phi psi rhot ne Te Br Bphi '//
.'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' .'Bz Jphi N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt'
write(608,*) ' #istep j k xt yt zt rt psin' write(608,*) ' #istep j k xt yt zt rt psin'
write(609,*) ' #istep j k xt yt zt rt psin' write(609,*) ' #istep j k xt yt zt rt psin'
write(617,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(617,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
write(633,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' write(633,*) ' #alpha0 beta0 i j k sst x y R z psi tauv Npl '//
.'alpha index_rt'
write(612,*) ' #i sst psi w1 w2' write(612,*) ' #i sst psi w1 w2'
write(607,*)'#Icd Pa Jphip dPdVp '// write(607,*)'#alpha0 beta0 Icd Pa Jphip dPdVp '//
.'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
.'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp' .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp'
write(648,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' write(648,*) '#alpha0 beta0 psi rhot Jphi Jcdb dPdV Icdins '//
.'Pins P% index_rt'
else else
@ -691,7 +696,6 @@ c
c read data for beam from file if ibeam>0 c read data for beam from file if ibeam>0
c c
if(ibeam.gt.0) then if(ibeam.gt.0) then
filenmbm='graybeam.data'
call read_beams(beamid,iox) call read_beams(beamid,iox)
else else
zrcsi=0.5d0*ak0*w0csi**2 zrcsi=0.5d0*ak0*w0csi**2
@ -921,7 +925,7 @@ c
. ypolygC, xpolygD, ypolygD . ypolygC, xpolygD, ypolygD
real*8 xvert(4), yvert(4) real*8 xvert(4), yvert(4)
c c
parameter(kspl=3,sspl=0.01) parameter(kspl=1,sspl=0.01)
c c
common/ibeam/ibeam common/ibeam/ibeam
common/filesn/filenmeqq,filenmprf,filenmbm common/filesn/filenmeqq,filenmprf,filenmbm
@ -2847,6 +2851,10 @@ c computation of flux surface averaged quantities
ffc(1)=fc ffc(1)=fc
rhot2q=0.0d0 rhot2q=0.0d0
dadrhotv(1)=0.0d0
dvdrhotv(1)=0.0d0
C rup=rmaxis C rup=rmaxis
C rlw=rmaxis C rlw=rmaxis
C zup=(zbmax+zmaxis)/2.0d0 C zup=(zbmax+zmaxis)/2.0d0
@ -4345,6 +4353,8 @@ c
common/tete/tekev common/tete/tekev
common/nplr/anpl,anpr common/nplr/anpl,anpr
common/parpl/brr,bphi,bzz,ajphi common/parpl/brr,bphi,bzz,ajphi
c
common/angles/alpha0,beta0
c c
ui=(0.0d0,1.0d0) ui=(0.0d0,1.0d0)
csth=anz0c csth=anz0c
@ -4567,11 +4577,12 @@ c
r0m=r0/1.0d2 r0m=r0/1.0d2
z0m=z0/1.0d2 z0m=z0/1.0d2
if(j.eq.nrayr) then if(j.eq.nrayr) then
write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, write(633,112) alpha0,beta0,izero,j,k,zero,x0m,y0m,r0m,z0m,
. psinv,zero,anpl,zero,one . psinv,zero,anpl,zero,one
end if end if
if(j.eq.1.and.k.eq.1) then if(j.eq.1.and.k.eq.1) then
write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, write(604,99) alpha0,beta0,zero,r0m,z0m,
. atan2(y0m,x0m)*180.0d0/pi,
. psinv,one,dens,tekev,brr,bphi,bzz, . psinv,one,dens,tekev,brr,bphi,bzz,
. ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero,
. zero,zero,zero,zero,zero,zero,zero,one . zero,zero,zero,zero,zero,zero,zero,one
@ -4594,6 +4605,7 @@ c
return return
99 format(24(1x,e16.8e3)) 99 format(24(1x,e16.8e3))
111 format(3i5,20(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3))
112 format(2(1x,e16.8e3),3i5,20(1x,e16.8e3))
end end
c c
c ray tracing initial conditions igrad=0 c ray tracing initial conditions igrad=0
@ -4618,6 +4630,8 @@ c
common/tete/tekev common/tete/tekev
common/nplr/anpl,anpr common/nplr/anpl,anpr
common/parpl/brr,bphi,bzz,ajphi common/parpl/brr,bphi,bzz,ajphi
c
common/angles/alpha0,beta0
c c
csth=anz0c csth=anz0c
snth=sqrt(1.0d0-csth**2) snth=sqrt(1.0d0-csth**2)
@ -4725,15 +4739,16 @@ c
r0m=r0/1.0d2 r0m=r0/1.0d2
z0m=z0/1.0d2 z0m=z0/1.0d2
if(j.eq.nrayr) then if(j.eq.nrayr) then
write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, write(633,112) alpha0,beta0,izero,j,k,zero,x0m,y0m,r0m,z0m,
. psinv,zero,anpl,zero,one . psinv,zero,anpl,zero,one
end if end if
if(j.eq.1.and.k.eq.1) then if(j.eq.1.and.k.eq.1) then
write(617,99) zero,zero,zero,zero write(617,99) zero,zero,zero,zero
write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, write(604,99) alpha0,beta0,zero,r0m,z0m,
. psinv,one,dens,tekev,brr,bphi,bzz, . atan2(y0m,x0m)*180.0d0/pi,
. ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, . psinv,one,dens,tekev,brr,bphi,bzz,
. zero,zero,zero,zero,zero,zero,zero,one . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero,
. zero,zero,zero,zero,zero,zero,zero,one
end if end if
end do end do
end do end do
@ -4749,6 +4764,7 @@ c
return return
99 format(24(1x,e16.8e3)) 99 format(24(1x,e16.8e3))
111 format(3i5,20(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3))
112 format(2(1x,e16.8e3),3i5,20(1x,e16.8e3))
end end
@ -4769,6 +4785,8 @@ c
common/tete/tekev common/tete/tekev
common/nplr/anpl,anpr common/nplr/anpl,anpr
common/parpl/brr,bphi,bzz,ajphi common/parpl/brr,bphi,bzz,ajphi
c
common/angles/alpha0,beta0
do j=1,nrayr do j=1,nrayr
do k=1,nrayth do k=1,nrayth
@ -4823,12 +4841,13 @@ c
r0m=r0/1.0d2 r0m=r0/1.0d2
z0m=z0/1.0d2 z0m=z0/1.0d2
if(j.eq.nrayr) then if(j.eq.nrayr) then
write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, write(633,112) alpha0,beta0,izero,j,k,zero,x0m,y0m,r0m,z0m,
. psinv,zero,anpl,zero,one . psinv,zero,anpl,zero,one
end if end if
if(j.eq.1.and.k.eq.1) then if(j.eq.1.and.k.eq.1) then
write(617,99) zero,zero,zero,zero write(617,99) zero,zero,zero,zero
write(604,99) strfl11,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, write(604,99) alpha0,beta0,strfl11,r0m,z0m,
. atan2(y0m,x0m)*180.0d0/pi,
. psinv,one,dens,tekev,brr,bphi,bzz, . psinv,one,dens,tekev,brr,bphi,bzz,
. ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero,
. zero,zero,zero,zero,zero,zero,zero,one . zero,zero,zero,zero,zero,zero,zero,one
@ -4847,6 +4866,7 @@ c
return return
99 format(24(1x,e16.8e3)) 99 format(24(1x,e16.8e3))
111 format(3i5,20(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3))
112 format(2(1x,e16.8e3),3i5,20(1x,e16.8e3))
end end
c c
c ray power weigth coefficient q(j) c ray power weigth coefficient q(j)
@ -6877,7 +6897,7 @@ c dPdV [MW/m^3], Jcd [MA/m^2]
. drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol,
. real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp
write(607,99) currtka,pabstot,ajphip,dpdvp, write(607,99) alpha0,beta0,currtka,pabstot,ajphip,dpdvp,
. rhotjfi,rhotjava,rhotp,rhotpav, . rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol,
. real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp
@ -6892,8 +6912,8 @@ c dPdV [MW/m^3], Jcd [MA/m^2]
end if end if
pinsr=0.0d0 pinsr=0.0d0
if(pabstot.gt.0) pinsr=pins(i)/pabstot if(pabstot.gt.0) pinsr=pins(i)/pabstot
write(648,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i), write(648,99) alpha0,beta0,psin,rhotv(i),ajphiv(i),ajcdbv(i),
. currins(i),pins(i),pinsr,real(index_rt) . dpdv(i),currins(i),pins(i),pinsr,real(index_rt)
end do end do
c fill output arguments only if ipec=-1 c fill output arguments only if ipec=-1
if (ipec.eq.-1) then if (ipec.eq.-1) then