diff --git a/src/gray-externals.f b/src/gray-externals.f index 37ae005..3049114 100644 --- a/src/gray-externals.f +++ b/src/gray-externals.f @@ -319,6 +319,8 @@ c common/nplr/anpl,anpr common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nprw/anprr,anpri +c + common/angles/alpha0,beta0 c x=ywrk(1,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 pt11=exp(-taujk) - write(604,99) stm,rrm,zzm,phideg, + write(604,99) alpha0,beta0,stm,rrm,zzm,phideg, . psi,rhot,dens11,tekev, . bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100, . 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(j.eq.nrayr) 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) end if c if(k.eq.nrayth) write(633,*) ' ' @@ -407,6 +409,7 @@ c return 99 format(30(1x,e16.8e3)) 111 format(3i5,16(1x,e16.8e3)) +112 format(2(1x,e16.8e3),3i5,16(1x,e16.8e3)) end c c @@ -418,17 +421,19 @@ c If(index_rt.eq.1) then - write(604,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '// - .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' + write(604,*)' #alpha0 beta0 sst R z phi psi rhot ne Te Br Bphi '// + .'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(609,*) ' #istep j k xt yt zt rt psin' 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(607,*)'#Icd Pa Jphip dPdVp '// + write(607,*)'#alpha0 beta0 Icd Pa Jphip dPdVp '// .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// .'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 @@ -691,7 +696,6 @@ c c read data for beam from file if ibeam>0 c if(ibeam.gt.0) then - filenmbm='graybeam.data' call read_beams(beamid,iox) else zrcsi=0.5d0*ak0*w0csi**2 @@ -921,7 +925,7 @@ c . ypolygC, xpolygD, ypolygD real*8 xvert(4), yvert(4) c - parameter(kspl=3,sspl=0.01) + parameter(kspl=1,sspl=0.01) c common/ibeam/ibeam common/filesn/filenmeqq,filenmprf,filenmbm @@ -2847,6 +2851,10 @@ c computation of flux surface averaged quantities ffc(1)=fc rhot2q=0.0d0 + + dadrhotv(1)=0.0d0 + dvdrhotv(1)=0.0d0 + C rup=rmaxis C rlw=rmaxis C zup=(zbmax+zmaxis)/2.0d0 @@ -4345,6 +4353,8 @@ c common/tete/tekev common/nplr/anpl,anpr common/parpl/brr,bphi,bzz,ajphi +c + common/angles/alpha0,beta0 c ui=(0.0d0,1.0d0) csth=anz0c @@ -4567,11 +4577,12 @@ c r0m=r0/1.0d2 z0m=z0/1.0d2 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 end if 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, . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one @@ -4594,6 +4605,7 @@ c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) +112 format(2(1x,e16.8e3),3i5,20(1x,e16.8e3)) end c c ray tracing initial conditions igrad=0 @@ -4618,6 +4630,8 @@ c common/tete/tekev common/nplr/anpl,anpr common/parpl/brr,bphi,bzz,ajphi +c + common/angles/alpha0,beta0 c csth=anz0c snth=sqrt(1.0d0-csth**2) @@ -4725,15 +4739,16 @@ c r0m=r0/1.0d2 z0m=z0/1.0d2 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 end if if(j.eq.1.and.k.eq.1) then write(617,99) zero,zero,zero,zero - write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, - . psinv,one,dens,tekev,brr,bphi,bzz, - . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, - . zero,zero,zero,zero,zero,zero,zero,one + write(604,99) alpha0,beta0,zero,r0m,z0m, + . atan2(y0m,x0m)*180.0d0/pi, + . psinv,one,dens,tekev,brr,bphi,bzz, + . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, + . zero,zero,zero,zero,zero,zero,zero,one end if end do end do @@ -4749,6 +4764,7 @@ c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) +112 format(2(1x,e16.8e3),3i5,20(1x,e16.8e3)) end @@ -4769,6 +4785,8 @@ c common/tete/tekev common/nplr/anpl,anpr common/parpl/brr,bphi,bzz,ajphi +c + common/angles/alpha0,beta0 do j=1,nrayr do k=1,nrayth @@ -4823,12 +4841,13 @@ c r0m=r0/1.0d2 z0m=z0/1.0d2 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 end if if(j.eq.1.and.k.eq.1) then 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, . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one @@ -4847,6 +4866,7 @@ c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) +112 format(2(1x,e16.8e3),3i5,20(1x,e16.8e3)) end c 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, . 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, . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp @@ -6892,8 +6912,8 @@ c dPdV [MW/m^3], Jcd [MA/m^2] end if pinsr=0.0d0 if(pabstot.gt.0) pinsr=pins(i)/pabstot - write(648,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i), - . currins(i),pins(i),pinsr,real(index_rt) + write(648,99) alpha0,beta0,psin,rhotv(i),ajphiv(i),ajcdbv(i), + . dpdv(i),currins(i),pins(i),pinsr,real(index_rt) end do c fill output arguments only if ipec=-1 if (ipec.eq.-1) then