From cb55af3857abb5d3488bba5ce1b1cb5a291e860f Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Sat, 23 Aug 2014 09:57:40 +0000 Subject: [PATCH 01/11] corrected psi grid leading dimension in JINTRAC interface. Committed new branch for updated wall reflection algorithm. --- src/gray.f | 504 +++++++++++++++++++++++++++++++++-------------------- 1 file changed, 311 insertions(+), 193 deletions(-) diff --git a/src/gray.f b/src/gray.f index c32421d..95d8f50 100644 --- a/src/gray.f +++ b/src/gray.f @@ -2,7 +2,6 @@ common/istop/istop common/ierr/ierr common/igrad/igrad - common/iovmin/iopmin,iowmin common/mode/sox common/p0/p0mw common/powrfl/powrfl @@ -35,7 +34,7 @@ c postprocessing currtott=currtot powtr=p0mw-pabstot - if (iowmin.eq.2.and.ipass.gt.1) then + if (ipass.gt.1) then c second pass into plasma p0mw1=p0mw igrad=0 @@ -169,14 +168,18 @@ c c subroutine after_onestep(i,istop) implicit real*8 (a-h,o-z) + complex*16 ext,eyt,extr,eytr,exin2,eyin2 parameter(jmx=31,kmx=36,nmx=8000) 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) - dimension iop(jmx,kmx),iow(jmx,kmx) - dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6) - dimension xv(3),anv(3),xvrfl(3),anvrfl(3) - + dimension iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx),istore(jmx,kmx) + 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 ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) +c common/pcjki/ppabs,ccci common/atjki/tauv,alphav common/tau1v/tau1v @@ -187,7 +190,8 @@ c common/istgr/istpr,istpl c common/iiv/iiv - common/iov/iop,iow + common/iov/iop,iow,ihcd,istore + common/refln/anwcl,jclosest common/psjki/psjki common/psival/psinv common/psinv11/psinv11 @@ -199,8 +203,8 @@ c c common/p0/p0mw common/pol0/psipol0,chipol0 - common/ipol/ipolc - common/iovmin/iopmin,iowmin + common/polcof/psipol,chipol + common/evt/ext,eyt common/densbnd/psdbnd common/yyrfl/yyrfl common/powrfl/powrfl @@ -209,9 +213,12 @@ c common/dsds/dst common/index_rt/index_rt common/ipass/ipass - common/rwallm/rwallm common/bound/zbmin,zbmax - + common/limiter/rlim,zlim,nlim +c + common/igrad/igrad + common/wrk/ywrk,ypwrk +c pabstot=0.0d0 currtot=0.0d0 taumn=1d+30 @@ -219,6 +226,7 @@ c psinv11=1.0d0 iopmin=100 iowmin=100 + iowmax=0 c do j=1,nrayr kkk=nrayth @@ -242,7 +250,7 @@ c exit zzm=xv(3)/100.d0 if(j.eq.1) rrm11=rrm - if (iwarm.gt.0.and.i.gt.1) then + if (iwarm.gt.0.and.i.gt.1.and.ihcd(j,k).gt.0) then if(psinv.ge.0.and.psinv.le.1.0d0.and. . zzm.ge.zbmin.and.zzm.le.zbmax) then call pabs_curr(i,j,k) @@ -257,87 +265,90 @@ c exit end if call print_output(i,j,k) - if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd) - . iop(j,k)=1 - if(iop(j,k).eq.1.and. - . (psinv.ge.psdbnd.or. - . (psinv.lt.1.0d0.and.(zzm.lt.zbmin.or.zzm.gt.zbmax)))) - . iop(j,k)=2 -c iov=0 initially, iov=1 first entrance in plasma, -c iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + if (mod(iop(j,k),2).eq.0 .and. 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))) + + if (ipass.gt.1 .and. index_rt.eq.1 .and. + . iowmax.gt.1 .and. istore(j,k).eq.0) then + istore(j,k)=istore(j,k)+1 + yyrfl(j,k,1:3)=xv + 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)) + iop(j,k)=iop(j,k)+1 + call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) + end if - if(index_rt.eq.1) then - if(j.eq.1) then - psinv11=psinv - if(ipolc.eq.0.and.iop(j,k).eq.1) then - call pol_limit(qqin,uuin,vvin) - ipolc=1 - qqa=cos(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr) - uua=sin(2.0d0*psipol0*cvdr)*cos(2.0d0*chipol0*cvdr) - vva=sin(2.0d0*chipol0*cvdr) - powa=0.5d0*(1.0d0+vva*vvin+uua*uuin+qqa*qqin) -c p0mw=p0mw*powa -c print*,' ' -c print*,'Coupled power fraction =',powa -c print*,' ' -c print*,'Input coupled power (MW) =',p0mw -c print*,' ' - end if - if (iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1 - . .and.ipolc.eq.1) then - call pol_limit(qqout,uuout,vvout) - ipolc=2 - call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) - strfl11=i*dst+dstvac - call pol_limit(qqin2,uuin2,vvin2) - if(irfl.gt.0) then - powrfl=0.5d0*(1.0d0+vvrfl*vvin2+uurfl*uuin2+qqrfl*qqin2) - else - powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2) + if (ipass.gt.1) then + 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)) + iow(j,k)=2 + if (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 - write(6,*) 'Reflected power fraction =',powrfl - iop(j,k)=3 - yyrfl(j,k,1)=xvrfl(1) - yyrfl(j,k,2)=xvrfl(2) - yyrfl(j,k,3)=xvrfl(3) - yyrfl(j,k,4)=anvrfl(1) - yyrfl(j,k,5)=anvrfl(2) - yyrfl(j,k,6)=anvrfl(3) + call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)), + . eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl) + istore(j,k)=istore(j,k)+1 + yyrfl(j,k,1:3)=xvrfl + yyrfl(j,k,4:6)=anvrfl tau1v(j,k)=tauv(j,k,iiv(j,k)) + ext(j,k,iop(j,k))=extr + eyt(j,k,iop(j,k))=eytr + if (j.lt.jclosest) then + jclosest=j + anwcl=anw + end if + xv=xvrfl + anv=anvrfl + rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2) + zzm=1.d-2*xv(3) + ywrk(1:3,j,k)=xv + ywrk(4:6,j,k)=anv + igrad=0 + call gwork(j,k) + if (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))) + if (index_rt.eq.1) ihcd(j,k)=0 + end if + end if end if - else - if(iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1) then - call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) - iop(j,k)=3 - yyrfl(j,k,1)=xvrfl(1) - yyrfl(j,k,2)=xvrfl(2) - yyrfl(j,k,3)=xvrfl(3) - yyrfl(j,k,4)=anvrfl(1) - yyrfl(j,k,5)=anvrfl(2) - yyrfl(j,k,6)=anvrfl(3) - tau1v(j,k)=tauv(j,k,iiv(j,k)) - 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 end do end do - if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1 + if(jclosest.le.nrayr) then + aknmin=1.0d0 + do j=1,nrayr + kkk=nrayth + if(j.eq.1) kkk=1 + do k=1,kkk + akdotn=dot_product(anvjk(j,k),anwcl) + if(akdotn.lt.aknmin) aknmin=akdotn + end do + end do + else + aknmin=-1.0d0 + end if - psimin=psjki(1,1,i) - if(nrayr.gt.1) - . psimin=min(psimin,minval(psjki(2:nrayr,1:nrayth,i))) - if(psimin.gt.1.0d0.and.rrm11.gt.rcen.and.index_rt.gt.1) - . istop=1 +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1 - if(iopmin.eq.2.and.ipass.eq.1) istop=1 - - if(iopmin.eq.3) istop=1 - c print ray positions for j=nrayr in local reference system istpr=istpr+1 @@ -353,7 +364,82 @@ c print*,'istep = ',i c if (istpl.eq.istpl0) istpl=0 istpl=istpl+1 - + + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c single pass is stopped when all the rays have crossed the plasma +c or complete absorption has occurred +c same for successive passes of multi-pass simulations (here exit +c from vessel is detected too +c first pass in multi-pass simulation is stopped when at least one +c ray has reflected and all rays are directed away from +c reflection point, or when no reflection has occurred and +c central ray re-enters the plasma + if((ipass.eq.1 .and. ((iopmin.gt.1) .or. + . (taumn.lt.1d+30.and.taumn.gt.taucr))) + . .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or. + . (taumn.lt.1d+30.and.taumn.gt.taucr)))) then + 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 +c flag second pass mode coupling as unset + powrfl=-1.0d0 + qqout=0.0d0 + uuout=0.0d0 + vvout=0.0d0 + do j=1,nrayr + kkk=nrayth + if(j.eq.1) kkk=1 + do k=1,kkk +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) + 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) +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) + call fwork(y,dery) + call pol_limit(exin2,eyin2) + call stokes(exin2,eyin2,qqin2,uuin2,vvin2) + powloop: do j1=1,nrayr + kkkk=nrayth + if(j1.eq.1) kkkk=1 + do k1=1,kkkk +c look for first ray which completed the first pass in the plasma + if (iop(j1,k1).gt.1) then +c if found, use its polarization state to compute mode coupling + call stokes(ext(j1,k1,2),eyt(j1,k1,2), + . qqout,uuout,vvout) + exit powloop + end if + end do + end do +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 + call stokes(ext(1,1,0),eyt(1,1,0),qqout,uuout,vvout) + end if + powrfl=0.5d0*(1.0d0+vvout*vvin2+ + . uuout*uuin2+qqout*qqin2) + end if + end if + end do + end do + strfl11=i*dst + write(6,*) 'Reflected power fraction =',powrfl + istop=1 + end if + return end c @@ -2605,16 +2691,18 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda parameter(jmx=31,kmx=36,nmx=8000) dimension psjki(jmx,kmx,nmx) dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) - dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) + dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),tau1v(jmx,kmx) dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),tau1v(jmx,kmx) + dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) + dimension istore(jmx,kmx),anwcl(3) parameter(tmax=5,npts=500) real*8 ttv(npts+1),extv(npts+1) + common/ttex/ttv,extv c common/warm/iwarm,ilarm common/iiv/iiv - common/iov/iop,iow + common/iov/iop,iow,ihcd,istore common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj @@ -2623,11 +2711,13 @@ c common/nray/nrayr,nrayth common/nstep/nstep common/tau1v/tau1v + common/refln/anwcl,jclosest c if(nstep.gt.nmx) nstep=nmx if(nrayr.gt.jmx) nrayr=jmx if(nrayth.gt.kmx) nrayth=kmx - + jclosest=nrayr+1 + anwcl(1:3)=0.0d0 c do i=1,nstep do k=1,nrayth @@ -2643,6 +2733,8 @@ c iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 + ihcd(j,k)=1 + istore(j,k)=0 tau1v(j,k)=0.0d0 end do end do @@ -2668,10 +2760,11 @@ 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),iop(jmx,kmx),iow(jmx,kmx) + dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) + dimension istore(jmx,kmx) common/iiv/iiv - common/iov/iop,iow + common/iov/iop,iow,ihcd,istore common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj @@ -2694,6 +2787,7 @@ c iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 + ihcd(j,k)=1 end do end do end do @@ -2710,14 +2804,12 @@ c common/istgr/istpr,istpl common/ierr/ierr common/istop/istop - common/ipol/ipolc -c + istpr=0 istpl=1 ierr=0 istep=0 istop=0 - ipolc=0 c return end @@ -3841,6 +3933,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy @@ -3862,6 +3955,7 @@ c common/gradjk/gri common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/evt/ext,eyt c ui=(0.0d0,1.0d0) csth=anz0c @@ -4060,6 +4154,8 @@ c ypwrk0(4,j,k) = dgr2x/an0/2.0d0 ypwrk0(5,j,k) = dgr2y/an0/2.0d0 ypwrk0(6,j,k) = dgr2z/an0/2.0d0 +c + call pol_limit(ext(j,k,0),eyt(j,k,0)) c grad2(j,k)=gr2 dgrad2v(1,j,k)=dgr2x @@ -4121,6 +4217,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) c common/nray/nrayr,nrayth common/rwmax/rwmax @@ -4135,6 +4232,7 @@ c common/gradjk/gri common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/evt/ext,eyt c csth=anz0c snth=sqrt(1.0d0-csth**2) @@ -4217,6 +4315,8 @@ c ypwrk0(4,j,k) = 0.0d0 ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 +c + call pol_limit(ext(j,k,0),eyt(j,k,0)) c do iv=1,3 gri(iv,j,k)=0.0d0 @@ -4277,6 +4377,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) c common/nray/nrayr,nrayth common/wrk/ywrk0,ypwrk0 @@ -4287,6 +4388,7 @@ c common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/yyrfl/yyrfl + common/evt/ext,eyt do j=1,nrayr do k=1,nrayth @@ -4316,6 +4418,8 @@ c ypwrk0(4,j,k) = 0.0d0 ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 +c + call pol_limit(ext(j,k,0),eyt(j,k,0)) c do iv=1,3 gri(iv,j,k)=0.0d0 @@ -6564,17 +6668,14 @@ c return end - subroutine pol_limit(qq,uu,vv) + subroutine pol_limit(ext,eyt) implicit none - integer*4 ipolc real*8 bv(3),anv(3) real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm real*8 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2 real*8 deno,denx,anpl2,dnl,del0 - real*8 uuom,vvom,qqom,uuxm,vvxm,qqxm,ellom,ellxm,qq,uu,vv - real*8 aaom,bbom,llmom,aaxm,bbxm,llmxm,psiom,psixm,chiom,chixm real*8 pi,beta0,alpha0,gam - real*8 sox,psipol,chipol + real*8 sox complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) c @@ -6582,11 +6683,7 @@ c common/nplr/anpl,anpr common/ygyg/yg common/bb/bv - common/angles/alpha0,beta0 common/mode/sox - common/polcof/psipol,chipol - common/ipol/ipolc - common/evt/ext,eyt c anx=anv(1) any=anv(2) @@ -6614,103 +6711,139 @@ c c exom=(ffo*csgam-ui*anpl*sngam)/sqrt(deno) eyom=(-ffo*sngam-ui*anpl*csgam)/sqrt(deno) - qqom=abs(exom)**2-abs(eyom)**2 - uuom=2.0d0*dble(exom*dconjg(eyom)) - vvom=2.0d0*dimag(exom*dconjg(eyom)) - llmom=sqrt(qqom**2+uuom**2) - aaom=sqrt((1+llmom)/2.0d0) - bbom=sqrt((1-llmom)/2.0d0) - ellom=bbom/aaom - psiom=0.5d0*atan2(uuom,qqom)*180.d0/pi - chiom=0.5d0*asin(vvom)*180.d0/pi c exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx) eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx) - qqxm=abs(exxm)**2-abs(eyxm)**2 - uuxm=2.0d0*dble(exxm*dconjg(eyxm)) - vvxm=2.0d0*dimag(exxm*dconjg(eyxm)) - llmxm=sqrt(qqxm**2+uuxm**2) - aaxm=sqrt((1+llmxm)/2.0d0) - bbxm=sqrt((1-llmxm)/2.0d0) - ellxm=bbxm/aaxm - psixm=0.5d0*atan2(uuxm,qqxm)*180.d0/pi - chixm=0.5d0*asin(vvxm)*180.d0/pi c if (sox.lt.0.0d0) then - psipol=psiom - chipol=chiom ext=exom eyt=eyom - qq=qqom - vv=vvom - uu=uuom else - psipol=psixm - chipol=chixm ext=exxm eyt=eyxm - qq=qqxm - vv=vvxm - uu=uuxm endif + gam=atan(sngam/csgam)*180.d0/pi return 111 format(20(1x,e12.5)) end - - - subroutine wall_refl(xvrfl,anvrfl,qqtr,uutr,vvtr,irfl) + subroutine stokes(ext,eyt,qq,uu,vv) implicit none - integer*4 ivac,irfl - real*8 anv(3),xv(3),xvrfl(3) - real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3),xvout(3) - real*8 uutr,vvtr,qqtr,qq,uu,vv + complex*16 ext,eyt + real*8 qq,uu,vv + qq=abs(ext)**2-abs(eyt)**2 + uu=2.0d0*dble(ext*dconjg(eyt)) + vv=2.0d0*dimag(ext*dconjg(eyt)) + end subroutine stokes + + subroutine polellipse(qq,uu,vv,psipol,chipol) + implicit none + real*8 qq,uu,vv,psipol,chipol +c real*8 llm,aa,bb,ell real*8 pi - real*8 psipol,chipol,psitr,chitr + parameter(pi=3.14159265358979d0) +c llm=sqrt(qq**2+uu**2) +c aa=sqrt((1+llm)/2.0d0) +c bb=sqrt((1-llm)/2.0d0) +c ell=bb/aa + psipol=0.5d0*atan2(uu,qq)*180.d0/pi + chipol=0.5d0*asin(vv)*180.d0/pi + end subroutine polellipse + + + logical function inside_plasma(rrm,zzm) + implicit none + real*8 rrm,zzm,psdbnd,psinv,zbmin,zbmax + integer iequil + + common/densbnd/psdbnd + common/bound/zbmin,zbmax + common/psival/psinv + common/iieq/iequil + + if(iequil.eq.1) then + call equian(rrm,zzm) + else + call equinum(rrm,zzm) + end if + + if (psinv.ge.0.0d0.and.psinv.lt.psdbnd) then + if (psinv.lt.1.0d0.and.zzm.lt.zbmin.or.zzm.gt.zbmax) then + inside_plasma=.false. + else + inside_plasma=.true. + end if + else + inside_plasma=.false. + end if + + end function inside_plasma + + + subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln, + . irfl) + 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 complex*16 ui,extr,eytr,eztr,ext,eyt complex*16 evin(3),evrfl(3) parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) + integer nbb,nlim + parameter(nbb=5000) + real*8 rlim(nbb),zlim(nbb) - common/xv/xv - common/anv/anv - common/polcof/psipol,chipol - common/evt/ext,eyt - common/wrefl/walln + common/limiter/rlim,zlim,nlim - anv=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) + anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) + rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2) + zzm=1.d-2*xv(3) c computation of reflection coordinates and normal to the wall + call inters_linewall(xv/1.d2,anv0,rlim(1:nlim),zlim(1:nlim), + . nlim,smax,walln) + smax=smax*1.d2 + xvrfl=xv+smax*anv0 irfl=1 - ivac=1 - call vacuum_rt(xv,xvout,ivac) - - if(ivac.lt.0) then - irfl=0 - xvrfl=xvout - xv=xvout - anvrfl=anv - return - end if + if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then + ! first wall interface is outside-inside + if (dot_product(walln,walln) Date: Sat, 23 Aug 2014 10:19:56 +0000 Subject: [PATCH 02/11] Syntax errors in branch new-refl fixed. --- src/gray.f | 40 +++++++++++++++++++++++----------------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/src/gray.f b/src/gray.f index 95d8f50..d4804e3 100644 --- a/src/gray.f +++ b/src/gray.f @@ -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 @@ -317,7 +322,6 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) if (index_rt.eq.1) ihcd(j,k)=0 - end if end if end if end if @@ -327,8 +331,8 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 From 8f06be1281dccda558e5068d28c1ee3e49451912 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Mon, 25 Aug 2014 13:40:44 +0000 Subject: [PATCH 03/11] fixed the computation of beam polarization at launch (used if the beam hits the wall before crossing the plasma) --- src/gray.f | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/gray.f b/src/gray.f index d4804e3..ed56d1c 100644 --- a/src/gray.f +++ b/src/gray.f @@ -3933,6 +3933,7 @@ c 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 ytmp(ndim),yptmp(ndim) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) @@ -4159,6 +4160,9 @@ c ypwrk0(5,j,k) = dgr2y/an0/2.0d0 ypwrk0(6,j,k) = dgr2z/an0/2.0d0 c + ytmp=ywrk0(:,j,k) + yptmp=ypwrk0(:,j,k) + call fwork(ytmp,yptmp) call pol_limit(ext(j,k,0),eyt(j,k,0)) c grad2(j,k)=gr2 @@ -4217,6 +4221,7 @@ c 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 ytmp(ndim),yptmp(ndim) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) @@ -4320,6 +4325,9 @@ c ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 c + ytmp=ywrk0(:,j,k) + yptmp=ypwrk0(:,j,k) + call fwork(ytmp,yptmp) call pol_limit(ext(j,k,0),eyt(j,k,0)) c do iv=1,3 @@ -4376,6 +4384,7 @@ c 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 ytmp(ndim),yptmp(ndim) dimension yyrfl(jmx,kmx,ndim) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) dimension dffiu(jmx),ddffiu(jmx) @@ -4423,6 +4432,9 @@ c ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 c + ytmp=ywrk0(:,j,k) + yptmp=ypwrk0(:,j,k) + call fwork(ytmp,yptmp) call pol_limit(ext(j,k,0),eyt(j,k,0)) c do iv=1,3 From 8e593fdb1a7ad3a76a5618c25349db3da957d7af Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Tue, 14 Oct 2014 14:14:31 +0000 Subject: [PATCH 04/11] corrected detection of turning point for rays missing inner wall --- src/gray.f | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/gray.f b/src/gray.f index ed56d1c..598c1c5 100644 --- a/src/gray.f +++ b/src/gray.f @@ -178,7 +178,7 @@ 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(3,jmx,kmx),anvjk(3,jmx,kmx) + dimension xwcl(3),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 @@ -195,7 +195,7 @@ c c common/iiv/iiv common/iov/iop,iow,ihcd,istore - common/refln/anwcl,jclosest + common/refln/anwcl,xwcl,jclosest common/psjki/psjki common/psival/psinv common/psinv11/psinv11 @@ -309,6 +309,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (j.lt.jclosest) then jclosest=j anwcl=anw + xwcl=xvrfl end if xv=xvrfl anv=anvrfl @@ -343,7 +344,11 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk - akdotn=dot_product(anvjk(:,j,k),anwcl) + anwclr=(anwcl(1)*xwcl(1)+anwcl(2)*xwcl(2)) + . /sqrt(xwcl(1)**2+xwcl(2)**2) + anvjkr=(anvjk(1,j,k)*xvjk(1,j,k)+anvjk(2,j,k)*xvjk(2,j,k)) + . /sqrt(xvjk(1,j,k)**2+xvjk(2,j,k)**2) + akdotn=anwclr*anvjkr+anwcl(3)*anvjk(3,j,k) if(akdotn.lt.aknmin) aknmin=akdotn end do end do @@ -2698,7 +2703,7 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),tau1v(jmx,kmx) dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) - dimension istore(jmx,kmx),anwcl(3) + dimension istore(jmx,kmx),anwcl(3),xwcl(3) parameter(tmax=5,npts=500) real*8 ttv(npts+1),extv(npts+1) @@ -2715,13 +2720,14 @@ c common/nray/nrayr,nrayth common/nstep/nstep common/tau1v/tau1v - common/refln/anwcl,jclosest + common/refln/anwcl,xwcl,jclosest c if(nstep.gt.nmx) nstep=nmx if(nrayr.gt.jmx) nrayr=jmx if(nrayth.gt.kmx) nrayth=kmx jclosest=nrayr+1 anwcl(1:3)=0.0d0 + xwcl(1:3)=0.0d0 c do i=1,nstep do k=1,nrayth From 3a798e9f4af0b78211557d989f1ee3bf317257c0 Mon Sep 17 00:00:00 2001 From: Daniela Farina Date: Mon, 22 Dec 2014 15:30:17 +0000 Subject: [PATCH 05/11] added ipol option and computation of polarization parameters at all steps, added case imx negative to disable convergence in dispersion --- Makefile | 3 +- src/gray.f | 190 ++++++++++++++++++++++++++++++++++++++++------------- 2 files changed, 145 insertions(+), 48 deletions(-) diff --git a/Makefile b/Makefile index 160f0ae..4b4c807 100644 --- a/Makefile +++ b/Makefile @@ -11,7 +11,8 @@ vpath %.f src # Fortran compiler name and flags FC=gfortran -FFLAGS=-O3 +#FFLAGS=-O3 +FFLAGS=-O3 -Wall -fcheck=all DIRECTIVES = -DREVISION="'$(shell svnversion src)'" diff --git a/src/gray.f b/src/gray.f index 598c1c5..75cd376 100644 --- a/src/gray.f +++ b/src/gray.f @@ -122,6 +122,10 @@ c common/taumnx/taumn,taumx,pabstot,currtot common/scal/iscal common/facttn/factt,factn + + common/pardens/dens0,aln1,aln2 + common/parqte/te0,dte0,alt1,alt2 + c c print all ray positions in local reference system c @@ -149,7 +153,7 @@ c if(iequil.eq.2) write(6,*) 'EQUILIBRIUM CASE : ',filenmeqq if(iequil.eq.1) write(6,*) 'ANALTYCAL EQUILIBRIUM' if(iprof.eq.1) write(6,*) 'PROFILE file : ',filenmprf - if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES' + if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES ne0,Te0',dens0,te0 if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm end if @@ -231,6 +235,11 @@ c iopmin=100 iowmin=100 iowmax=0 + + if(i.eq.1) then + psipol=psipol0 + chipol=chipol0 + end if c do j=1,nrayr kkk=nrayth @@ -445,7 +454,9 @@ c initial polarization (possibly reflected) end do end do strfl11=i*dst + write(6,*) ' ' write(6,*) 'Reflected power fraction =',powrfl + write(66,*) psipol0,chipol0,powrfl istop=1 end if @@ -611,9 +622,11 @@ c .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp' write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' + write(66,*) "# psipol0 chipol0 powrfl" else + c If(index_rt.eq.3) then write(4,*) ' ' write(8,*) ' ' @@ -675,7 +688,9 @@ c common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/pol0/psipol0,chipol0 + common/ipol/ipol c + common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/pardens/dens0,aln1,aln2 common/parban/b0,rr0m,zr0m,rpam common/parqq/q0,qa,alq @@ -743,7 +758,7 @@ c iwarm=1 :weakly relativistic absorption c iwarm=2 :relativistic absorption, n<1 asymptotic expansion c iwarm=3 :relativistic absorption, numerical integration c ilarm :order of larmor expansion -c + read(2,*) iwarm,ilarm c if(iwarm.gt.2) iwarm=2 c @@ -768,10 +783,12 @@ c from center of mirror and with angular spread c ipass=1/2 1 or 2 passes into plasma c iox=1/2 OM/XM c idst=0/1/2 0 integration in s, 1 integr. in ct, 2 integr. in Sr +c psipol0,chipol0 polarization angles +c ipol=0 compute mode polarization at antenna, ipol=1 use polariz angles c read(2,*) dst,nstep,istpr0,istpl0,idst read(2,*) igrad,ipass,rwallm - read(2,*) iox, psipol0,chipol0 + read(2,*) iox, psipol0,chipol0,ipol c c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004 c sspl spline parameter for psi interpolation @@ -935,7 +952,7 @@ c versus psi, rhop, rhot c set simple limiter as two cylindrical walls at rwallm and r00 nlim=5 rlim(1)=rwallm - rlim(2)=r00*1.d-2 + rlim(2)=max(r00*1.d-2,rmxm) rlim(3)=rlim(2) rlim(4)=rlim(1) rlim(5)=rlim(1) @@ -982,7 +999,7 @@ c set simple limiter as two cylindrical walls at rwallm and r00 return 900 format('# Nray_r Nray_th rwmax : ',2i5,1x,es12.5) -901 format('# igrad iwarm ilarm ieccd idst : ',6i5) +901 format('# igrad iwarm ilarm ieccd idst ipol: ',7i5) 902 format('# X0 Y0 Z0 launching point (cm) : ',3(1x,es12.5)) 903 format('# w0xi w0eta d0xi d0eta (cm) phiw (deg) : ',5(1x,es12.5)) 904 format('# GRAY revision : ',a) @@ -1992,9 +2009,10 @@ c common/crhotq/crhotq rpsi=sqrt(psi) - ip=int((nintp-1)*rpsi+1) - if(ip.eq.0) ip=1 - if(ip.eq.nintp) ip=nintp-1 +c ip=int((nintp-1)*rpsi+1) +c if(ip.eq.0) ip=1 +c if(ip.eq.nintp) ip=nintp-1 + ip=min(max(1,ip),nintp-1) dps=rpsi-rpstab(ip) frhotor_av=spli(crhotq,nintp,ip,dps) return @@ -2672,9 +2690,10 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda common/pstab/rpstab common/eqnn/nr,nz,npp,nintp common/cdadrhot/cdadrhot - ip=int((nintp-1)*rpsi+1) - if(ip.eq.0) ip=1 - if(ip.eq.nintp) ip=nintp-1 +c ip=int((nintp-1)*rpsi+1) +c if(ip.eq.0) ip=1 +c if(ip.eq.nintp) ip=nintp-1 + ip=min(max(1,ip),nintp-1) dps=rpsi-rpstab(ip) fdadrhot=spli(cdadrhot,nintp,ip,dps) return @@ -2688,8 +2707,9 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda common/eqnn/nr,nz,npp,nintp common/cdvdrhot/cdvdrhot ip=int((nintp-1)*rpsi+1) - if(ip.eq.0) ip=1 - if(ip.eq.nintp) ip=nintp-1 + if((ip.le.0).or.(ip.ge.nintp)) print*,rpsi, ip +c if(ip.eq.nintp) ip=nintp-1 + ip=min(max(1,ip),nintp-1) dps=rpsi-rpstab(ip) fdvdrhot=spli(cdvdrhot,nintp,ip,dps) return @@ -3671,8 +3691,8 @@ c call fpbisp(tr,nsr,tz,nsz,ccspl,3,3, . rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2)) psinv=(ffspl(1)-psinop)/psiant - if(psinv.lt.0.0d0) - . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim +c if(psinv.lt.0.0d0) +c . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim c nur=1 nuz=0 @@ -3950,7 +3970,6 @@ c complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy complex*16 catand external catand -c parameter(ui=(0.0d0,1.0d0)) c common/nray/nrayr,nrayth common/rwmax/rwmax @@ -3967,6 +3986,9 @@ c common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/evt/ext,eyt + common/pol0/psipol0,chipol0 + common/ipol/ipol + c ui=(0.0d0,1.0d0) csth=anz0c @@ -4169,7 +4191,26 @@ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) - call pol_limit(ext(j,k,0),eyt(j,k,0)) + + if(ipol.eq.0) then + call pol_limit(ext(j,k,0),eyt(j,k,0)) + qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 + uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) + vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) + call polellipse(qq,uu,vv,psipol0,chipol0) + else + qq=cos(2.0d0*chipol0*cvdr)*cos(2.0d0*psipol0*cvdr) + uu=cos(2.0d0*chipol0*cvdr)*sin(2.0d0*psipol0*cvdr) + vv=sin(2.0d0*chipol0*cvdr) + if(qq**2.lt.1) then + deltapol=asin(vv/sqrt(1.0d0-qq**2)) + ext(j,k,0)= sqrt((1.0d0+qq)/2) + eyt(j,k,0)= sqrt((1.0d0-qq)/2)*exp(-ui*deltapol) + else + ext(j,k,0)= 1.0d0 + eyt(j,k,0)= 0.0d0 + end if + endif c grad2(j,k)=gr2 dgrad2v(1,j,k)=dgr2x @@ -4223,9 +4264,11 @@ c ray tracing initial conditions igrad=0 c subroutine ic_rt implicit real*8 (a-h,o-z) + complex*16 ui parameter(ndim=6,ndimm=3) parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0) parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) + parameter(ui=(0.0d0,1.0d0)) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension ytmp(ndim),yptmp(ndim) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) @@ -4248,6 +4291,8 @@ c common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/evt/ext,eyt + common/pol0/psipol0,chipol0 + common/ipol/ipol c csth=anz0c snth=sqrt(1.0d0-csth**2) @@ -4334,7 +4379,32 @@ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) - call pol_limit(ext(j,k,0),eyt(j,k,0)) + + if(ipol.eq.0) then + call pol_limit(ext(j,k,0),eyt(j,k,0)) + qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 + uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) + vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) + call polellipse(qq,uu,vv,psipol0,chipol0) + else + qq=cos(2.0d0*chipol0*cvdr)*cos(2.0d0*psipol0*cvdr) + uu=cos(2.0d0*chipol0*cvdr)*sin(2.0d0*psipol0*cvdr) + vv=sin(2.0d0*chipol0*cvdr) + if(qq**2.lt.1.0d0) then +c deltapol=phix-phiy, phix =0 + deltapol=atan2(vv,uu) + ext(j,k,0)= sqrt((1.0d0+qq)/2) + eyt(j,k,0)= sqrt((1.0d0-qq)/2)*exp(-ui*deltapol) + else + if(qq.gt.0.0d0) then + ext(j,k,0)= 1.0d0 + eyt(j,k,0)= 0.0d0 + else + eyt(j,k,0)= 1.0d0 + ext(j,k,0)= 0.0d0 + end if + end if + endif c do iv=1,3 gri(iv,j,k)=0.0d0 @@ -4408,6 +4478,7 @@ c common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/yyrfl/yyrfl common/evt/ext,eyt + common/pol0/psipol0,chipol0 do j=1,nrayr do k=1,nrayth @@ -4441,7 +4512,12 @@ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) + call pol_limit(ext(j,k,0),eyt(j,k,0)) + qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 + uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) + vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) + call polellipse(qq,uu,vv,psipol0,chipol0) c do iv=1,3 gri(iv,j,k)=0.0d0 @@ -4541,9 +4617,10 @@ c common/pstab/rpstab common/eqnn/nr,nz,npp,nintp c - ip=int((nintp-1)*rpsi+1) - if(ip.eq.0) ip=1 - if(ip.eq.nintp) ip=nintp-1 +c ip=int((nintp-1)*rpsi+1) +c if(ip.eq.0) ip=1 +c if(ip.eq.nintp) ip=nintp-1 + ip=min(max(1,ip),nintp-1) c dps=rpsi-rpstab(ip) c @@ -4572,9 +4649,10 @@ c common/pstab/rpstab common/eqnn/nr,nz,npp,nintp common/cratj/cratja,cratjb,cratjpl - ip=int((nintp-1)*rpsi+1) - if(ip.eq.0) ip=1 - if(ip.eq.nintp) ip=nintp-1 +c ip=int((nintp-1)*rpsi+1) +c if(ip.eq.0) ip=1 +c if(ip.eq.nintp) ip=nintp-1 + ip=min(max(1,ip),nintp-1) dps=rpsi-rpstab(ip) ratjai=spli(cratja,nintp,ip,dps) ratjbi=spli(cratjb,nintp,ip,dps) @@ -4770,15 +4848,14 @@ c complex*16 e33,e21,e31,e32 complex*16 a13,a31,a23,a32,a33 c common/ygyg/yg + common/xgxg/xg common/nplr/anpl,anprf -c common/mode/sox common/warm/iwarm,ilarm -c common/nprw/anprr,anpri common/epolar/ex,ey,ez common/amut/amu -c + errnpr=1.0d0 anpr2a=anprf**2 anpl2=anpl*anpl @@ -4789,7 +4866,10 @@ c call diel_tens_fr(e330,epsl,lrm) end if c - do i=1,imx + imxx=abs(imx) +c loop to disable convergence if failure detected + do + do i=1,imxx c do j=1,3 do k=1,3 @@ -4816,8 +4896,6 @@ c e33=e330+anpr2a*a33 c e21=-e12 c e31=e13 c e32=-e23 -c - if(i.gt.2.and.errnpr.lt.1.0d-3) go to 999 c cc4=(e11-anpl2)*(1.0d0-a33)+(a13+anpl)*(a31+anpl) cc2=-e12*e12*(1.0d0-a33) @@ -4837,25 +4915,38 @@ c end if c anpr2=(-cc2+s*sqrt(rr))/(2.0d0*cc4) -c - if(dble(anpr2).lt.0.0d0.and.dimag(anpr2).lt.0.0d0) then - anpr2=0.0d0 - print*,' Y =',yg,' nperp2 < 0' -c ierr=99 - go to 999 - end if c errnpr=abs(1.0d0-abs(anpr2)/abs(anpr2a)) + if(i.gt.1.and.errnpr.lt.1.0d-5) exit anpr2a=anpr2 - end do -c - 999 continue - if(i.gt.imx) print*,' i>imx ',yg,errnpr,i + end do + if(i.gt.imxx .and. imxx.gt.1) then + if (imx.lt.0) then + write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,"// + ."': convergence disabled.',e12.5)") xg,yg,sqrt(abs(anpr2)),anpl + imxx=1 + else + write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,"// + ."': convergence failed.',e12.5)") xg,yg,sqrt(abs(anpr2)),anpl + exit + end if + else + exit + end if + print*,yg,imx,imxx + end do +c end loop to disable convergence + if(sqrt(dble(anpr2)).lt.0.0d0 .or. anpr2.ne.anpr2 + . .or. abs(anpr2).eq.huge(one) .or. abs(anpr2).le.tiny(one)) then + write(*,"(' X =',f7.4,' Y =',f7.4,"// + . "' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(anpr2)) + ierr=99 + anpr2=(0.0d0,0.0d0) + end if c anpr=sqrt(anpr2) anprr=dble(anpr) anpri=dimag(anpr) -99 format(20(1x,e12.5)) c ex=dcmplx(0.0d0,0.0d0) ey=dcmplx(0.0d0,0.0d0) @@ -4887,6 +4978,7 @@ c end if c return +99 format(20(1x,e12.5)) end c c Fully relativistic case @@ -6279,7 +6371,9 @@ c do k=1,kkk ise0=0 ii=iiv(j,k) - if (ii.lt.nmx.and.psjki(j,k,ii+1).ne.0.0d0) ii=ii+1 + if (ii.lt.nmx) then + if(psjki(j,k,ii+1).ne.0.0d0) ii=ii+1 + end if idecr=-1 is=1 do i=1,ii @@ -6422,20 +6516,20 @@ c of gaussian profile rhpp=frhopol(rhotpav) rhpj=frhopol(rhotjava) dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp)) - ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) - call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) - call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx, . drhotp,drhopp) + if(ieccd.ne.0) then + ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) + call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi, . drhotjfi,drhopfi) xps=rhopfi -c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx) else rhotjfi=0.0d0 rhopfi=0.0d0 ajmxfi=0.0d0 + ajphip=0.0d0 drhotjfi=0.0d0 drhopfi=0.0d0 xps=rhopp @@ -6482,6 +6576,8 @@ c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx) drhotjava=0.0d0 drhotp=0.0d0 drhotpav=0.0d0 + ratjamx=0.0d0 + ratjbmx=0.0d0 taumn=0 taumx=0 stmx=stf From 61f31fce7f513421016c172a0d880f1e7b2042a4 Mon Sep 17 00:00:00 2001 From: Daniela Farina Date: Thu, 8 Jan 2015 15:49:38 +0000 Subject: [PATCH 06/11] corrected interpolation bug and updated outputs --- src/gray.f | 98 +++++++++++++++++++++++++++++++----------------------- 1 file changed, 56 insertions(+), 42 deletions(-) diff --git a/src/gray.f b/src/gray.f index 75cd376..dcacff4 100644 --- a/src/gray.f +++ b/src/gray.f @@ -2009,7 +2009,7 @@ c common/crhotq/crhotq rpsi=sqrt(psi) -c ip=int((nintp-1)*rpsi+1) + ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 ip=min(max(1,ip),nintp-1) @@ -2049,9 +2049,9 @@ c spline interpolation of rhopol versus rhotor . trp,crp,rp,wrkp,lwrkp,iwrkp,ier) c print*,ier call splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier) - do i=1,nnr - write(54,*) rhop(i),rhot(i),rhopi(i) - end do +c do i=1,nnr +c write(54,*) rhop(i),rhot(i),rhopi(i) +c end do return end @@ -2690,7 +2690,7 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda common/pstab/rpstab common/eqnn/nr,nz,npp,nintp common/cdadrhot/cdadrhot -c ip=int((nintp-1)*rpsi+1) + ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 ip=min(max(1,ip),nintp-1) @@ -2707,8 +2707,6 @@ c if(ip.eq.nintp) ip=nintp-1 common/eqnn/nr,nz,npp,nintp common/cdvdrhot/cdvdrhot ip=int((nintp-1)*rpsi+1) - if((ip.le.0).or.(ip.ge.nintp)) print*,rpsi, ip -c if(ip.eq.nintp) ip=nintp-1 ip=min(max(1,ip),nintp-1) dps=rpsi-rpstab(ip) fdvdrhot=spli(cdvdrhot,nintp,ip,dps) @@ -3988,6 +3986,11 @@ c common/evt/ext,eyt common/pol0/psipol0,chipol0 common/ipol/ipol + common/nplr/anpl,anpr + common/psival/psinv + common/parpl/brr,bphi,bzz,ajphi + common/dens/dens,ddens + common/tete/tekev c ui=(0.0d0,1.0d0) @@ -4225,20 +4228,21 @@ c vgradi=anxt*gxt+anyt*gyt+anzt*gzt ddi=2.0d0*vgradi c + 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.nrayr) then - r0=sqrt(x0**2+y0**2) - x0m=x0/1.0d2 - y0m=y0/1.0d2 - r0m=r0/1.0d2 - z0m=z0/1.0d2 - write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . -1.0d0,zero,zero,zero,one + write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, + . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then 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 + . psinv,one,dens,tekev,brr,bphi,bzz, + . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, + . zero,zero,zero,zero,zero,zero,zero,one + ddr110=dd end if if(j.eq.nrayr.and.k.eq.1) then @@ -4246,7 +4250,7 @@ c end if end do end do -c + call pweigth c if(nrayr.gt.1) then @@ -4293,6 +4297,11 @@ c common/evt/ext,eyt common/pol0/psipol0,chipol0 common/ipol/ipol + common/nplr/anpl,anpr + common/psival/psinv + common/parpl/brr,bphi,bzz,ajphi + common/dens/dens,ddens + common/tete/tekev c csth=anz0c snth=sqrt(1.0d0-csth**2) @@ -4420,25 +4429,25 @@ c vgradi=0.0d0 ddi=2.0d0*vgradi c + 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.nrayr) then - r0=sqrt(x0**2+y0**2) - x0m=x0/1.0d2 - y0m=y0/1.0d2 - r0m=r0/1.0d2 - z0m=z0/1.0d2 - write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . -1.0d0,zero,zero,zero,one + write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, + . psinv,zero,anpl,zero,one end if 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 + . psinv,one,dens,tekev,brr,bphi,bzz, + . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, + . zero,zero,zero,zero,zero,zero,zero,one end if end do end do -c + call pweigth c if(nrayr.gt.1) then @@ -4479,6 +4488,11 @@ c common/yyrfl/yyrfl common/evt/ext,eyt common/pol0/psipol0,chipol0 + common/nplr/anpl,anpr + common/psival/psinv + common/parpl/brr,bphi,bzz,ajphi + common/dens/dens,ddens + common/tete/tekev do j=1,nrayr do k=1,nrayth @@ -4533,21 +4547,21 @@ c vgradi=0.0d0 ddi=2.0d0*vgradi c + 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.nrayr) then - r0=sqrt(x0**2+y0**2) - x0m=x0/1.0d2 - y0m=y0/1.0d2 - r0m=r0/1.0d2 - z0m=z0/1.0d2 - write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . -1.0d0,zero,zero,zero,one + write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, + . psinv,zero,anpl,zero,one end if 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 + write(4,99) strfl11,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, + . psinv,one,dens,tekev,brr,bphi,bzz, + . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, + . zero,zero,zero,zero,zero,zero,zero,one end if end do end do @@ -4617,7 +4631,7 @@ c common/pstab/rpstab common/eqnn/nr,nz,npp,nintp c -c ip=int((nintp-1)*rpsi+1) + ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 ip=min(max(1,ip),nintp-1) @@ -4649,7 +4663,7 @@ c common/pstab/rpstab common/eqnn/nr,nz,npp,nintp common/cratj/cratja,cratjb,cratjpl -c ip=int((nintp-1)*rpsi+1) + ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 ip=min(max(1,ip),nintp-1) From 94212dba99b4eb4b4402c75baa2e3d85690bc86f Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Fri, 9 Jan 2015 15:05:20 +0000 Subject: [PATCH 07/11] added few missing dummy values in case of zero absorption --- src/gray.f | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/gray.f b/src/gray.f index dcacff4..c59513b 100644 --- a/src/gray.f +++ b/src/gray.f @@ -4631,7 +4631,7 @@ c common/pstab/rpstab common/eqnn/nr,nz,npp,nintp c - ip=int((nintp-1)*rpsi+1) + ip=int((nintp-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.nintp) ip=nintp-1 ip=min(max(1,ip),nintp-1) @@ -6580,6 +6580,8 @@ c of gaussian profile else ajmxfi=0.0d0 + ajphip=0.0d0 + dpdvp=0.0d0 dpdvmx=0.0d0 rhotjfi=1.0d0 rhotjav=1.0d0 @@ -6590,8 +6592,8 @@ c of gaussian profile drhotjava=0.0d0 drhotp=0.0d0 drhotpav=0.0d0 - ratjamx=0.0d0 - ratjbmx=0.0d0 + ratjamx=1.0d0 + ratjbmx=1.0d0 taumn=0 taumx=0 stmx=stf From 608d63acfe5b28082db4c676b76a9cf70ebc4f18 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Fri, 23 Jan 2015 15:08:41 +0000 Subject: [PATCH 08/11] fixed few out of bounds checks. added imx (dispersion) read from gray.data --- Makefile | 3 +-- src/gray.f | 18 +++++++++++++----- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/Makefile b/Makefile index 4b4c807..0c94f4a 100644 --- a/Makefile +++ b/Makefile @@ -11,8 +11,7 @@ vpath %.f src # Fortran compiler name and flags FC=gfortran -#FFLAGS=-O3 -FFLAGS=-O3 -Wall -fcheck=all +FFLAGS=-Wall -fcheck=all DIRECTIVES = -DREVISION="'$(shell svnversion src)'" diff --git a/src/gray.f b/src/gray.f index c59513b..f05ecff 100644 --- a/src/gray.f +++ b/src/gray.f @@ -663,6 +663,7 @@ c common/warm/iwarm,ilarm common/ieccd/ieccd common/idst/idst + common/imx/imx c common/filesn/filenmeqq,filenmprf,filenmbm c @@ -758,8 +759,10 @@ c iwarm=1 :weakly relativistic absorption c iwarm=2 :relativistic absorption, n<1 asymptotic expansion c iwarm=3 :relativistic absorption, numerical integration c ilarm :order of larmor expansion +c imx :max n of iterations in dispersion, imx<0 uses 1st +c iteration in case of failure after |imx| iterations - read(2,*) iwarm,ilarm + read(2,*) iwarm,ilarm,imx c if(iwarm.gt.2) iwarm=2 c c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models @@ -1987,8 +1990,9 @@ c common/crhot/crhot c irt=int((nr-1)*psi+1) - if(irt.eq.0) irt=1 - if(irt.eq.nr) irt=nr-1 +c if(irt.eq.0) irt=1 +c if(irt.eq.nr) irt=nr-1 + irt=min(max(1,irt),nr-1) dps=psi-psinr(irt) frhotor_eq=spli(crhot,nr,irt,dps) return @@ -4852,7 +4856,6 @@ c c subroutine dispersion(lrm) implicit real*8(a-h,o-z) - parameter(imx=20) complex*16 cc0,cc2,cc4,rr complex*16 anpr2a,anpra complex*16 anpr,anpr2,ex,ey,ez,den @@ -4869,6 +4872,7 @@ c common/nprw/anprr,anpri common/epolar/ex,ey,ez common/amut/amu + common/imx/imx errnpr=1.0d0 anpr2a=anprf**2 @@ -6681,7 +6685,11 @@ c rte1=0.0d0 end if call locatex(yy,nd,ipk,nd,yye,ie2) - call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2) + if(ie2.gt.0.and.ie2.lt.nd) then + call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2) + else + rte2=0.0d0 + end if else ipk=2 xpk=xx(2) From 1df304f221854859e15d2541aa708f2426fc7b4b Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Wed, 18 Feb 2015 16:17:21 +0000 Subject: [PATCH 09/11] few fixes in pec plus other minor changes --- src/gray.f | 66 +++++++++++++++++++++++++++++------------------------- 1 file changed, 36 insertions(+), 30 deletions(-) diff --git a/src/gray.f b/src/gray.f index f05ecff..4d73c90 100644 --- a/src/gray.f +++ b/src/gray.f @@ -69,7 +69,12 @@ c second pass into plasma subroutine gray_integration - implicit real*8 (a-h,o-z) + implicit none + integer istep,istop,index_rt + integer nstep + real*8 st,dst,strfl11 + integer i + real*8 st0 common/ss/st common/dsds/dst @@ -137,24 +142,24 @@ c c c print final results on screen c - print*,' ' - print'(a,f9.4)','final step (s, ct, Sr) = ',st + write(*,*) + write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st if(iwarm.gt.0) then - print '(a,2e12.5)','taumn, taumx = ', taumn,taumx + write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx else - print '(a,2f9.4)','taumn, taumx = ', zero,zero + write(*,'(a,2f9.4)') 'taumn, taumx = ', zero,zero end if c - print'(a,f9.4)','Pabs_tot (MW) = ',pabstot + write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot currtka =currtot*1.0d3 - print'(a,f9.4)','I_tot (kA) = ',currtka + write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka c if (index_rt.eq.1) then - if(iequil.eq.2) write(6,*) 'EQUILIBRIUM CASE : ',filenmeqq - if(iequil.eq.1) write(6,*) 'ANALTYCAL EQUILIBRIUM' - if(iprof.eq.1) write(6,*) 'PROFILE file : ',filenmprf - if(iprof.eq.0) write(6,*) 'ANALTYCAL PROFILES ne0,Te0',dens0,te0 - if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm + if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq + if(iequil.eq.1) write(*,*) 'ANALTYCAL EQUILIBRIUM' + if(iprof.eq.1) write(*,*) 'PROFILE file : ',filenmprf + if(iprof.eq.0) write(*,*) 'ANALTYCAL PROFILES ne0,Te0',dens0,te0 + if(ibeam.ge.1) write(*,*) 'LAUNCHER CASE : ',filenmbm end if c @@ -248,7 +253,7 @@ c call gwork(j,k) c if(ierr.gt.0) then - print*,' IERR = ', ierr + write(*,*) ' IERR = ', ierr if(ierr.eq.97) then c igrad=0 c ierr=0 @@ -959,9 +964,9 @@ c set simple limiter as two cylindrical walls at rwallm and r00 rlim(3)=rlim(2) rlim(4)=rlim(1) rlim(5)=rlim(1) - zlim(1)=-dst*nmx*1.d-2 + zlim(1)=(z00-dst*nmx)*1.d-2 zlim(2)=zlim(1) - zlim(3)=dst*nmx*1.d-2 + zlim(3)=(z00+dst*nmx)*1.d-2 zlim(4)=zlim(3) zlim(5)=zlim(1) ipass=abs(ipass) @@ -1568,8 +1573,8 @@ c compute B_toroidal on axis 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) + print'(a,f8.4)','BT_centr= ',btrcen + print'(a,f8.4)','BT_axis = ',btaxis c compute normalized rho_tor from eqdsk q profile call rhotor(nr) @@ -3693,8 +3698,8 @@ c call fpbisp(tr,nsr,tz,nsz,ccspl,3,3, . rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2)) psinv=(ffspl(1)-psinop)/psiant -c if(psinv.lt.0.0d0) -c . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim + if(psinv.lt.0.0d0) + . print'(a,3e12.4)', ' psin < 0 , R, z ',psinv,rpsim,zpsim c nur=1 nuz=0 @@ -3945,8 +3950,7 @@ c fzeff=zeff else call locate(psrad,npp,ps,k) - if(k.eq.0) k=1 - if(k.eq.npp) k=npp-1 + k=max(1,min(k,npp-1)) dps=ps-psrad(k) fzeff=spli(cz,npmx,k,dps) endif @@ -6420,7 +6424,8 @@ c else if(xxi(i).gt.rtbc) exit if(xxi(i).lt.xxi(i-1)) then - isev(is)=i +! isev(is)=i !!!!!!!!!! it should be isev(is)=i-1 + isev(is)=i-1 is=is+1 idecr=-1 end if @@ -6446,12 +6451,13 @@ c ind2=nd iind=1 end if - do ind=ind1,ind2,iind - iind=ind - if (idecr.eq.-1) iind=ind-1 - rt1=rtab1(iind) + do ind=ind1,ind2,iind !!!!!!!!!! is it safe? +! iind=ind !!!!!!!!!! iind reused in the loop! + indi=ind !!!!!!!!!! iind reused in the loop! + if (idecr.eq.-1) indi=ind-1 + rt1=rtab1(indi) call locatex(xxi,iise,iise0,iise,rt1,itb1) - if(itb1.gt.iise0.and.itb1.lt.iise) then + if(itb1.ge.iise0.and.itb1.lt.iise) then call intlin(xxi(itb1),ypt(itb1),xxi(itb1+1), . ypt(itb1+1),rt1,ppa2) call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1), @@ -6534,12 +6540,12 @@ c of gaussian profile rhpp=frhopol(rhotpav) rhpj=frhopol(rhotjava) dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp)) + ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) + call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) + call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx, . drhotp,drhopp) - if(ieccd.ne.0) then - ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) - call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi, . drhotjfi,drhopfi) xps=rhopfi From 6b3a4a807d1f982ecef5a1102e2fcd2f32bb7cbc Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 21 May 2015 09:09:16 +0000 Subject: [PATCH 10/11] changed dimension of ext eyt polarization arrays to avoid out of bounds error. to be checked. --- src/gray.f | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/gray.f b/src/gray.f index 4d73c90..dc26ba8 100644 --- a/src/gray.f +++ b/src/gray.f @@ -185,7 +185,7 @@ c dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) dimension iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx),istore(jmx,kmx) dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6) - dimension ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) + dimension ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) dimension xv(3),anv(3),xvrfl(3),anvrfl(3),xvvac(3),anw(3),anwcl(3) dimension xwcl(3),xvjk(3,jmx,kmx),anvjk(3,jmx,kmx) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) @@ -3970,7 +3970,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) - complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) + complex*16 ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy @@ -4287,7 +4287,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) - complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) + complex*16 ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) c common/nray/nrayr,nrayth common/rwmax/rwmax @@ -4483,7 +4483,7 @@ c dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) - complex*16 ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) + complex*16 ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) c common/nray/nrayr,nrayth common/wrk/ywrk0,ypwrk0 From 707342a64503906a5c5abfd1b0a70399223e746f Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 21 May 2015 10:52:04 +0000 Subject: [PATCH 11/11] grayvg branch integrated in gaussfit. other branches synced