diff --git a/src/gray.f b/src/gray.f index 3d48d4a..bc15cd2 100644 --- a/src/gray.f +++ b/src/gray.f @@ -2,7 +2,7 @@ common/istop/istop common/ierr/ierr common/igrad/igrad - common/iovmin/iovmin + common/iovmin/iopmin,iowmin common/mode/sox common/p0/p0mw common/powrfl/powrfl @@ -19,7 +19,10 @@ c read data plus initialization call vectinit if(igrad.eq.0) call ic_rt if(igrad.gt.0) call ic_gb - if(ierr.gt.0) go to 999 + if(ierr.gt.0) then + print*,' IERR = ', ierr + stop + end if c beam/ray propagation call gray_integration @@ -32,40 +35,37 @@ c postprocessing currtott=currtot powtr=p0mw-pabstot - if (iovmin.eq.3.and.istop.eq.1.and.ipass.gt.1) then + if (iowmin.eq.2.and.ipass.gt.1) then c second pass into plasma - p0mw1=p0mw - igrad=0 + p0mw1=p0mw + igrad=0 - index_rt=2 - p0mw=p0mw1*powrfl - call prfile - call vectinit2 - call paraminit - call ic_rt2 - call gray_integration - call after_gray_integration - pabstott=pabstott+pabstot - currtott=currtott+currtot + index_rt=2 + p0mw=p0mw1*powrfl + call prfile + call vectinit2 + call paraminit + call ic_rt2 + call gray_integration + call after_gray_integration + pabstott=pabstott+pabstot + currtott=currtott+currtot - index_rt=3 - sox=-sox - p0mw=p0mw1*(1.0d0-powrfl) - call prfile - call vectinit2 - call paraminit - call ic_rt2 - call gray_integration - call after_gray_integration - pabstott=pabstott+pabstot - currtott=currtott+currtot + index_rt=3 + sox=-sox + p0mw=p0mw1*(1.0d0-powrfl) + call prfile + call vectinit2 + call paraminit + call ic_rt2 + call gray_integration + call after_gray_integration + pabstott=pabstott+pabstot + currtott=currtott+currtot end if print*,' ' write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3 -999 continue - print*,' IERR = ', ierr - stop end @@ -173,7 +173,8 @@ c 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 iov(jmx,kmx),tau1v(jmx,kmx),yyrfl(jmx,kmx,6) + dimension iop(jmx,kmx),iow(jmx,kmx) + dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6) dimension xv(3),anv(3),xvrfl(3),anvrfl(3) common/pcjki/ppabs,ccci @@ -186,7 +187,7 @@ c common/istgr/istpr,istpl c common/iiv/iiv - common/iov/iov + common/iov/iop,iow common/psjki/psjki common/psival/psinv common/psinv11/psinv11 @@ -199,7 +200,7 @@ c common/p0/p0mw common/pol0/psipol0,chipol0 common/ipol/ipolc - common/iovmin/iovmin + common/iovmin/iopmin,iowmin common/densbnd/psdbnd common/yyrfl/yyrfl common/powrfl/powrfl @@ -216,7 +217,8 @@ c taumn=1d+30 taumx=-1d+30 psinv11=1.0d0 - iovmin=100 + iopmin=100 + iowmin=100 c do j=1,nrayr kkk=nrayth @@ -256,18 +258,18 @@ c ierr=0 call print_output(i,j,k) if(i.gt.1.and.psinv.ge.0.and.psinv.lt.psdbnd) - . iov(j,k)=1 - if(iov(j,k).eq.1.and. + . 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)))) - . iov(j,k)=2 + . 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 if(index_rt.eq.1) then if(j.eq.1) then psinv11=psinv - if(ipolc.eq.0.and.iov(j,k).eq.1) then + 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) @@ -281,7 +283,7 @@ c print*,' ' c print*,'Input coupled power (MW) =',p0mw c print*,' ' end if - if (iov(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1 + 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 @@ -294,7 +296,7 @@ c print*,' ' powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2) end if write(6,*) 'Reflected power fraction =',powrfl - iov(j,k)=3 + iop(j,k)=3 yyrfl(j,k,1)=xvrfl(1) yyrfl(j,k,2)=xvrfl(2) yyrfl(j,k,3)=xvrfl(3) @@ -304,9 +306,9 @@ c print*,' ' tau1v(j,k)=tauv(j,k,iiv(j,k)) end if else - if(iov(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1) then + if(iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1) then call wall_refl(xvrfl,anvrfl,qqrfl,uurfl,vvrfl,irfl) - iov(j,k)=3 + iop(j,k)=3 yyrfl(j,k,1)=xvrfl(1) yyrfl(j,k,2)=xvrfl(2) yyrfl(j,k,3)=xvrfl(3) @@ -318,7 +320,7 @@ c print*,' ' end if end if - if(iov(j,k).lt.iovmin) iovmin=iov(j,k) + if(iop(j,k).lt.iopmin) iopmin=iop(j,k) end do end do @@ -332,9 +334,9 @@ c print*,' ' . istop=1 if(rrm11.lt.rwallm.and.ipass.eq.1) istop=1 - if(iovmin.eq.2.and.ipass.eq.1) istop=1 + if(iopmin.eq.2.and.ipass.eq.1) istop=1 - if(iovmin.eq.3) istop=1 + if(iopmin.eq.3) istop=1 c print ray positions for j=nrayr in local reference system @@ -2483,14 +2485,14 @@ 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),iov(jmx,kmx),tau1v(jmx,kmx) + dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),tau1v(jmx,kmx) 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/iov + common/iov/iop,iow common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj @@ -2517,7 +2519,8 @@ c ccci(j,k,i)=0.0d0 currj(j,k,i)=0.0d0 iiv(j,k)=1 - iov(j,k)=0 + iop(j,k)=0 + iow(j,k)=0 tau1v(j,k)=0.0d0 end do end do @@ -2543,10 +2546,10 @@ 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),iov(jmx,kmx) + dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx) common/iiv/iiv - common/iov/iov + common/iov/iop,iow common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj @@ -2567,7 +2570,8 @@ c ccci(j,k,i)=0.0d0 currj(j,k,i)=0.0d0 iiv(j,k)=1 - iov(j,k)=0 + iop(j,k)=0 + iow(j,k)=0 end do end do end do