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