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

Syntax errors in branch new-refl fixed.

Browse Source
This commit is contained in:
Lorenzo Figini 2014-08-23 10:19:56 +00:00
parent cb55af3857
commit d5a7ec1f80
1 changed files with 23 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
src/gray.f
View File

@ -167,9 +167,10 @@ c
c
c
subroutine after_onestep(i,istop)
use reflections, only : inside
implicit real*8 (a-h,o-z)
complex*16 ext,eyt,extr,eytr,exin2,eyin2
parameter(jmx=31,kmx=36,nmx=8000)
parameter(jmx=31,kmx=36,nmx=8000,nbb=5000)
parameter(taucr=12.0d0,pi=3.14159265358979d0,cvdr=pi/180.0d0)
dimension psjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
@ -177,8 +178,11 @@ c
dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6)
dimension ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3)
dimension xv(3),anv(3),xvrfl(3),anvrfl(3),xvvac(3),anw(3),anwcl(3)
dimension xvjk(jmx,kmx),anvjk(jmx,kmx)
dimension xvjk(3,jmx,kmx),anvjk(3,jmx,kmx)
dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
real*8 rlim(nbb),zlim(nbb)
logical inside_plasma
external inside_plasma
c
common/pcjki/ppabs,ccci
common/atjki/tauv,alphav
@ -278,7 +282,8 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
yyrfl(j,k,4:6)=anv
ihcd(j,k)=0
end if
else if (mod(iop(j,k),2).eq.1.and..not.inside_plasma(rrm,zzm))
else if (mod(iop(j,k),2).eq.1.and.
. .not.inside_plasma(rrm,zzm)) then
iop(j,k)=iop(j,k)+1
call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k)))
end if
@ -287,7 +292,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if (iow(j,k).eq.0 .and. inside(rlim,zlim,nlim,rrm,zzm)) then
iow(j,k)=1
else if (iow(j,k).eq.1 .and.
. .not.inside(rlim,zlim,nlim,rrm,zzm))
. .not.inside(rlim,zlim,nlim,rrm,zzm)) then
iow(j,k)=2
if (inside_plasma(rrm,zzm)) then
iop(j,k)=iop(j,k)+1
@ -320,15 +325,14 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
end if
end if
end if
end if
if(index_rt.eq.1 .and. j.eq.1) psinv11=psinv
if(iop(j,k).lt.iopmin) iopmin=iop(j,k)
if(iow(j,k).lt.iowmin) iowmin=iow(j,k)
if(iow(j,k).gt.iowmax) iowmax=iow(j,k)
xvjk(j,k)=xv
anvjk(j,k)=anv
xvjk(:,j,k)=xv
anvjk(:,j,k)=anv
end do
end do
@ -339,7 +343,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
kkk=nrayth
if(j.eq.1) kkk=1
do k=1,kkk
akdotn=dot_product(anvjk(j,k),anwcl)
akdotn=dot_product(anvjk(:,j,k),anwcl)
if(akdotn.lt.aknmin) aknmin=akdotn
end do
end do
@ -383,7 +387,7 @@ c central ray re-enters the plasma
istop=1
else if(ipass.gt.1 .and. index_rt.eq.1 .and.
. ((iowmin.gt.1 .and. aknmin.gt.0) .or.
. (iowmax.le.1 .and. iop(1,1).gt.2)) then
. (iowmax.le.1 .and. iop(1,1).gt.2))) then
c flag second pass mode coupling as unset
powrfl=-1.0d0
qqout=0.0d0
@ -396,18 +400,18 @@ c flag second pass mode coupling as unset
c store missing initial conditions for the second pass
if (istore(j,k).eq.0) then
istore(j,k)=istore(j,k)+1
yyrfl(j,k,1:3)=xvjk(j,k)
yyrfl(j,k,4:6)=anvjk(j,k)
yyrfl(j,k,1:3)=xvjk(:,j,k)
yyrfl(j,k,4:6)=anvjk(:,j,k)
tau1v(j,k)=tauv(j,k,iiv(j,k))
end if
c determine mode coupling at the plasma boundary
if (powrfl.lt.0.0d0) then
call vacuum_rt(xvjk(j,k),anvjk(j,k),xvvac,ivac)
call vacuum_rt(xvjk(:,j,k),anvjk(:,j,k),xvvac,ivac)
c look for first ray hitting the plasma, starting from the central
c and evaluate polarization
if (ivac.eq.1) then
y(1:3)=xvjk(j,k)
y(4:6)=anvjk(j,k)
y(1:3)=xvjk(:,j,k)
y(4:6)=anvjk(:,j,k)
call fwork(y,dery)
call pol_limit(exin2,eyin2)
call stokes(exin2,eyin2,qqin2,uuin2,vvin2)
@ -423,7 +427,7 @@ c if found, use its polarization state to compute mode coupling
exit powloop
end if
end do
end do
end do powloop
c if no ray completed a first pass in the plasma, use central ray
c initial polarization (possibly reflected)
if (qqout.le.0.0d0) then
@ -6784,11 +6788,12 @@ c ell=bb/aa
subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,
. irfl)
use reflections, only : inters_linewall,inside
implicit none
integer*4 irfl
real*8 anv(3),anv0(3),xv(3),xvrfl(3)
real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3)
real*8 pi,smax
real*8 pi,smax,rrm,zzm
complex*16 ui,extr,eytr,eztr,ext,eyt
complex*16 evin(3),evrfl(3)
parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
@ -6873,7 +6878,8 @@ c wave vector and electric field after reflection in lab frame
real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6)
real*8 xv0(3),anv0(3)
real*8 rlim(nbb),zlim(nbb)
logical plfound
logical plfound,inside_plasma
external inside_plasma
common/limiter/rlim,zlim,nlim
common/dsds/dst