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

GRAY.f bugs fixed to avoid NaN when rays do not enter into plasma

Browse Source 
GRAY.f "ad-hoc" fix in subroutine contour_psi, that MUST be checked
This commit is contained in:
Daniela Farina 2012-07-13 16:24:03 +00:00
parent 6153060056
commit b7c4630657
1 changed files with 46 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
src/gray.f
View File

@ -208,7 +208,8 @@ c
common/dsds/dst
common/index_rt/index_rt
common/ipass/ipass
c
common/rwallm/rwallm
pabstot=0.0d0
currtot=0.0d0
taumn=1d+30
@ -324,6 +325,9 @@ c print*,' '
if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1)
. istop=1
if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1
if(iovmin.eq.2.and.ipass.eq.1) istop=1
if(iovmin.eq.3) istop=1
c print ray positions for j=nray in local reference system
@ -1389,31 +1393,29 @@ c
end if
print*,' '
c
c compute B_toroidal on axis
c
btaxis=fpol(1)/rmaxis
btrcen=abs(btrcen)*factb
print'(a,f8.4)','factb = ',factb
print'(a,f8.4)','|BT_centr|= ',btrcen
print'(a,f8.4)','|BT_axis| = ',abs(btaxis)
c
c compute normalized rho_tor from eqdsk q profile
c
call rhotor(nr)
c phitedge=deltapsi*rhotsx*2*pi
c rrtor=sqrt(phitedge/abs(btrcen)/pi)
c
c compute flux surface averaged quantities
c
rup=rmaxis
rlw=rmaxis
zup=(zbmax+zmaxis)/2.0d0
zlw=(zmaxis+zbmin)/2.0d0
zup=zmaxis+(zbmax-zmaxis)/10.0d0
zlw=zmaxis-(zmaxis-zbmin)/10.0d0
call flux_average
c
c locate psi surface for q=1.5 and q=2
c
rup=rmaxis
rlw=rmaxis
zup=(zbmax+zmaxis)/2.0d0
@ -1809,7 +1811,7 @@ c
function frhotor_av(psi)
implicit real*8(a-h,o-z)
parameter(nnintp=101)
parameter(nnintp=41)
dimension rpstab(nnintp),crhotq(nnintp,4)
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
@ -2055,18 +2057,17 @@ c
if(ier.gt.0) print*,' profil =',ier
val=h+psiaxis0/psia
call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
call sort (mest, zeroc)
if (zeroc(1).gt.rwallm) then
rcn(ic)=zeroc(1)
zcn(ic)=zc
rcn(2*np+2-ic)=zeroc(2)
zcn(2*np+2-ic)=zc
iic=1
else
rcn(ic)=zeroc(2)
zcn(ic)=zc
rcn(2*np+2-ic)=zeroc(3)
zcn(2*np+2-ic)=zc
iic=2
if (zeroc(2).lt.rwallm) iic=3
end if
rcn(ic)=zeroc(iic)
zcn(ic)=zc
rcn(2*np+2-ic)=zeroc(iic+1)
zcn(2*np+2-ic)=zc
end do
if (ipr.gt.0) then
do ii=1,2*np+1
@ -2077,7 +2078,22 @@ c
end if
return
111 format(i6,12(1x,e12.5))
99 format(2i6,12(1x,e16.8e3))
99 format(12(1x,e12.5))
end
subroutine sort(n,a)
implicit none
integer n,i,j
double precision a(n),temp
do i = 2, n
j = i - 1
temp = a(i)
do while (j.ge.1.and.a(j).gt.temp)
a(j+1) = a(j)
j = j - 1
end do
a(j+1) = temp
end do
end
c
c
@ -2086,7 +2102,7 @@ c
implicit real*8 (a-h,o-z)
real*8 lam
parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=101)
parameter(nnintp=41,ncnt=100,ncntt=2*ncnt+1,nlam=41)
parameter(zero=0.0d0,one=1.0d0)
parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi))
parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
@ -2179,10 +2195,10 @@ c computation of flux surface averaged quantities
ffc(1)=fc
rhot2q=0.0d0
rup=rmaxis
rlw=rmaxis
zup=(zbmax+zmaxis)/2.0d0
zlw=(zmaxis+zbmin)/2.0d0
C rup=rmaxis
C rlw=rmaxis
C zup=(zbmax+zmaxis)/2.0d0
C zlw=(zmaxis+zbmin)/2.0d0
do jp=2,nintp
height=dble(jp-1)/dble(nintp-1)
@ -4208,7 +4224,7 @@ c
subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi,
. bmxi,bmni,fci,intp)
implicit real*8 (a-h,o-z)
parameter(nnintp=101)
parameter(nnintp=41)
dimension rpstab(nnintp),cbmx(nnintp,4),cbmn(nnintp,4)
dimension cvol(nnintp,4),crri(nnintp,4),crbav(nnintp,4)
dimension carea(nnintp,4),cfc(nnintp,4)
@ -4242,7 +4258,7 @@ c
c
subroutine ratioj(rpsi,ratj1i,ratj2i)
implicit real*8 (a-h,o-z)
parameter(nnintp=101)
parameter(nnintp=41)
dimension rpstab(nnintp),cratj1(nnintp,4),cratj2(nnintp,4)
common/pstab/rpstab
common/eqnn/nr,nz,npp,nintp
@ -5550,7 +5566,7 @@ c gg=F(u)/u with F(u) as in Cohen paper
subroutine vlambda(alam,psi,fv,dfv)
implicit real*8 (a-h,o-z)
parameter (nnintp=101,nlam=101)
parameter (nnintp=41,nlam=41)
parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+
. njest*nnintp+nlest+54)