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

grayvg branch integrated in gaussfit. other branches synced

Browse Source
This commit is contained in:
Lorenzo Figini 2015-05-21 10:52:04 +00:00
commit 4cc5013d20
1 changed files with 137 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

145
src/gray.f
View File

@ -168,6 +168,7 @@ c
if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES'
if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm
call print_tox
end if
c
c compute power and current density profiles for all rays
@ -282,6 +283,7 @@ c iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma
if(index_rt.eq.1) then
if(j.eq.1) then
psinv11=psinv
call calc_tox(i)
if(ipolc.eq.0.and.iov(j,k).eq.1) then
call pol_limit(qqin,uuin,vvin)
ipolc=1
@ -370,6 +372,77 @@ c
end
c
c
c
subroutine calc_tox(i)
implicit real*8 (a-h,o-z)
c computation of O-X transmission function
c calculated only for xg=(wpl/w)**2>0.5
c values stored at maximum xg and at maximum transmission
parameter(pi=3.14159265358979d0)
dimension derxg(3),deryg(3),xv(3),anv(3)
c
common/xgxg/xg
common/ygyg/yg
common/xv/xv
common/anv/anv
common/nplr/anpl,anpr
common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11
common/nprw/anprr,anpri
common/parwv/ak0,akinv,fhz
common/dxgyg/derxg,deryg
common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx
! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx
anplopt=sqrt(yg/(yg+1.0d0))
thkb=1.8d2/pi*acos(anpl/sqrt(an2))
dlen=sqrt(derxg(1)**2+derxg(2)**2+derxg(3)**2)
if (dlen.gt.0.0d0) then
dlen=xg/dlen
tox=exp(-pi*ak0*dlen*sqrt(0.5d0*yg)*
. (2*(1.0d0+yg)*(anplopt-abs(anpl))**2+anpr**2))
else
tox=0.0d0
end if
if (xg.gt.0.5d0.and.xg.le.1.0d0) then
if (xg.gt.xgmax) then
ixmx=i
toxxmx=tox
xgmax=xg
ygxmx=yg
anplxmx=anpl
anprxmx=anpr
thkbxmx=thkb
end if
! if (tox.gt.toxmax) then
! itmx=i
! toxmax=tox
! xgtmx=xg
! ygtmx=yg
! anpltmx=anpl
! anprtmx=anpr
! thkbtmx=thkb
! end if
! write(46,98) i,tox,xg,yg,anpl,anpr,thkb,xv/1.d2,anv,derxg,deryg
end if
c
return
98 format(1x,i8,30(1x,e16.8e3))
end
c
c
c
subroutine print_tox
implicit real*8 (a-h,o-z)
common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx
! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx
write(44,98) ixmx,toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx
! write(45,98) itmx,toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx
c
return
98 format(1x,i8,30(1x,e16.8e3))
end
c
c
c
subroutine print_output(i,j,k)
implicit real*8 (a-h,o-z)
@ -518,8 +591,8 @@ 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 psin'
write(9,*) ' #istep j k xt yt zt rt psin'
write(8,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn'
write(9,*) ' #istep j k xt yt zt rt psin modvg modn vgdotn'
write(82,*) ' #istep j k xt yt zwspl zwparab srspl,etre,etim'
write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
@ -528,6 +601,7 @@ c
write(7,*)'#Icd Pa Jphimx dPdVmx '//
.'rhotj rhotjava rhotp rhotpav '//
.'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt'
write(44,*) '#i Tox X Y Npl Nperp thetakB'
write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
else
@ -852,7 +926,11 @@ c versus psi, rhop, rhot
end if
if (iequil.eq.1) call surf_anal
c
c print psi,X,Y on the R,z grid
c
call print_2d
c
if (iequil.ne.2.or.ipass.lt.0) then
c set simple limiter as two cylindrical walls at rwallm and r00
nlim=5
@ -1584,6 +1662,34 @@ c
end
c
c
c
subroutine print_2d
implicit real*8 (a-h,o-z)
parameter(nnw=501,nnh=501)
dimension rv(nnw),zv(nnh),psin(nnw,nnh)
dimension btotal(nnw,nnh)
c
common/cpsin/rv,zv,psin
common/psival/psinv
common/eqnn/nr,nz,npp,nintp
common/parpl/brr,bphi,bzz,ajphi
common/xgxg/xg
common/ygyg/yg
c
do j=1,nr
rjcm=rv(j)*1.d2
do k=1,nz
zkcm=zv(k)*1.d2
call plas_deriv(rjcm,0.d0,zkcm)
write(91,111) rv(j),zv(k),psinv,brr,bphi,bzz,xg,yg
enddo
write(91,*) ' '
enddo
return
111 format(12(1x,e12.5))
end
c
c
c
subroutine print_prof
implicit real*8 (a-h,o-z)
@ -2607,6 +2713,9 @@ c
common/ierr/ierr
common/istop/istop
common/ipol/ipolc
c used only for computation of O-X transmission function
common/toxmxx/toxxmx,xgmax,ygxmx,anplxmx,anprxmx,thkbxmx,ixmx
! common/toxmxt/toxmax,xgtmx,ygtmx,anpltmx,anprtmx,thkbtmx,itmx
c
istpr=0
istpl=1
@ -2615,6 +2724,11 @@ c
s11=0
istop=0
ipolc=0
c
ixmx=0
! itmx=0
xgmax=0.0d0
! toxmax=0.0d0
c
return
end
@ -2878,6 +2992,7 @@ c
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
dimension vgv(3),vgv11(3)
c
common/nray/nrayr,nrayth
common/dsds/dst
@ -2890,6 +3005,9 @@ c
common/gr/gr2
common/dgr/dgr2,dgr,ddgr
common/igrad/igrad
common/vgv11/vgv,vgv11
c
h=dst
hh=h*0.5d0
@ -2930,6 +3048,7 @@ c
. +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq))
end do
end do
if(j.eq.1) vgv11=vgv
end do
c
call updatepos
@ -2994,11 +3113,12 @@ c
implicit real*8 (a-h,o-z)
parameter(ndim=6)
dimension y(ndim),dery(ndim)
dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3)
dimension xv(3),anv(3),bv(3),derbv(3,3),derxg(3),deryg(3)
dimension derdxv(3),danpldxv(3),derdnv(3)
dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3)
dimension yy11(6,3),derdxv11(3,3)
dimension dery0(3)
dimension vgv(3),vgv11(3)
c
save yy11,derdxv11
@ -3019,6 +3139,7 @@ c
common/anv/anv
common/xv/xv
common/idst/idst
common/vgv11/vgv,vgv11
c
common/iplane/iplane
c
@ -5760,6 +5881,7 @@ c parameter(nxmax=2*jmx-1)
dimension txs(nxest),tys(nxest)
parameter(lwrkbsp=8*nxmax,liwrkbsp=2*nxmax)
dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)
dimension vgv(3),vgv11(3)
c
common/nray/nrayr,nrayth
common/wrk/ywrk,ypwrk
@ -5778,6 +5900,7 @@ c
common/dery0/dery0,dery0mod
common/iplane/iplane
common/gradjk/gri
common/vgv11/vgv,vgv11
c
x4m=0.0d0
x3ym=0.0d0
@ -5947,10 +6070,16 @@ c dirx=dery0(1)
c diry=dery0(2)
c dirz=dery0(3)
c
dirx=ywrk(4,j,k)
diry=ywrk(5,j,k)
dirz=ywrk(6,j,k)
dirxn=ywrk(4,j,k)
diryn=ywrk(5,j,k)
dirzn=ywrk(6,j,k)
dirn=sqrt(dirxn*dirxn+diryn*diryn+dirzn*dirzn)
dirx=vgv11(1)
diry=vgv11(2)
dirz=vgv11(3)
dir=sqrt(dirx*dirx+diry*diry+dirz*dirz)
c
vgdotn=dirxn*dirx+diryn*diry+dirzn*dirz
c
if(j.eq.1.and.k.eq.1) then
csth1=dirz/dir
@ -6024,7 +6153,7 @@ c . snth1*snps1*dx/dr+snth1*csps1*dy/dr+csth1*dz/dr
if((iproj.eq.1).or.(iproj.eq.0.and.j.eq.nrayr))
. write(nfile,111) istep,j,k,xti,yti,zti,rti,
. xdpsi,ydpsi,zdpsi
. xdpsi,ydpsi,zdpsi,psinv11,dir,dirn,vgdotn
c
x4m=x4m+xti**4
x3ym=x3ym+xti**3*yti