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

corrected a bug in fwork for dbgr and in ic_gb for ypwrk0

Browse Source
This commit is contained in:
Daniela Farina 2015-11-04 10:21:50 +00:00
parent 75b77b01c1
commit 9eb901015e
2 changed files with 62 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile
View File

@ -11,7 +11,7 @@ vpath %.f src
# Fortran compiler name and flags
FC=gfortran
FFLAGS=-Wall -fcheck=all
FFLAGS=-O3 #-Wall -g -fcheck=all
DIRECTIVES = -DREVISION="'$(shell svnversion src)'"

88
src/gray.f
View File

@ -565,13 +565,13 @@ c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1
. psi,rhot,dens11,tekev,
. bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100,
. alpha11,tau11,pt11,ajcd11,dids11,
. dble(nhn),dble(iohkawa),dble(index_rt)
. dble(nhn),dble(iohkawa),dble(index_rt),ddr
c call polarcold(exf,eyif,ezf,elf,etf)
end if
c central ray only end
if(k.eq.1.and.j.eq.nrayr) write(17,99) st,ddr11,ddr,ddi
if(k.eq.nrayth.and.j.eq.nrayr) write(17,99) st,ddr,ddi
c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
if(istpl.eq.istpl0) then
@ -598,10 +598,10 @@ c
If(index_rt.eq.1) then
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'
.'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt Dr_11'
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(17,*) ' #sst Dr_Nr1 Di_Nr1'
write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
write(12,*) ' #i sst psi w1 w2'
write(7,*)'#Icd Pa Jphip dPdVp '//
@ -654,6 +654,7 @@ c
common/filesn/filenmeqq,filenmprf,filenmbm
c
common/nray/nrayr,nrayth
common/jray1/jray1
common/rwmax/rwmax
common/dsds/dst
common/igrad/igrad
@ -708,10 +709,16 @@ c
c nrayr number of rays in radial direction
c nrayth number of rays in angular direction
c rwmax normalized maximum radius of beam power
c rwmax=1 -> last ray at radius = waist
c
c nray=jray1 -> ray at radius = waist
c
read(2,*) nrayr,nrayth,rwmax
if(nrayr.eq.1) nrayth=1
if(nrayr.eq.1) then
nrayth= 1
jray1 = 1
else
jray1 = 1+nint((nrayr-1)/rwmax)
rwmax = dble(nrayr-1)/dble(jray1-1)
end if
c
c x00,y00,z00 coordinates of launching point
c
@ -3879,7 +3886,7 @@ c
bdotgr=bdotgr+bv(iv)*dgr(iv)
dbgr(iv)=0.0d0
do jv=1,3
dbgr(iv)=dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv)
dbgr(iv)=dbgr(iv)+dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv)
end do
end do
c
@ -4614,6 +4621,7 @@ c
common/parpl/brr,bphi,bzz,ajphi
common/dens/dens,ddens
common/tete/tekev
common/idst/idst
c
ui=(0.0d0,1.0d0)
@ -4806,13 +4814,25 @@ c
ywrk0(4,j,k)=anx0
ywrk0(5,j,k)=any0
ywrk0(6,j,k)=anz0
select case(idst)
case(1)
! integration variable: c*t
denom=1.0d0
case(2)
! integration variable: Sr
denom=an20
case default ! idst=0
! integration variable: s
denom=an0
end select
c
ypwrk0(1,j,k) = anx0/an0
ypwrk0(2,j,k) = any0/an0
ypwrk0(3,j,k) = anz0/an0
ypwrk0(4,j,k) = dgr2x/an0/2.0d0
ypwrk0(5,j,k) = dgr2y/an0/2.0d0
ypwrk0(6,j,k) = dgr2z/an0/2.0d0
ypwrk0(1,j,k) = anx0/denom
ypwrk0(2,j,k) = any0/denom
ypwrk0(3,j,k) = anz0/denom
ypwrk0(4,j,k) = dgr2x/(2.0d0*denom)
ypwrk0(5,j,k) = dgr2y/(2.0d0*denom)
ypwrk0(6,j,k) = dgr2z/(2.0d0*denom)
c
ytmp=ywrk0(:,j,k)
yptmp=ypwrk0(:,j,k)
@ -6547,12 +6567,18 @@ c
duu=abs(uu1-uu2)
c
if(duu.gt.1.d-6) then
call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ier,
call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ierp,
. liw,lw,last,iw,w)
if (ier.gt.0) ierr=90
if (ierp.gt.0) ierr=90
end if
c
rdu2t=cst2*anpl2+ygn2-1.0d0
nevalj0=0
nevalj1=0
nevalj2=0
lastj0=0
lastj1=0
lastj2=0
c
if (rdu2t.lt.0.0d0.or.cst2.eq.0.0d0) then
c
@ -6561,7 +6587,7 @@ c
iokhawa=0
if(duu.gt.1.d-4) then
call dqags(fcur,uu1,uu2,epsa,epsr,
. resj,ej,neval,ier,liw,lw,last,iw,w)
. resj,ej,nevalj0,ier,liw,lw,lastj0,iw,w)
if (ier.gt.0) ierr=91
end if
else
@ -6583,7 +6609,7 @@ c
duu=abs(uu1-uu2)
if(duu.gt.1.d-6) then
call dqags(fcur,uu1,uu2,epsa,epsr,
. resj1,ej1,neval,ier,liw,lw,last,iw,w)
. resj1,ej1,nevalj1,ier,liw,lw,lastj1,iw,w)
if (ier.gt.0) then
if (abs(resj1).lt.1.0d-10) then
resj1=0.0d0
@ -6604,7 +6630,7 @@ c
duu=abs(uu1-uu2)
if(duu.gt.1.d-6) then
call dqags(fcur,uu1,uu2,epsa,epsr,
. resj2,ej2,neval,ier,liw,lw,last,iw,w)
. resj2,ej2,nevalj2,ier,liw,lw,lastj2,iw,w)
if (ier.gt.0) then
if(ier.ne.2) ierr=93
end if
@ -6613,8 +6639,11 @@ c
resj=resj1+resj2
end if
end if
write(90,1) yg,nhn,ierp,last,neval,lastj0,nevalj0,
. lastj1,nevalj1,lastj2,nevalj2
return
1 format(e12.5,20i5)
end
c
c computation of integral for power density, integrand function fpp
@ -6820,6 +6849,7 @@ c gg=F(u)/u with F(u) as in Cohen paper
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
c
common/nray/nrayr,nrayth
common/jray1/jray1
common/wrk/ywrk,ypwrk
common/psinv11/psinv11
common/istep/istep
@ -6829,8 +6859,12 @@ c
rtimx=-1.d-30
c
jd=1
if(iproj.eq.0) jd=nrayr-1
do j=1,nrayr,jd
jf=nrayr
if(iproj.eq.0) then
jd=jray1-1
jf=jray1
end if
do j=1,jf,jd
kkk=nrayth
if(j.eq.1) kkk=1
do k=1,kkk
@ -6873,8 +6907,8 @@ C end if
if(.not.(iproj.eq.0.and.j.eq.1))
. write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
c
if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti
if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti
if(rti.ge.rtimx.and.j.eq.jray1) rtimx=rti
if(rti.le.rtimn.and.j.eq.jray1) rtimn=rti
c
end do
c
@ -7039,9 +7073,9 @@ c
ind2=nd
iind=1
end if
do ind=ind1,ind2,iind !!!!!!!!!! is it safe?
! iind=ind !!!!!!!!!! iind reused in the loop!
indi=ind !!!!!!!!!! iind reused in the loop!
do ind=ind1,ind2,iind
indi=ind
if (idecr.eq.-1) indi=ind-1
rt1=rtab1(indi)
call locatex(xxi,iise,iise0,iise,rt1,itb1)
@ -7056,7 +7090,7 @@ c
ajphiv(ind)=ajphiv(ind)+didst
ppa1=ppa2
cci1=cci2
end if
end if
end do
end do
end do