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

output files revised

Browse Source
This commit is contained in:
Lorenzo Figini 2012-06-29 13:53:20 +00:00
parent 5985462447
commit ee2832a2d7
2 changed files with 116 additions and 116 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Makefile
View File

@ -12,6 +12,8 @@ vpath %.f src
FC=gfortran
FFLAGS=-O3
DIRECTIVES = -DREVISION="'$(shell svnversion)'"
all: $(EXE)
# Build executable from object files
@ -28,7 +30,7 @@ itm_constants.o: itm_types.o
$(FC) $(FFLAGS) -c $<
gray.o:gray.f green_func_p.o
$(FC) $(FFLAGS) -c $<
$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $<
grayl.o:grayl.f
$(FC) $(FFLAGS) -c $^

228
src/gray.f
View File

@ -42,7 +42,7 @@ c second pass into plasma
call prfile
call vectinit2
call paraminit
call ic_rt0
call ic_rt2
call gray_integration
call after_gray_integration
pabstott=pabstott+pabstot
@ -54,7 +54,7 @@ c second pass into plasma
call prfile
call vectinit2
call paraminit
call ic_rt0
call ic_rt2
call gray_integration
call after_gray_integration
pabstott=pabstott+pabstot
@ -64,7 +64,7 @@ c second pass into plasma
999 continue
print*,' '
print*,' IERR = ', ierr
print*,'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3
write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3
stop
end
@ -235,6 +235,7 @@ c ierr=0
psjki(j,k,i)=psinv
rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0
if(j.eq.1) rrm11=rrm
if (iwarm.gt.0.and.i.gt.1) then
if(psinv.ge.0.and.psinv.le.1.0d0) then
@ -285,7 +286,7 @@ c print*,' '
else
powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2)
end if
print*,'Reflected power fraction =',powrfl
write(6,*) 'Reflected power fraction =',powrfl
iov(j,k)=3
yyrfl(j,k,1)=xvrfl(1)
yyrfl(j,k,2)=xvrfl(2)
@ -320,7 +321,7 @@ c print*,' '
psimin=psjki(1,1,i)
if(nray.gt.1) psimin=min(psimin,minval(psjki(2:nray,1:ktx,i)))
if(psimin.gt.1.0d0.and.rrm.gt.rcen.and.index_rt.gt.1)
if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1)
. istop=1
if(iovmin.eq.3) istop=1
@ -492,15 +493,16 @@ c
write(4,*)' #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(8,*) ' #istep j k xt yt zt rt '
write(9,*) ' #istep j k xt yt zt rt '
write(17,*) ' #sst ddr_11 ddr_Nr1 ddi_Nr1'
write(8,*) ' #istep j k xt yt zt rt psin'
write(9,*) ' #istep j k xt yt zt rt psin'
write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
write(33,*) ' #i j k sst x y R z psi tauv Npl alpha'
write(12,*) ' #i sst rhop w1 w2 w1a w2a'
write(12,*) ' #i sst psi w1 w2'
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'
c write(29,*) '#beta alfa qqom uuom vvom psiom chiom'//
c .' qqxm uuxm vvxm psixm chixm'
c write(28,*) '#beta alfa anpl -gam ffo ffx xe2om xe2xm'
else
If(index_rt.eq.3) then
@ -509,6 +511,7 @@ c write(28,*) '#beta alfa anpl -gam ffo ffx xe2om xe2xm'
write(9,*) ' '
write(17,*) ' '
write(12,*) ' '
write(48,*) ' '
end if
end if
@ -521,6 +524,8 @@ c
use green_func_p, only:Setup_SpitzFunc
implicit real*8 (a-h,o-z)
real*8 me
character*8 wdat
character*10 wtim
character*24 filenmeqq,filenmprf,filenmbm
parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
@ -816,47 +821,58 @@ c versus psi, rhop, rhot
if (iequil.eq.1) call surf_anal
if (iequil.eq.2.and.iprof.eq.1) then
nfil=78
open(nfil,file='tit.txt',status= 'unknown')
write(nfil,905) filenmeqq
write(nfil,907) filenmprf
write(nfil,911) fghz
write(nfil,914) btrcen, sgnbphi,sgniphi,icocos
write(nfil,900) nray, ktx, rmx
open(nfil,file='headers.txt',status= 'unknown')
call date_and_time(wdat,wtim)
write(nfil,916) wdat(1:4),wdat(5:6),wdat(7:8),
. wtim(1:2),wtim(3:4),wtim(5:6)
write(nfil,904) REVISION
if (iequil.eq.2) then
write(nfil,905) filenmeqq
else
if (iequil.eq.1) write(nfil,912) rr0,zr0,rpa,B0,q0,qa,alq
if (iequil.eq.0) write(nfil,'("# VACUUM CASE")')
end if
if (iprof.eq.1) then
write(nfil,907) filenmprf
else
write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeff
end if
write(nfil,911) fghz,p0mw,iox
write(nfil,902) x00,y00,z00
write(nfil,908) alfac,betac
if(ibeam.ge.1) write(nfil,909) filenmbm
if(ibeam.eq.0) write(nfil,903) w0xt,w0yt,pw0xt,pw0yt,awr
write(nfil,900) nray, ktx, rmx
write(nfil,901) igrad,iwarm,ilarm,ieccd,idst
if(ieccd.eq.1) write(nfil,912)
if(ieccd.eq.11) write(nfil,913)
write(nfil,904) p0mw,iox
write(nfil,915) dst,nstep
write(nfil,914) sgnbphi,sgniphi,icocos
write(nfil,906) factb,factt,factn,iscal
write(nfil,910) ipec,nnd,ipsinorm,sspl,psdbnd
c write(nfil,915) psi15,sqrt(psi15),rhot15
c write(nfil,920) psi2,sqrt(psi2),rhot2
write(nfil,910) sspl,psdbnd,nnd,ipec,ipsinorm
write(nfil,*) ' '
close(nfil)
end if
return
900 format('# Nray ktx rmx',2i5,1x,e12.5)
901 format('# igrad iwarm ilarm ieccd idst',6i5)
902 format('# X Y Z launching point (cm) : ',3(1x,e12.5))
903 format('# w0xt w0yt pw0xt pw0yt (cm) delta (deg) : ',5(1x,e12.5))
904 format('# P0 IOX ',(1x,e12.5,i5))
906 format('# fact_B fact_T fact_n iscal',(3(1x,e12.5),i5))
905 format('# EQUILIBRIUM CASE : ',a24)
907 format('# PROFILE file : ',a24)
909 format('# LAUNCHER CASE : ',a24)
910 format('# ipec nd ipsinorm sspl psdbnd : ',3i5,2(1x,e12.5))
914 format('# |BT0| , sgnB_phi, sgnI_phi, icocos : ',3(1x,e12.5),i5)
915 format('# psi_1.5 rhop_1.5 rhot_1.5 : ',3(1x,e12.5))
920 format('# psi_2 rhop_2 rhot_2 : ',3(1x,e12.5))
911 format('# fghz : ',(1x,e12.5))
912 format('# Cohen model')
913 format('# Momentum conservation model')
900 format('# Nray ktx rmx : ',2i5,1x,es12.5)
901 format('# igrad iwarm ilarm ieccd idst : ',6i5)
902 format('# X Y Z launching point (cm) : ',3(1x,es12.5))
903 format('# w0xt w0yt pw0xt pw0yt (cm) delta (deg) : ',5(1x,es12.5))
904 format('# GRAY revision : ',a)
905 format('# EQUILIBRIUM file : ',a24)
906 format('# fact_B fact_T fact_n iscal : ',(3(1x,es12.5),i5))
907 format('# PROFILES file : ',a24)
908 format('# alpha beta launching angles (deg) CYL : ',2(1x,es12.5))
909 format('# LAUNCHER file : ',a24)
910 format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5)
911 format('# fghz P0 IOX : ',2(1x,es12.5),i5)
912 format('# AN. EQUILIBRIUM R0 z0 a B0 q0 qa alq : ',7(1x,es12.5))
913 format('# AN. PROFILES ne0 aln1 aln2 Te0 Tea alt1 alt2 Zeff : '
. ,8(1x,es12.5))
914 format('# sgnB_phi sgnI_phi icocos : ',2(1x,es12.5),i5)
915 format('# dst nstep : ',1x,es12.5,i5)
916 format('# Date : ',a4,2('/',a2),1x,a2,2(':',a2))
99 format(20(1x,e16.8e3))
end
@ -871,7 +887,7 @@ c
c
c print circular magnetic surfaces iequil=1
c
write(71,*) '#i psi R z'
write(71,*) '#i psin R z'
dal=2.0d-2*pi
drn=0.1d0
do i=1,10
@ -1305,7 +1321,7 @@ c
call points_ox(rmaxis,zmaxis,rmop,zmop,psinop,info)
rmaxis=rmop
zmaxis=zmop
print'(a,2f8.4,e12.5)','O-point',rmop,zmop,psinop
print'(a,2f8.4,es12.5)','O-point',rmop,zmop,psinop
c
c search for X-point if ixp not = 0
c
@ -1315,7 +1331,7 @@ c
z10=zbmin
call points_ox(r10,z10,rxp,zxp,psinxp,info)
if(psinxp.ne.-1.0d0) then
print'(a,2f8.4,e12.5)','X-point',rxp,zxp,psinxp
print'(a,2f8.4,es12.5)','X-point',rxp,zxp,psinxp
rbmin=rxp
zbmin=zxp
delpsinox=(psinxp-psinop)
@ -1481,7 +1497,6 @@ c
tol = sqrt(dpmpar(1))
call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
if(info.gt.1) then
c print*,' subr points_tgo: info=',info
end if
rf=xvec(1)
zf=xvec(2)
@ -1522,7 +1537,7 @@ c
common/eqnn/nr,nz,npp,nintp
common/dens/dens,ddens
c
write(55,*) ' #psi rhop rhot ne Te q Jphi'
write(55,*) ' #psi rhot ne Te q Jphi'
psin=0.0d0
rhop=0.0d0
rhot=0.0d0
@ -1531,7 +1546,7 @@ c
te=temp(psin)
qq=qpsi(1)
c
write(55,111) psin,rhop,rhot,dens,te,qq,ajphi*1.d-6
write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6
c
nst=nr
do i=2,nst
@ -1550,7 +1565,7 @@ c
end if
rhot=frhotor(psin)
call tor_curr_psi(psin,ajphi)
write(55,111) psin,rhop,rhot,dens,te,qq,ajphi*1.d-6
write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6
end do
c
return
@ -2117,7 +2132,7 @@ c
c computation of flux surface averaged quantities
write(71,*)' #i rhop R z'
write(71,*)' #i psin R z'
nintp=nnintp
ninpr=(nintp-1)/10
@ -2333,7 +2348,7 @@ c
end do
end do
write(56,*)' #psi rhop rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
write(56,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
.' Area Vol |I_pl| <J_phi> qq fc'
qqv(1)=qqv(2)
@ -2346,7 +2361,7 @@ c
pstab(jp)=1.0d0
end if
rhot_eq=frhotor_eq(pstab(jp))
write(56,99) pstab(jp),rpstab(jp),rhot_eq,rhotqv(jp),
write(56,99) pstab(jp),rhot_eq,rhotqv(jp),
. bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp),
. vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)
end do
@ -2465,7 +2480,7 @@ c
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension iiv(jmx,kmx),iov(jmx,kmx),tau1v(jmx,kmx)
dimension iiv(jmx,kmx),iov(jmx,kmx)
common/iiv/iiv
common/iov/iov
@ -2938,7 +2953,6 @@ c
anpl=anv(1)*bv(1)+anv(2)*bv(2)+anv(3)*bv(3)
c
if(abs(anpl).gt.0.99d0) then
c print*,' |Nparallel| > 1 !', sqrt(an2),anpl
if(abs(anpl).le.1.05d0) then
ierr=97
else
@ -3632,7 +3646,7 @@ c
subroutine ic_gb
implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
@ -3706,13 +3720,11 @@ c
qi1=rci1-ui*ww1
qi2=rci2-ui*ww2
end if
c w01=sqrt(2.0d0*akinv/ww1)
c z01=-rci1/(rci1**2+ww1**2)
c w02=sqrt(2.0d0*akinv/ww2)
c z02=-rci2/(rci2**2+ww2**2)
c
c spostare nel do????
c
qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2
qqyy=qi1*sin(phic)**2+qi2*cos(phic)**2
@ -3724,15 +3736,6 @@ c
rciyy=dble(qqyy)
rcixy=dble(qqxy)/2.0d0
C d2ww1=-2.0d0*(dww1*rci1+ww1*drci1)
C d2ww2=-2.0d0*(dww2*rci2+ww2*drci2)
C d2rci1=2.0d0*(ww1*dww1-rci1*drci1)
C d2rci2=2.0d0*(ww2*dww2-rci2*drci2)
C dqi1=drci1-ui*dww1
C dqi2=drci2-ui*dww2
C d2qi1=d2rci1-ui*d2ww1
C d2qi2=d2rci2-ui*d2ww2
c
dqi1=-qi1**2
dqi2=-qi2**2
d2qi1=2*qi1**3
@ -3743,7 +3746,6 @@ c
d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2
d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2
d2qqxy=-(d2qi1-d2qi2)*sin(2.0d0*phic)
c
dwwxx=-dimag(dqqxx)
dwwyy=-dimag(dqqyy)
dwwxy=-dimag(dqqxy)/2.0d0
@ -3782,7 +3784,6 @@ c
dx0= x0t*csps+snps*(y0t*csth+z0t*snth)
dy0=-x0t*snps+csps*(y0t*csth+z0t*snth)
dz0= z0t*csth-y0t*snth
c
x0=x00+dx0
y0=y00+dy0
z0=z00+dz0
@ -3885,7 +3886,13 @@ c
write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
. -1.0d0,zero,zero,zero
end if
if(j.eq.1) write(17,99) zero,zero,zero,zero
if(j.eq.1.and.k.eq.1) then
write(17,99) zero,zero,zero,zero
write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
. zero,zero,zero,zero,
. zero,zero,zero,zero,one,zero,zero,
. zero,zero,one,zero,zero,zero,zero,one
end if
end do
end do
c
@ -4034,7 +4041,13 @@ c
write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
. -1.0d0,zero,zero,zero
end if
if(j.eq.1) write(17,99) zero,zero,zero,zero
if(j.eq.1.and.k.eq.1) then
write(17,99) zero,zero,zero,zero
write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
. zero,zero,zero,zero,
. zero,zero,zero,zero,one,zero,zero,
. zero,zero,one,zero,zero,zero,zero,one
end if
end do
end do
c
@ -4053,7 +4066,7 @@ c
subroutine ic_rt0
subroutine ic_rt2
implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0)
@ -4118,20 +4131,22 @@ c
vgradi=0.0d0
ddi=2.0d0*vgradi
c
if(j.eq.nray.or.j.eq.1) then
if(j.eq.nray) then
r0=sqrt(x0**2+y0**2)
x0m=x0/1.0d2
y0m=y0/1.0d2
r0m=r0/1.0d2
z0m=z0/1.0d2
if(j.eq.nray) write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m
* ,zero,zero,zero,zero,zero
if(j.eq.1.and.k.eq.1)
* write(31,111) izero,j,k,zero,x0m,y0m,r0m,z0m,zero,zero
* ,zero,zero,zero
write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
. -1.0d0,zero,zero,zero
end if
if(k.eq.1.and.j.eq.nray)
* write(27,99) zero,x0,y0,z0,r0,dd,vgradi
if(j.eq.1.and.k.eq.1) then
write(17,99) zero,zero,zero,zero
write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
. zero,zero,zero,zero,
. zero,zero,zero,zero,one,zero,zero,
. zero,zero,one,zero,zero,zero,zero,one
end if
end do
end do
c
@ -5627,8 +5642,7 @@ c
end if
if(.not.(iproj.eq.0.and.j.eq.1))
. write(nfile,111) istep,j,k,xti,yti,zti,rti,
. sqrt(psinv11)
. write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
c
if(rti.ge.rtimx.and.j.eq.nray) rtimx=rti
if(rti.le.rtimn.and.j.eq.nray) rtimn=rti
@ -5636,14 +5650,13 @@ c
end do
c
if(.not.(iproj.eq.0.and.j.eq.1))
. write(nfile,111) istep,j,k,xti1,yti1,zti1,rti1,
. sqrt(psinv11)
. write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
if(iproj.eq.1) write(nfile,*) ' '
end do
c
write(nfile,*) ' '
c
write(12,99) istep,st,sqrt(abs(psinv11)),rtimn,rtimx
write(12,99) istep,st,psinv11,rtimn,rtimx
return
99 format(i5,12(1x,e16.8e3))
111 format(3i5,12(1x,e16.8e3))
@ -5956,41 +5969,32 @@ c dPdV [MW/m^3], Jcd [MA/m^2]
if(ieccd.eq.0) currt=0.0d0
currtka=currt*1.0d3
if(index_rt.eq.1) then
write(7,*)'#B0 beta alpha Icd Pa Jphi '//
.'rhotj drhotj rhotjav rhotjava drhotjava dpdvmx rhotp drhotp '//
.'rhotpav drhotpav taumn taumx stmx polpsi polchi index_rt'
write(48,*) '#beta alpha rhop rhot dPdV Jphi Jcda P% Pins'//
. ' Icdins index_rt'
end if
write(6,*)' '
write(6,*)'#beta alpha Icd Pa Jphi rhotj drhotj '//
.'rhotjav rhotjava drhotjava dpdvmx rhotp drhotp rhotpav '//
.'drhotpav taumn taumx stmx Pins_02 Pins_05 Pins_085'
write(6,99) betac,alfac,currtka,pabstot,ajmxfi,rhotjfi,drhotjfi,
. rhotjav,rhotjava,drhotjava,
. dpdvmx,rhotp,drhotp,rhotpav,drhotpav,
. taumn,taumx,stmx,pins_02,pins_05,pins_085
write(6,*)'#beta alpha Icd Pa Jphimx dPdVmx '//
.'rhotj rhotjava rhotp rhotpav '//
.'drhotjava drhotpav stmx polpsi polchi index_rt'
write(6,99) betac,alfac,currtka,pabstot,ajmxfi,dpdvmx,
. rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav,
. stmx,psipol,chipol,real(index_rt)
write(7,99) btrcen,betac,alfac,currtka,pabstot,
. ajmxfi,rhotjfi,
. drhotjfi,rhotjav,rhotjava,drhotjava,
. dpdvmx,rhotp,drhotp,rhotpav,drhotpav,
. taumn,taumx,stmx,psipol,chipol,real(index_rt)
write(7,99) currtka,pabstot,ajmxfi,dpdvmx,
. rhotjfi,rhotjava,rhotp,rhotpav,
. drhotjava,drhotpav,
. stmx,psipol,chipol,real(index_rt)
do i=1,nd
if (ipec.eq.0) then
psip=rtab(i)
psin=rtab(i)
rhop=sqrt(rtab(i))
else
psip=rtab(i)**2
psin=rtab(i)**2
rhop=rtab(i)
end if
pinsr=0.0d0
if(pabstot.gt.0) pinsr=pins(i)/pabstot
write(48,99) betac,alfac,rhop,rhotv(i),dpdv(i),ajphiv(i)
. ,ajcdv(i),pinsr,pins(i),currins(i),real(index_rt)
write(48,99) psin,rhotv(i),ajphiv(i),ajcdv(i),dpdv(i),
. currins(i),pins(i),pinsr,real(index_rt)
end do
return
@ -6246,12 +6250,6 @@ c
endif
gam=atan(sngam/csgam)*180.d0/pi
c if(ipolc.eq.0.or.ipolc.eq.1) then
c write(28,111) beta,alfa,anpl,-gam,ffo,ffx,xe2om,xe2xm
c write(29,111) beta,alfa,qqom,uuom,vvom,psiom,chiom,
c . qqxm,uuxm,vvxm,psixm,chixm
c end if
return
111 format(20(1x,e12.5))
end