diff --git a/src/gray-externals.f b/src/gray-externals.f index 4b03a35..37ae005 100644 --- a/src/gray-externals.f +++ b/src/gray-externals.f @@ -58,6 +58,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 @@ -85,7 +89,7 @@ c 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' + if(iprof.eq.0) write(*,*) 'ANALTYCAL PROFILES ne0,Te0',dens0,te0 if(ibeam.ge.1) write(*,*) 'LAUNCHER CASE : ',filenmbm end if @@ -445,14 +449,15 @@ c c subroutine read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, - . dne, zeff, qsf, powin) + . dne, zeff, qsf, powin, alphain, betain, beamid) use green_func_p, only:Setup_SpitzFunc use const_and_precisions, only : pi use beamdata, only : nrayr,nrayth,nstep implicit real*8 (a-h,o-z) - integer ijetto, mr, mz, nrho, nbnd + integer ijetto, mr, mz, nrho, nbnd, beamid real*8 r(mr), z(mz), psin(mr,mz) real*8 psiax, psibnd, rax, zax, powin + real*8 alphain, betain real*8 rbnd(nbnd), zbnd(nbnd) real*8 psijet(nrho), f(nrho), te(nrho), dne(nrho), . zeff(nrho), qsf(nrho) @@ -474,6 +479,7 @@ c common/warm/iwarm,ilarm common/ieccd/ieccd common/idst/idst + common/imx/imx c common/filesn/filenmeqq,filenmprf,filenmbm c @@ -514,6 +520,7 @@ c common/mode/sox common/angles/alpha0,beta0 common/scal/iscal + common/fghz/fghz c open(602,file='gray.data',status= 'unknown') c @@ -521,7 +528,10 @@ c alpha0, beta0 (cartesian) launching angles c fghz wave frequency (GHz) c p0mw injected power (MW) c - read(602,*) alpha0,beta0 + read(602,*) dummy, dummy !alpha0,beta0 +c angles read from values passed as arguments (in rad) + alpha0 = alphain + beta0 = betain read(602,*) fghz read(602,*) p0mw c power value overwritten with value passed as argument (in W) @@ -568,8 +578,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 c - read(602,*) iwarm,ilarm + read(602,*) iwarm,ilarm,imx c if(iwarm.gt.2) iwarm=2 c c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models @@ -675,14 +687,12 @@ c c xg=xgcn*dens19 c xgcn=1.0d-5*qe**2/(pi*me*fghz**2) -c - sox=-1.0d0 - if(iox.eq.2) sox=1.0d0 c c read data for beam from file if ibeam>0 c if(ibeam.gt.0) then - call read_beams + filenmbm='graybeam.data' + call read_beams(beamid,iox) else zrcsi=0.5d0*ak0*w0csi**2 zreta=0.5d0*ak0*w0eta**2 @@ -705,6 +715,9 @@ c phir=phiw end if end if +c + sox=-1.0d0 + if(iox.eq.2) sox=1.0d0 c write(*,'(a,2f8.3)') 'alpha0, beta0 = ',alpha0,beta0 c @@ -769,13 +782,13 @@ 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) - zlim(1)=-dst*nstep*2.d-2 + zlim(1)=(z00-dst*nstep*2)*1.d-2 zlim(2)=zlim(1) - zlim(3)=dst*nstep*2.d-2 + zlim(3)=(z00+dst*nstep*2)*1.d-2 zlim(4)=zlim(3) zlim(5)=zlim(1) ipass=abs(ipass) @@ -878,45 +891,68 @@ c 111 format(i4,20(1x,e16.8e3)) end c -c - subroutine read_beams +c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ + subroutine read_beams(beamid,iox) + use reflections, only : inside implicit real*8(a-h,o-z) +c character*255 filenmeqq,filenmprf,filenmbm - parameter(nstrmx=50) -c - dimension alphastv(nstrmx),betastv(nstrmx),cbeta(nstrmx,4) - dimension x00v(nstrmx),y00v(nstrmx),z00v(nstrmx) - dimension cx0(nstrmx,4),cy0(nstrmx,4),cz0(nstrmx,4) - dimension waist1v(nstrmx),waist2v(nstrmx) - dimension rci1v(nstrmx),rci2v(nstrmx) - dimension cwaist1(nstrmx,4),cwaist2(nstrmx,4) - dimension crci1(nstrmx,4),crci2(nstrmx,4) - dimension phi1v(nstrmx),phi2v(nstrmx) - dimension cphi1(nstrmx,4),cphi2(nstrmx,4) -c + character*20 beamname +c + integer beamid, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta, + . iopt, incheck, nxcoord, nycoord, nxest, nyest, lwrk, kwrk, + . nxwaist1, nywaist1, nxwaist2, nywaist2, nxrci1, nyrci1, nxrci2, + . nyrci2, nxphi1, nyphi1, nxphi2, nyphi2, nxx0, nyx0, nxy0, nyy0, + . nxz0, nyz0, kx, ky, ii, nfbeam, iox, npolyg, nmax, lwrk2, in + integer, DIMENSION(:), ALLOCATABLE :: iwrk +c + real*8 alphast, betast, x00, y00, z00, waist1, waist2, dvir1, + . dvir2, delta, phi1, phi, zr1, zr2, ako, w1, w2, rci1, rci2, fghz, + . fhz, wcso, weta, rcicsi, rcieta, phiw, phir, xcoord0, ycoord0, + . alpha0, beta0, fp, minx, maxx, miny, maxy, eps, dmin + real*8, DIMENSION(:), ALLOCATABLE :: x00v, y00v, z00v, alphastv, + . betastv, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, xcoord, + . ycoord, ycoord2, wrk, txwaist1, tywaist1, txwaist2, tywaist2, + . txrci1, tyrci1, txrci2, tyrci2, txphi1, typhi1, txphi2, typhi2, + . txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1, + . cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2, + . xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC, + . ypolygC, xpolygD, ypolygD + real*8 xvert(4), yvert(4) +c + parameter(kspl=3,sspl=0.01) +c common/ibeam/ibeam common/filesn/filenmeqq,filenmprf,filenmbm common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/mirr/x00,y00,z00 common/angles/alpha0,beta0 common/parwv/ak0,akinv,fhz + common/fghz/fghz c c for given alpha0 -> beta0 + beam parameters c c ibeam=1 simple astigmatic beam c ibeam=2 general astigmatic beam -c -c initial beam data are measured in mm -> transformed to cm c nfbeam=603 lbm=length(filenmbm) open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam) c - read(nfbeam,*) nisteer - do i=1,nisteer - if(ibeam.eq.1) then - read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm, - . waist1,dvir1,waist2,dvir2,delta +c####################################################################################### + if(ibeam.eq.1) then + read(nfbeam,*) nisteer + naplha = nisteer + nbeta = 1 + allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), + . waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer), + . phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer)) +c +c c==================================================================================== + do i=1,nisteer + read(nfbeam,*) alphast,betast,x00,y00,z00, waist1,dvir1, + . waist2,dvir2,delta phi1=delta phi2=delta zr1=0.5d-1*ak0*waist1**2 @@ -925,71 +961,667 @@ c w2=waist2*sqrt(1.0d0+(dvir2/zr2)**2) rci1=-dvir1/(dvir1**2+zr1**2) rci2=-dvir2/(dvir2**2+zr2**2) - else - read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm, - . w1,w2,rci1,rci2,phi1,phi2 - end if - x00v(i)=0.1d0*x00mm - y00v(i)=0.1d0*y00mm - z00v(i)=0.1d0*z00mm - alphastv(i)=alphast - betastv(i)=betast - waist1v(i)=0.1d0*w1 - rci1v(i)=1.0d1*rci1 - waist2v(i)=0.1d0*w2 - rci2v(i)=1.0d1*rci2 - phi1v(i)=phi1 - phi2v(i)=phi2 - end do -c - iopt=0 - call difcsg(alphastv,betastv,nisteer,iopt,cbeta,ier) - call difcsg(alphastv,waist1v,nisteer,iopt,cwaist1,ier) - call difcsg(alphastv,rci1v,nisteer,iopt,crci1,ier) - call difcsg(alphastv,waist2v,nisteer,iopt,cwaist2,ier) - call difcsg(alphastv,rci2v,nisteer,iopt,crci2,ier) - call difcsg(alphastv,phi1v,nisteer,iopt,cphi1,ier) - call difcsg(alphastv,phi2v,nisteer,iopt,cphi2,ier) - call difcsg(alphastv,x00v,nisteer,iopt,cx0,ier) - call difcsg(alphastv,y00v,nisteer,iopt,cy0,ier) - call difcsg(alphastv,z00v,nisteer,iopt,cz0,ier) -c - if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then - call vlocate(alphastv,nisteer,alpha0,k) - dal=alpha0-alphastv(k) - betst=fspli(cbeta,nstrmx,k,dal) - x00=fspli(cx0,nstrmx,k,dal) - y00=fspli(cy0,nstrmx,k,dal) - z00=fspli(cz0,nstrmx,k,dal) - wcsi=fspli(cwaist1,nstrmx,k,dal) - weta=fspli(cwaist2,nstrmx,k,dal) - rcicsi=fspli(crci1,nstrmx,k,dal) - rcieta=fspli(crci2,nstrmx,k,dal) - phiw=fspli(cphi1,nstrmx,k,dal) - phir=fspli(cphi2,nstrmx,k,dal) +c +c initial beam data (x00, y00, z00) are measured in mm -> transformed to cm + x00v(i)=0.1d0*x00 + y00v(i)=0.1d0*y00 + z00v(i)=0.1d0*z00 + alphastv(i)=alphast + betastv(i)=betast + waist1v(i)=0.1d0*w1 + rci1v(i)=1.0d1*rci1 + waist2v(i)=0.1d0*w2 + rci2v(i)=1.0d1*rci2 + phi1v(i)=phi1 + phi2v(i)=phi2 + end do +c c==================================================================================== +c +c free variable = alpha + fdeg = 1 +c======================================================================================= else - write(*,*) ' alpha0 outside table range !!!' - if(alpha0.ge.alphastv(nisteer)) ii=nisteer - if(alpha0.le.alphastv(1)) ii=1 - betst=betastv(ii) - x00=x00v(ii) - y00=y00v(ii) - z00=z00v(ii) - wcsi=waist1v(ii) - weta=waist2v(ii) - rcicsi=rci1v(ii) - rcieta=rci2v(ii) - phiw=phi1v(ii) - phir=phi2v(ii) +c======================================================================================= +c # of beams + read(nfbeam,*) nbeam +c +c unused beams' data + jumprow=0 +c c==================================================================================== + do i=1,beamid-1 + read(nfbeam,*) beamname, iox, fghz, nalpha, nbeta + jumprow = jumprow+nalpha*nbeta + end do +c c==================================================================================== +c +c beam of interest + read(nfbeam,*) beamname, iox, fghz, nalpha, nbeta + fhz=fghz*1.0d9 +c +c c==================================================================================== +c unused beams' data grids + do i=1,(jumprow + (nbeam - beamid)) + read(nfbeam,*) + end do +c c==================================================================================== +c +c # of elements in beam data grid + nisteer = nalpha*nbeta +c + allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), + . waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer), + . phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer)) +c +c c==================================================================================== +c beam data grid reading + do i=1,nisteer + read(nfbeam,*) alphast,betast,x00,y00,z00, + . waist1,waist2,rci1,rci2,phi1,phi2 +c +c initial beam data (x00, y00, z00) are measured in mm -> transformed to cm + x00v(i)=0.1d0*x00 + y00v(i)=0.1d0*y00 + z00v(i)=0.1d0*z00 + alphastv(i)=alphast + betastv(i)=betast + waist1v(i)=0.1d0*waist1 + rci1v(i)=1.0d1*rci1 + waist2v(i)=0.1d0*waist2 + rci2v(i)=1.0d1*rci2 + phi1v(i)=phi1 + phi2v(i)=phi2 + end do +c c==================================================================================== +c +c fdeg = 0 alpha, beta free variables +c 1 alpha free variable +c 2 beta free variable +c 3 no free variables + fdeg = 2*1/nalpha + 1/nbeta end if - beta0=betst -c +c####################################################################################### +c close(nfbeam) +c +c####################################################################################### +c no free variables + if(fdeg.eq.3) then + alpha0=alphastv(1) + beta0=betastv(1) + x00=x00v(1) + y00=y00v(1) + z00=z00v(1) + wcsi=waist1v(1) + weta=waist2v(1) + rcicsi=rci1v(1) + rcieta=rci2v(1) + phiw=phi1v(1) + phir=phi2v(1) + return + end if +c####################################################################################### +c +c +c####################################################################################### + if(fdeg.eq.2) then +c beta = independent variable +c alpha = dependent variable + xcoord = betastv + ycoord = alphastv + xcoord0 = beta0 + ycoord0 = alpha0 + kx=min(nbeta-1,kspl) +c c==================================================================================== + else +c c==================================================================================== +c alpha = independent variable +c beta = dependent/independent (fdeg = 1/0) + xcoord = alphastv + ycoord = betastv + xcoord0 = alpha0 + ycoord0 = beta0 + nxcoord = nalpha + nycoord = nbeta + kx=min(nalpha-1,kspl) + ky=min(nbeta-1,kspl) + end if +c####################################################################################### +c + iopt = 0 + incheck = 0 +c +c####################################################################################### + if(fdeg.ne.0) then + nxest = kx + nxcoord + 1 + lwrk = (nxcoord*(kx+1)+nxest*(7+3*kx)) + kwrk = nxest + allocate(cycoord(nxest), txycoord(nxest), cwaist1(nxest), + . txwaist1(nxest), cwaist2(nxest), txwaist2(nxest), + . crci1(nxest), txrci1(nxest), crci2(nxest), txrci2(nxest), + . cphi1(nxest), txphi1(nxest), cphi2(nxest), txphi2(nxest), + . cx0(nxest), txx0(nxest), cy0(nxest), txy0(nxest), + . cz0(nxest), txz0(nxest), w(nxcoord), wrk(lwrk), iwrk(kwrk)) +c + do i=1,nxcoord + w(i) = 1.d0 + end do +c +c 2D interpolation + call dierckx_curfit(iopt,nxcoord,xcoord,ycoord,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxycoord, + . txycoord,cycoord,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,waist1v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist1, + . txwaist1,cwaist1,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,waist2v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist2, + . txwaist2,cwaist2,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,rci1v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci1, + . txrci1,crci1,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,rci2v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci21, + . txrci2,crci2,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,phi1v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi1, + . txphi1,cphi1,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,phi2v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi2, + . txphi2,cphi2,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,x00v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxx0, + . txx0,cx0,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,y00v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxy0, + . txy0,cy0,fp,wrk,lwrk,iwrk,ier) +c + call dierckx_curfit(iopt,nxcoord,xcoord,z00v,w, + . xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxz0, + . txz0,cz0,fp,wrk,lwrk,iwrk,ier) +c +c check if xcoord0 is out of table range +c incheck = 1 inside / 0 outside + if(xcoord0.gt.xcoord(1).and.xcoord0.lt.xcoord(nisteer)) + . incheck=1 +c c==================================================================================== + else +c c==================================================================================== + npolyg = 2*(nxcoord+nycoord-2) + minx = minval(xcoord) + maxx = maxval(xcoord) + miny = minval(ycoord) + maxy = maxval(ycoord) + nxest = kx + 1 + sqrt(nisteer/2.) + nyest = ky + 1 + sqrt(nisteer/2.) + nmax = max(nxest,nyest) + eps = 10.**(-8) + lwrk = (nmax-2)**2*(7*nmax-2)+18*nmax+8*nisteer-19 + lwrk2 = (nmax-2)**2*(4*nmax-1)+4*nmax-2 + kwrk = nisteer+(nmax-3)*(nmax-3) + allocate(cwaist1(nxest*nyest), txwaist1(nmax), tywaist1(nmax), + . cwaist2(nxest*nyest), txwaist2(nmax), tywaist2(nmax), + . crci1(nxest*nyest), txrci1(nmax), tyrci1(nmax), + . crci2(nxest*nyest), txrci2(nmax), tyrci2(nmax), + . cphi1(nxest*nyest), txphi1(nmax), typhi1(nmax), + . cphi2(nxest*nyest), txphi2(nmax), typhi2(nmax), + . cx0(nxest*nyest), txx0(nmax), tyx0(nmax), + . cy0(nxest*nyest), txy0(nmax), tyy0(nmax), + . cz0(nxest*nyest), txz0(nmax), tyz0(nmax), + . wrk(lwrk), wrk2(lwrk2), iwrk(kwrk), + . xpolyg(npolyg), ypolyg(npolyg), w(nisteer)) +c + do i=1,nisteer + w(i) = 1.d0 + end do +c +c 3D interpolation + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,waist1v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxwaist1,txwaist1,nywaist1,tywaist1,cwaist1, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,waist2v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxwaist2,txwaist2,nywaist2,tywaist2,cwaist2, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,rci1v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxrci1,txrci1,nyrci1,tyrci1,crci1, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,rci2v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxrci2,txrci2,nyrci2,tyrci2,crci2, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,phi1v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxphi1,txphi1,nyphi1,typhi1,cphi1, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,phi2v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxphi2,txphi2,nyphi2,typhi2,cphi2, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,x00v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxx0,txx0,nyx0,tyx0,cx0, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,y00v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxy0,txy0,nyy0,tyy0,cy0, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c + call dierckx_surfit(iopt,nisteer,xcoord,ycoord,z00v,w, + . minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, + . nxz0,txz0,nyz0,tyz0,cz0, + . fp,wrk,lwrk,wrk2,lwkr2,iwrk,kwrk,ier) +c data range polygon + xpolyg(1:nxcoord) = xcoord(1:nxcoord) + ypolyg(1:nxcoord) = ycoord(1:nxcoord) +c +c c==================================================================================== + do i=1,nycoord-2 + xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord) + xpolyg(2*nxcoord+nycoord-2+i) = + . xcoord((nycoord-i-1)*nxcoord+1) + ypolyg(nxcoord+i) = ycoord((i+1)*nxcoord) + ypolyg(2*nxcoord+nycoord-2+i) = + . ycoord((nycoord-i-1)*nxcoord+1) + end do +c c==================================================================================== + do i=1,nxcoord + xpolyg(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1) + ypolyg(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1) + end do +c c==================================================================================== +c +c check if (xcoord0, ycoord0) is out of table range +c incheck = 1 inside / 0 outside + if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0)) + . incheck = 1 + deallocate(wrk,iwrk) + end if +c####################################################################################### +c +c +c####################################################################################### + if(fdeg.ne.0) then +c c==================================================================================== + if(incheck.eq.1) then + call dierckx_splev(txycoord,nxycoord,cycoord,kx,xcoord0, + . ycoord0,1,ier) +c + call dierckx_splev(txycoord,nxycoord,cwaist1,kx,xcoord0,wcsi, + . 1,ier) +c + call dierckx_splev(txwaist2,nxwaist2,cwaist2,kx,xcoord0,weta, + . 1,ier) +c + call dierckx_splev(txrci1,nxrci1,crci1,kx,xcoord0,rcicsi, + . 1,ier) +c + call dierckx_splev(txrci2,nxrci2,crci2,kx,xcoord0,rcieta, + . 1,ier) +c + call dierckx_splev(txphi1,nxphi1,cphi1,kx,xcoord0,phiw,1,ier) +c + call dierckx_splev(txphi2,nxphi2,cphi2,kx,xcoord0,phir,1,ier) +c + call dierckx_splev(txx0,nxx0,cx0,kx,xcoord0,x00,1,ier) +c + call dierckx_splev(txy0,nxy0,cy0,kx,xcoord0,y00,1,ier) +c + call dierckx_splev(txz0,nxz0,cz0,kx,xcoord0,z00,1,ier) +c c---------------------------------------------------------------------------------- + else +c c---------------------------------------------------------------------------------- + write(*,*) ' alpha0 or beta0 outside table range !!!' + if(xcoord0.ge.xcoord(nisteer)) ii=nisteer + if(xcoord0.le.xcoord(1)) ii=1 +c + xcoord0=xcoord(ii) + ycoord0=ycoord(ii) + x00=x00v(ii) + y00=y00v(ii) + z00=z00v(ii) + wcsi=waist1v(ii) + weta=waist2v(ii) + rcicsi=rci1v(ii) + rcieta=rci2v(ii) + phiw=phi1v(ii) + phir=phi2v(ii) + end if +c c==================================================================================== + else +c c==================================================================================== + if(incheck.eq.0) then + allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord), + . ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord), + . xpolygD(nycoord), ypolygD(nycoord)) +c coordinates of vertices v1,v2,v3,v4 + xvert(1) = xpolyg(1) + xvert(2) = xpolyg(nxcoord) + xvert(3) = xpolyg(nxcoord+nycoord-1) + xvert(4) = xpolyg(2*nxcoord+nycoord-2) + yvert(1) = ypolyg(1) + yvert(2) = ypolyg(nxcoord) + yvert(3) = ypolyg(nxcoord+nycoord-1) + yvert(4) = ypolyg(2*nxcoord+nycoord-2) +c coordinates of side A,B,C,D + xpolygA = xpolyg(1:nxcoord) + ypolygA = ypolyg(1:nxcoord) + xpolygB = xpolyg(nxcoord:nxcoord+nycoord-1) + ypolygB = ypolyg(nxcoord:nxcoord+nycoord-1) + xpolygC = xpolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2) + ypolygC = ypolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2) + xpolygD(1:nycoord-1) = xpolyg(2*nxcoord+nycoord-2:npolyg) + xpolygD(nycoord) = xpolyg(1) + ypolygD(1:nycoord-1) = ypolyg(2*nxcoord+nycoord-2:npolyg) + ypolygD(nycoord) = ypolyg(1) +c c---------------------------------------------------------------------------------- +c search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid +c +c | | +c (6) (5) (4) +c | | +c _ _ _ v4 _________________v3_ _ _ _ +c | C | (1)->(8) outside regions +c | | +c (7) D | | B (3) v1->v4 grid vertices +c | | +c _ _ _ _ |_________________| _ _ _ _ A-D grid sides +c v1 A v2 +c | | +c (8) (1) (2) +c | | +c + if(xcoord0.gt.xvert(1).and.xcoord0.lt.xvert(2).and. + . ycoord0.le.maxval(ypolygA)) then + in=1 + else if(ycoord0.gt.yvert(2).and.ycoord0.lt.yvert(3).and. + . xcoord0.ge.minval(xpolygB)) then + in=3 + else if(xcoord0.lt.xvert(3).and.xcoord0.gt.xvert(4).and. + . ycoord0.ge.minval(ypolygC)) then + in=5 + else if(ycoord0.lt.yvert(4).and.ycoord0.gt.yvert(1).and. + . xcoord0.le.maxval(xpolygD)) then + in=7 + else if(xcoord0.ge.xvert(2).and.ycoord0.le.yvert(2)) then + in=2 + else if(xcoord0.ge.xvert(3).and.ycoord0.ge.yvert(3)) then + in=4 + else if(xcoord0.le.xvert(4).and.ycoord0.ge.yvert(4)) then + in=6 + else if(xcoord0.le.xvert(1).and.ycoord0.le.yvert(1)) then + in=8 + endif +c c---------------------------------------------------------------------------------- +c +c c---------------------------------------------------------------------------------- +c (xcoord0,ycoord0) is set to its nearest point on (alpha, beta) grid border +c depending on the region +c 1: xcoord0 unchanged, ycoord0 moved on side A +c 3: xcoord0 moved on side B, ycoord0 unchanged +c 5: xcoord0 unchanged, ycoord0 moved on side C +c 7: xcoord0 moved on side D, ycoord0 unchanged +c 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates +c in 1,3,5,7 incheck is set back to 1 to evaluate x00,y00,z00,waist,rci,phi in +c new (xcoord0,ycoord0) +c in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the +c (xcoord0,ycoord0) vertex are used + alpha0 = xcoord0 + beta0 = ycoord0 + SELECT CASE (in) + CASE (1) + write(*,*) ' beta0 outside table range !!!' +c locate position of xcoord0 with respect to x coordinates of side A + call vlocate(xpolygA,nxcoord,xcoord0,ii) +c find corresponding y value on side A for xcoord position + call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1), + . ypolygA(ii+1),xcoord0,ycoord0) + incheck = 1 + CASE (2) + write(*,*) ' alpha0 and beta0 outside table range !!!' +c xcoord0, ycoord0 set + xcoord0 = xvert(2) + ycoord0 = yvert(2) + ii = nxcoord !indice per assegnare valori waist, rci, phi + CASE (3) + write(*,*) ' alpha0 outside table range !!!' + call vlocate(ypolygB,nycoord,ycoord0,ii) + call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1), + . xpolygB(ii+1),ycoord0,xcoord0) + incheck = 1 + CASE (4) + write(*,*) ' alpha0 and beta0 outside table range !!!' + xcoord0 = xvert(3) + ycoord0 = yvert(3) + ii = nxcoord+nycoord-1 + CASE (5) + write(*,*) ' beta0 outside table range !!!' + call vlocate(xpolygC,nxcoord,xcoord0,ii) + call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii), + . ypolygC(ii),xcoord0,ycoord0) + incheck = 1 + CASE (6) + write(*,*) ' alpha0 and beta0 outside table range !!!' + xcoord0 = xvert(4) + ycoord0 = yvert(4) + ii = 2*nxcoord+nycoord-2 + CASE (7) + write(*,*) ' alpha0 outside table range !!!' + call vlocate(ypolygD,nycoord,ycoord0,ii) + call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1), + . xpolygD(ii+1),ycoord0,xcoord0) + incheck = 1 + CASE (8) + write(*,*) ' alpha0 and beta0 outside table range !!!!' + xcoord0 = xvert(1) + ycoord0 = yvert(1) + ii = 1 + END SELECT +c c---------------------------------------------------------------------------------- +c +c print new (alpha0, beta0) + write(*,*) ' (alpha0, beta0) values changed' + write(*,10) ' old (alpha0,beta0) = (', alpha0, + . ',', beta0, ')' + write(*,10) ' new (alpha0,beta0) = (', xcoord0, + . ',', ycoord0, ')' + 10 format(a27,f10.5,a1,f10.5,a1) + end if +c c==================================================================================== +c +c c==================================================================================== + if(incheck.eq.1) then + lwrk = 2*(kx+ky+2) + kwrk = 4 + allocate(wrk(lwrk),iwrk(kwrk)) +c + call dierckx_bispev(txwaist1,nxwaist1,tywaist1,nywaist1, + . cwaist1,kx,ky,(/xcoord0/),1,(/ycoord0/),1,wcsi, + . wrk,lwrk,iwrk,kwrk,ier) +c + call dierckx_bispev(txwaist2,nxwaist2,tywaist2,nywaist2, + . cwaist2,kx,ky,(/xcoord0/),1,(/ycoord0/),1,weta, + . wrk,lwrk,iwrk,kwrk,ier) +c + call dierckx_bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1,kx, + . ky,(/xcoord0/),1,(/ycoord0/),1,rcicsi,wrk,lwrk,iwrk,kwrk,ier) +c + call dierckx_bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2,kx, + . ky,(/xcoord0/),1,(/ycoord0/),1,rcieta,wrk,lwrk,iwrk,kwrk,ier) +c + call dierckx_bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1,kx, + . ky,(/xcoord0/),1,(/ycoord0/),1,phiw,wrk,lwrk,iwrk,kwrk,ier) +c + call dierckx_bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2,kx, + . ky,(/xcoord0/),1,(/ycoord0/),1,phir,wrk,lwrk,iwrk,kwrk,ier) +c + call dierckx_bispev(txx0,nxx0,tyx0,nyx0,cx0,kx,ky, + . (/xcoord0/),1,(/ycoord0/),1,x00,wrk,lwrk,iwrk,kwrk,ier) +c + call dierckx_bispev(txy0,nxy0,tyy0,nyy0,cy0,kx,ky, + . (/xcoord0/),1,(/ycoord0/),1,y00,wrk,lwrk,iwrk,kwrk,ier) +c + call dierckx_bispev(txz0,nxz0,tyz0,nyz0,cz0,kx,ky, + . (/xcoord0/),1,(/ycoord0/),1,z00,wrk,lwrk,iwrk,kwrk,ier) +c c---------------------------------------------------------------------------------- + else +c c---------------------------------------------------------------------------------- + x00=x00v(ii) + y00=y00v(ii) + z00=z00v(ii) + wcsi=waist1v(ii) + weta=waist2v(ii) + rcicsi=rci1v(ii) + rcieta=rci2v(ii) + phiw=phi1v(ii) + phir=phi2v(ii) + end if +c c==================================================================================== + end if +c####################################################################################### +c +c +c####################################################################################### +c set correct values for alpha, beta + if(fdeg.eq.2) then + alpha0 = ycoord0 + beta0 = xcoord0 + else + alpha0 = xcoord0 + beta0 = ycoord0 + end if +c####################################################################################### c return end -c -c +c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ + subroutine read_beams_OLD +C implicit real*8(a-h,o-z) +C character*255 filenmeqq,filenmprf,filenmbm +C parameter(nstrmx=50) +Cc +C dimension alphastv(nstrmx),betastv(nstrmx),cbeta(nstrmx,4) +C dimension x00v(nstrmx),y00v(nstrmx),z00v(nstrmx) +C dimension cx0(nstrmx,4),cy0(nstrmx,4),cz0(nstrmx,4) +C dimension waist1v(nstrmx),waist2v(nstrmx) +C dimension rci1v(nstrmx),rci2v(nstrmx) +C dimension cwaist1(nstrmx,4),cwaist2(nstrmx,4) +C dimension crci1(nstrmx,4),crci2(nstrmx,4) +C dimension phi1v(nstrmx),phi2v(nstrmx) +C dimension cphi1(nstrmx,4),cphi2(nstrmx,4) +Cc +C common/ibeam/ibeam +C common/filesn/filenmeqq,filenmprf,filenmbm +C common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir +C common/mirr/x00,y00,z00 +C common/angles/alpha0,beta0 +C common/parwv/ak0,akinv,fhz +Cc +Cc for given alpha0 -> beta0 + beam parameters +Cc +Cc ibeam=1 simple astigmatic beam +Cc ibeam=2 general astigmatic beam +Cc +Cc initial beam data are measured in mm -> transformed to cm +Cc +C nfbeam=603 +C lbm=length(filenmbm) +C open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam) +Cc +C read(nfbeam,*) nisteer +C do i=1,nisteer +C if(ibeam.eq.1) then +C read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm, +C . waist1,dvir1,waist2,dvir2,delta +C phi1=delta +C phi2=delta +C zr1=0.5d-1*ak0*waist1**2 +C zr2=0.5d-1*ak0*waist2**2 +C w1=waist1*sqrt(1.0d0+(dvir1/zr1)**2) +C w2=waist2*sqrt(1.0d0+(dvir2/zr2)**2) +C rci1=-dvir1/(dvir1**2+zr1**2) +C rci2=-dvir2/(dvir2**2+zr2**2) +C else +C read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm, +C . w1,w2,rci1,rci2,phi1,phi2 +C end if +C x00v(i)=0.1d0*x00mm +C y00v(i)=0.1d0*y00mm +C z00v(i)=0.1d0*z00mm +C alphastv(i)=alphast +C betastv(i)=betast +C waist1v(i)=0.1d0*w1 +C rci1v(i)=1.0d1*rci1 +C waist2v(i)=0.1d0*w2 +C rci2v(i)=1.0d1*rci2 +C phi1v(i)=phi1 +C phi2v(i)=phi2 +C end do +Cc +C iopt=0 +C call difcsg(alphastv,betastv,nisteer,iopt,cbeta,ier) +C call difcsg(alphastv,waist1v,nisteer,iopt,cwaist1,ier) +C call difcsg(alphastv,rci1v,nisteer,iopt,crci1,ier) +C call difcsg(alphastv,waist2v,nisteer,iopt,cwaist2,ier) +C call difcsg(alphastv,rci2v,nisteer,iopt,crci2,ier) +C call difcsg(alphastv,phi1v,nisteer,iopt,cphi1,ier) +C call difcsg(alphastv,phi2v,nisteer,iopt,cphi2,ier) +C call difcsg(alphastv,x00v,nisteer,iopt,cx0,ier) +C call difcsg(alphastv,y00v,nisteer,iopt,cy0,ier) +C call difcsg(alphastv,z00v,nisteer,iopt,cz0,ier) +Cc +C if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then +C call vlocate(alphastv,nisteer,alpha0,k) +C dal=alpha0-alphastv(k) +C betst=fspli(cbeta,nstrmx,k,dal) +C x00=fspli(cx0,nstrmx,k,dal) +C y00=fspli(cy0,nstrmx,k,dal) +C z00=fspli(cz0,nstrmx,k,dal) +C wcsi=fspli(cwaist1,nstrmx,k,dal) +C weta=fspli(cwaist2,nstrmx,k,dal) +C rcicsi=fspli(crci1,nstrmx,k,dal) +C rcieta=fspli(crci2,nstrmx,k,dal) +C phiw=fspli(cphi1,nstrmx,k,dal) +C phir=fspli(cphi2,nstrmx,k,dal) +C else +C write(*,*) ' alpha0 outside table range !!!' +C if(alpha0.ge.alphastv(nisteer)) ii=nisteer +C if(alpha0.le.alphastv(1)) ii=1 +C betst=betastv(ii) +C x00=x00v(ii) +C y00=y00v(ii) +C z00=z00v(ii) +C wcsi=waist1v(ii) +C weta=waist2v(ii) +C rcicsi=rci1v(ii) +C rcieta=rci2v(ii) +C phiw=phi1v(ii) +C phir=phi2v(ii) +C end if +C beta0=betst +Cc +C close(nfbeam) +Cc +C return + end +c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ c subroutine equidata(ijetto,mr,mz,r,z,psin2d,psiaxis,psiedge, . rax,zax,nbnd,rbnd,zbnd,mrho,psijet,fpjet,qjet) @@ -1349,8 +1981,8 @@ c compute B_toroidal on axis btaxis=fpol(1)/rmaxis btrcen=btrcen*factb write(*,'(a,f8.4)') 'factb = ',factb - write(*,'(a,f8.4)') '|BT_centr|= ',abs(btrcen) - write(*,'(a,f8.4)') '|BT_axis| = ',abs(btaxis) + write(*,'(a,f8.4)') 'BT_centr= ',btrcen + write(*,'(a,f8.4)') 'BT_axis = ',btaxis c compute normalized rho_tor from eqdsk q profile call rhotor(nrho) @@ -3694,6 +4326,7 @@ c implicit real*8 (a-h,o-z) parameter(zero=0.0d0,izero=0,one=1.0d0) parameter(cvdr=pi/180.0d0) + dimension ytmp(6),yptmp(6) 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 @@ -3707,6 +4340,11 @@ c common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/psinv/psinv + common/dddens/dens,ddens + common/tete/tekev + common/nplr/anpl,anpr + common/parpl/brr,bphi,bzz,ajphi c ui=(0.0d0,1.0d0) csth=anz0c @@ -3906,6 +4544,10 @@ 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) + grad2(j,k)=gr2 dgrad2v(1,j,k)=dgr2x dgrad2v(2,j,k)=dgr2y @@ -3919,20 +4561,20 @@ 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(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . -1.0d0,zero,zero,zero,one + write(633,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(604,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 @@ -3964,12 +4606,18 @@ c implicit real*8 (a-h,o-z) parameter(zero=0.0d0,izero=0,one=1.0d0) parameter(cvdr=pi/180.0d0) + dimension ytmp(6),yptmp(6) c common/rwmax/rwmax common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/psinv/psinv + common/dddens/dens,ddens + common/tete/tekev + common/nplr/anpl,anpr + common/parpl/brr,bphi,bzz,ajphi c csth=anz0c snth=sqrt(1.0d0-csth**2) @@ -4053,6 +4701,10 @@ 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) + do iv=1,3 gri(iv,j,k)=0.0d0 dgrad2v(iv,j,k)=0.0d0 @@ -4067,21 +4719,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(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . -1.0d0,zero,zero,zero,one + write(633,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(617,99) zero,zero,zero,zero write(604,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 @@ -4109,8 +4761,14 @@ c implicit real*8 (a-h,o-z) parameter(zero=0.0d0,izero=0,one=1.0d0) parameter(cvdr=pi/180.0d0) + dimension ytmp(6),yptmp(6) c common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 + common/psinv/psinv + common/dddens/dens,ddens + common/tete/tekev + common/nplr/anpl,anpr + common/parpl/brr,bphi,bzz,ajphi do j=1,nrayr do k=1,nrayth @@ -4141,6 +4799,10 @@ 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) + do iv=1,3 gri(iv,j,k)=0.0d0 dgrad2v(iv,j,k)=0.0d0 @@ -4155,21 +4817,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(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m, - . -1.0d0,zero,zero,zero,one + write(633,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(617,99) zero,zero,zero,zero - write(604,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(604,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 @@ -4237,8 +4899,9 @@ c common/eqnn/nr,nz,nrho,npp,nintp c ip=int((nintp-1)*rpsi+1) - if(ip.eq.0) ip=1 - if(ip.eq.nintp) ip=nintp-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 @@ -4268,8 +4931,9 @@ c common/eqnn/nr,nz,nrho,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 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=fspli(cratja,nnintp,ip,dps) ratjbi=fspli(cratjb,nnintp,ip,dps) @@ -4428,7 +5092,6 @@ c ratiovgr=2.0d0*anpr/derdnm*vgm alpha=2.0d0*akim*ratiovgr if(alpha.lt.0.0d0) then ierr=94 - write(*,*) ' IERR = ', ierr,' alpha negative' end if c ithn=1 @@ -4443,7 +5106,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 @@ -4453,15 +5115,15 @@ 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 + common/imx/imx + errnpr=1.0d0 anpr2a=anprf**2 anpl2=anpl*anpl @@ -4472,7 +5134,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 @@ -4499,8 +5164,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) @@ -4520,25 +5183,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 - write(*,*) ' 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) write(*,*) ' 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) @@ -4570,6 +5246,7 @@ c end if c return +99 format(20(1x,e12.5)) end c c Fully relativistic case @@ -5968,7 +6645,9 @@ c radial coordinate of i-(i+1) interval mid point do k=1,kkk ise0=0 ii=iiv(j,k) - if (ii.lt.nstep.and.psjki(j,k,ii+1).ne.0.0d0) ii=ii+1 + if (ii.lt.nstep) then + if(psjki(j,k,ii+1).ne.0.0d0) ii=ii+1 + end if idecr=-1 is=1 do i=1,ii @@ -5987,7 +6666,7 @@ c radial coordinate of i-(i+1) interval mid point if(ise0.eq.0) then if(xxi(i).lt.rtbc) then ise0=i - isev(is)=i-1 + isev(is)=max(i-1,1) is=is+1 end if else @@ -6008,7 +6687,7 @@ c radial coordinate of i-(i+1) interval mid point end if end if end do - isev(is)=i-1 + isev(is)=max(i-1,1) c ppa1=0.0d0 cci1=0.0d0 @@ -6124,6 +6803,7 @@ c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx) rhotjfi=0.0d0 rhopfi=0.0d0 ajmxfi=0.0d0 + ajphip=0.0d0 drhotjfi=0.0d0 drhopfi=0.0d0 xps=rhopp @@ -6160,6 +6840,8 @@ c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx) else ajmxfi=0.0d0 + ajphip=0.0d0 + dpdvp=0.0d0 dpdvmx=0.0d0 rhotjfi=1.0d0 rhotjav=1.0d0 @@ -6170,6 +6852,8 @@ c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx) drhotjava=0.0d0 drhotp=0.0d0 drhotpav=0.0d0 + ratjamx=1.0d0 + ratjbmx=1.0d0 taumn=0 taumx=0 stmx=stf @@ -6267,7 +6951,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) diff --git a/src/gray.f b/src/gray.f index 2e07b0d..58fe4ee 100644 --- a/src/gray.f +++ b/src/gray.f @@ -1,15 +1,18 @@ subroutine gray(ijetto, mr, mz, mrd, r, z, psin, psiax, psibnd, - . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, - . zeff, qsf, powin, dpdv, jcd, pec, icd, ierr) + . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, + . qsf, nbeam, powin, alphin, betain, dpdv, jcd, pec, icd, ierr) c input arguments - integer ijetto, mr, mz, nbnd, nrho + integer ijetto, mr, mz, nbnd, nrho, nbeam real*8 r(mr), z(mz), psin(mrd,mz) real*8 psiax, psibnd, rax, zax real*8 rbnd(nbnd), zbnd(nbnd) real*8 psijet(nrho), f(nrho), qsf(nrho), te(nrho), dne(nrho) real*8 zeff(nrho) + real*8 powin(nbeam), alphin(nbeam), betain(nbeam) c output arguments - real*8 dpdv(nrho), jcd(nrho), powin, pec, icd + real*8 dpdv(nrho), jcd(nrho), pec, icd +c gray_main output arguments + real*8 dpdvloop(nrho), jcdloop(nrho), pecloop, icdloop integer ierr c local variables c real*8 fgray(nrho), qgray(nrho), jcdgry(nrho), icdgry @@ -43,7 +46,10 @@ c dne Electron density on JETTO radial grid [m-3] c zeff Effective nuclear charge Zeff on JETTO radial grid c qsf Safety factor on JETTO radial grid c -c powin Input ECRH power [W] +c nbeam Total number of injected beams +c powin Input ECRH power array [W] (powin(i) =< 0 means i-th beam is unused) +c alphin Beams poloidal injection angles array [rad] +c betain Beams toroidal injection angles array [rad] c c output arguments c @@ -61,23 +67,43 @@ c c JETTO coordinate system assumes toroidal angle increasing CW c in GRAY toroidal angle increases CCW --> adapt signs on input data c - do i=1,nrho - f(i)=-f(i) - qsf(i)=-qsf(i) +c do i=1,nrho +c f(i)=-f(i) +c qsf(i)=-qsf(i) +c end do +c +c set output variables to 0 +c + dpdv = 0.d0 + jcd = 0.d0 + pec = 0.d0 + icd = 0.d0 + + do j=1,nbeam +c +c call main subroutine for the j-th beam +c + if (powin(j).gt.0.0d0) then + call gray_main(ijetto, mr, mz, r, z, psin(1:mr,:), psiax, + . psibnd, rax, zax, nbnd, rbnd, zbnd, nrho, psijet, -f, te, + . dne, zeff, -qsf, j, powin(j), alphin(j), betain(j), + . dpdvloop, jcdloop, pecloop, icdloop, ierr) +c +c add contribution of j-th beam to the total +c + dpdv = dpdv + dpdvloop + jcd = jcd + jcdloop + pec = pec + pecloop + icd = icd + icdloop + end if end do c -c call main subroutine -c - call gray_main(ijetto, mr, mz, r, z, psin(1:mr,:), psiax, psibnd, - . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, - . zeff, qsf, powin, dpdv, jcd, pec, icd, ierr) -c -c adapt output data to JETTO convention on toroidal angle +c adapt output data to JETTO convention on toroidal angle c do i=1,nrho jcd(i)=-jcd(i) - f(i)=-f(i) - qsf(i)=-qsf(i) +c f(i)=-f(i) +c qsf(i)=-qsf(i) end do icd=-icd c diff --git a/src/gray_main.f90 b/src/gray_main.f90 index 9165fcf..2059770 100644 --- a/src/gray_main.f90 +++ b/src/gray_main.f90 @@ -1,16 +1,16 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, & rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, & - powin, dpdv, jcd, pec, icd, ierr) + beamid, powin, alphain, betain, dpdv, jcd, pec, icd, ierr) use const_and_precisions, only : r8 implicit none - integer, intent(in) :: ijetto, mr, mz, nrho, nbnd + integer, intent(in) :: ijetto, mr, mz, nrho, nbnd, beamid real(r8), intent(in) :: r(mr), z(mz), psin(mr,mz) real(r8), intent(in) :: psiax, psibnd, rax, zax real(r8), intent(in), dimension(nbnd) :: rbnd, zbnd real(r8), intent(in), dimension(nrho) :: psijet, f, qsf, te, dne, zeff - real(r8), intent(in) :: powin + real(r8), intent(in) :: powin, alphain, betain real(r8), intent(out), dimension(nrho) :: dpdv, jcd real(r8), intent(out) :: pec, icd integer, intent(out) :: ierr @@ -30,13 +30,14 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, & common/index_rt/index_rt common/taumnx/taumn,taumx,pabstot,currtot common/ipass/ipass - + ! read data plus initialization index_rt=1 call prfile call paraminit call read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, rax, zax, & - nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin) + nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin, alphain, betain, & + beamid) call vectinit if(iercom.eq.0) then if(igrad.eq.0) call ic_rt