implicit real*8 (a-h,o-z) common/istop/istop common/ierr/ierr common/igrad/igrad common/iovmin/iopmin,iowmin common/mode/sox common/p0/p0mw common/powrfl/powrfl common/index_rt/index_rt common/taumnx/taumn,taumx,pabstot,currtot common/ipass/ipass c read data plus initialization index_rt=1 call prfile call paraminit call read_data call vectinit if(igrad.eq.0) call ic_rt if(igrad.gt.0) call ic_gb if(ierr.gt.0) then print*,' IERR = ', ierr stop end if c beam/ray propagation call gray_integration c postprocessing call after_gray_integration pabstott=pabstot currtott=currtot powtr=p0mw-pabstot if (iowmin.eq.2.and.ipass.gt.1) then c second pass into plasma 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=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 end subroutine gray_integration implicit real*8 (a-h,o-z) common/ss/st common/dsds/dst common/istep/istep common/nstep/nstep common/istop/istop common/strfl11/strfl11 common/index_rt/index_rt c ray integration: begin st0=0.0d0 if(index_rt.gt.1) st0=strfl11 do i=1,nstep istep=i st=i*dst+st0 c advance one step call rkint4 c calculations after one step: call after_onestep(istep,istop) if(istop.eq.1) exit c end do c ray integration: end return end subroutine after_gray_integration implicit real*8 (a-h,o-z) parameter(zero=0.0d0) character*24 filenmeqq,filenmprf,filenmbm c common/ss/st common/ibeam/ibeam common/warm/iwarm,ilarm common/filesn/filenmeqq,filenmprf,filenmbm common/nray/nrayr,nrayth common/iieq/iequil common/iipr/iprof common/index_rt/index_rt c common/p0/p0mw common/factb/factb common/taumnx/taumn,taumx,pabstot,currtot common/scal/iscal common/facttn/factt,factn c c print all ray positions in local reference system c if(nrayr.gt.1) then iproj=1 nfilp=9 call projxyzt(iproj,nfilp) end if c c print final results on screen c print*,' ' print'(a,f9.4)','final step (s, ct, Sr) = ',st if(iwarm.gt.0) then print '(a,2e12.5)','taumn, taumx = ', taumn,taumx else print '(a,2f9.4)','taumn, taumx = ', zero,zero end if c print'(a,f9.4)','Pabs_tot (MW) = ',pabstot currtka =currtot*1.0d3 print'(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' if(ibeam.ge.1) write(6,*) 'LAUNCHER CASE : ',filenmbm end if c c compute power and current density profiles for all rays c pabs=pabstot currt=currtot call pec(pabs,currt) c return 99 format(20(1x,e16.8e3)) end c c c subroutine after_onestep(i,istop) implicit real*8 (a-h,o-z) 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) common/pcjki/ppabs,ccci common/atjki/tauv,alphav common/tau1v/tau1v c common/warm/iwarm,ilarm common/nray/nrayr,nrayth common/ist/istpr0,istpl0 common/istgr/istpr,istpl c common/iiv/iiv common/iov/iop,iow common/psjki/psjki common/psival/psinv common/psinv11/psinv11 common/ierr/ierr common/taumnx/taumn,taumx,pabstot,currtot common/xv/xv common/anv/anv common/cent/btrcen,rcen c common/p0/p0mw common/pol0/psipol0,chipol0 common/ipol/ipolc common/iovmin/iopmin,iowmin common/densbnd/psdbnd common/yyrfl/yyrfl common/powrfl/powrfl common/dstvac/dstvac common/strfl11/strfl11 common/dsds/dst common/index_rt/index_rt common/ipass/ipass common/rwallm/rwallm common/bound/zbmin,zbmax pabstot=0.0d0 currtot=0.0d0 taumn=1d+30 taumx=-1d+30 psinv11=1.0d0 iopmin=100 iowmin=100 c do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk call gwork(j,k) c if(ierr.gt.0) then ! print*,' IERR = ', ierr if(ierr.eq.97) then c igrad=0 c ierr=0 else istop=1 exit end if end if psjki(j,k,i)=psinv rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0 zzm=xv(3)/100.d0 if(j.eq.1) rrm11=rrm if (iwarm.gt.0.and.i.gt.1) 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) iiv(j,k)=i else if(iiv(j,k).gt.1) tauv(j,k,i)=tauv(j,k,i-1) end if if(tauv(j,k,i).le.taumn) taumn=tauv(j,k,i) if(tauv(j,k,i).ge.taumx) taumx=tauv(j,k,i) pabstot=pabstot+ppabs(j,k,iiv(j,k)) currtot=currtot+ccci(j,k,iiv(j,k)) 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 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) 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) tau1v(j,k)=tauv(j,k,iiv(j,k)) 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(iop(j,k).lt.iopmin) iopmin=iop(j,k) end do end do if(iwarm.gt.0.and.taumn.lt.1d+30.and.taumn.gt.taucr) istop=1 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 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 if (istpr.eq.istpr0) then c print*,'istep = ',i if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) end if istpr=0 end if c if (istpl.eq.istpl0) istpl=0 istpl=istpl+1 return end c c c subroutine print_output(i,j,k) implicit real*8 (a-h,o-z) parameter(ndim=6,jmx=31,kmx=36,nmx=8000) parameter(taucr=12.0d0) parameter(pi=3.14159265358979d0) dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) dimension psjki(jmx,kmx,nmx) 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 q(jmx),tau1v(jmx,kmx) complex*16 ex,ey,ez c common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj common/pcjki/ppabs,ccci common/dcjki/didst common/nhn/nhn common/iokh/iohkawa common/p0/p0mw common/tau1v/tau1v common/qw/q common/index_rt/index_rt common/ss/st common/nray/nrayr,nrayth common/istgr/istpr,istpl common/ist/istpr0,istpl0 common/iieq/iequil c common/parpl/brr,bphi,bzz,ajphi common/btot/btot common/xgxg/xg common/ygyg/yg common/dens/dens,ddens common/tete/tekev common/absor/alpha,effjcd,akim,tau0 common/densbnd/psdbnd common/epolar/ex,ey,ez c common/nplr/anpl,anpr common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nprw/anprr,anpri c common/wrk/ywrk,ypwrk c x=ywrk(1,j,k) y=ywrk(2,j,k) z=ywrk(3,j,k) rr=sqrt(x*x+y*y) c anx=ywrk(4,j,k) any=ywrk(5,j,k) anz=ywrk(6,j,k) anr=(anx*x+any*y)/rr anphi=(any*x-anx*y)/rr c cnphi=(any*x-anx*y) c rrm=rr*1.0d-2 zzm=z*1.0d-2 stm=st*1.0d-2 xxm=x*1.0d-2 yym=y*1.0d-2 c if(index_rt.gt.1) then taujk=tauv(j,k,i)+tau1v(j,k) else taujk=tauv(j,k,i) end if c central ray only begin if(j.eq.1) then phi=acos(x/sqrt(y*y+x*x)) if(y.lt.0.0d0) phi=-phi phideg=phi*180.0d0/pi psi=psjki(j,k,i) rhot=1.0d0 bbr=0.0d0 bbz=0.0d0 bpol=0.0d0 rhot=1.0d0 dens11=0.0d0 if(psi.ge.0.0d0) then if(iequil.eq.2) then if (psi.le.1.0d0) rhot=frhotor(psi) bbr=brr bbz=bzz bpol=sqrt(bbr**2+bbz**2) else rhot=sqrt(psi) end if else tekev=0.0d0 akim=0.0d0 end if ddr11=ddr c cutoff=xg-(1-yg)*(1-anpl**2) c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1 if(psi.le.psdbnd.and.psi.ge.0.0d0) dens11=dens dids11=didst(j,k,i)*1.0d2/(p0mw*q(j)) dpdv11=pdjki(j,k,i)/(p0mw*q(j)) ajcd11=currj(j,k,i)/(p0mw*q(j)) alpha11=alphav(j,k,i) tau11=taujk pt11=exp(-taujk) write(4,99) stm,rrm,zzm,phideg, . psi,rhot,dens11,tekev, . bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100, . alpha11,tau11,pt11,ajcd11,dids11, . dble(nhn),dble(iohkawa),dble(index_rt) c call polarcold(exf,eyif,ezf,elf,etf) end if c central ray only end if(k.eq.1.and.j.eq.nrayr) write(17,99) st,ddr11,ddr,ddi c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3 if(istpl.eq.istpl0) then if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then write(33,111) i,j,k,stm,xxm,yym,rrm,zzm, . psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt) end if c if(k.eq.nrayth) write(33,*) ' ' end if c return 99 format(30(1x,e16.8e3)) 111 format(3i5,16(1x,e16.8e3)) end c c c subroutine prfile implicit none integer*4 index_rt common/index_rt/index_rt If(index_rt.eq.1) then write(4,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '// .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' write(8,*) ' #istep j k xt yt zt rt psin' write(9,*) ' #istep j k xt yt zt rt psin' write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' write(12,*) ' #i sst psi w1 w2' write(7,*)'#Icd Pa Jphimx dPdVmx '// .'rhotj rhotjava rhotp rhotpav '// .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' else c If(index_rt.eq.3) then write(4,*) ' ' write(8,*) ' ' write(9,*) ' ' write(17,*) ' ' write(12,*) ' ' write(48,*) ' ' c end if end if return end c c c subroutine read_data use green_func_p, only:Setup_SpitzFunc implicit real*8 (a-h,o-z) real*8 me character*8 wdat character*10 wtim character*24 filenmeqq,filenmprf,filenmbm parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) parameter(nmx=8000,nbb=1000) real*8 rlim(nbb),zlim(nbb) c common/xgcn/xgcn common/ipec/ipec,nnd common/nstep/nstep common/ibeam/ibeam common/ist/istpr0,istpl0 common/warm/iwarm,ilarm common/ieccd/ieccd common/idst/idst c common/filesn/filenmeqq,filenmprf,filenmbm c common/nray/nrayr,nrayth common/rwmax/rwmax common/dsds/dst common/igrad/igrad common/ipass/ipass common/rwallm/rwallm common/limiter/rlim,zlim,nlim common/iieq/iequil common/icocos/icocos common/ixp/ixp common/ipsn/ipsinorm common/sspl/sspl common/iipr/iprof common/factb/factb common/facttn/factt,factn common/cent/btrcen,rcen c common/parwv/ak0,akinv,fhz common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/pol0/psipol0,chipol0 c common/pardens/dens0,aln1,aln2 common/parban/b0,rr0m,zr0m,rpam common/parqq/q0,qa,alq common/parqte/te0,dte0,alt1,alt2 common/zz/Zeff c common/parbres/bres common/densbnd/psdbnd common/nfile/neqdsk,nprof common/sgnib/sgnbphi,sgniphi common/p0/p0mw c common/mode/sox common/angles/alpha0,beta0 common/scal/iscal c open(2,file='gray.data',status= 'unknown') c c alpha0, beta0 (cartesian) launching angles c fghz wave frequency (GHz) c p0mw injected power (MW) c read(2,*) alpha0,beta0 read(2,*) fghz read(2,*) p0mw c c nrayr number of rays in radial direction c nrayth number of rays in angular direction c rwmax normalized maximum radius of beam power c rwmax=1 -> last ray at radius = waist c read(2,*) nrayr,nrayth,rwmax if(nrayr.eq.1) nrayth=1 c c x00,y00,z00 coordinates of launching point c read(2,*) x00,y00,z00 c c beams parameters in local reference system c w0 -> waist, d0 -> waist distance from launching point c phiw angle of beam ellipse c read(2,*) w0csi,w0eta,d0csi,d0eta,phiw c c ibeam=0 :read data for beam as above c ibeam=1 :read data from file simple astigmatic beam c ibeam=2 :read data from file general astigmatic beam read(2,*) ibeam read(2,*) filenmbm c c iequil=0 :vacuum c iequil=1 :analytical equilibrium c iequil=2 :read eqdsk c ixp=0,-1,+1 : no X point , bottom/up X point c read(2,*) iequil,ixp c c iprof=0 :analytical density and temp. profiles c iprof>0 :numerical density and temp. profiles c read(2,*) iprof c c iwarm=0 :no absorption and cd 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 c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models c read(2,*) ieccd if (ieccd.ge.10) call Setup_SpitzFunc c c ipec=0/1 :pec profiles grid in psi/rhop c nnd :number of grid steps for pec profiles +1 c read(2,*) ipec,nnd c c dst integration step c nstep maximum number of integration steps c istpr0 projection step = dsdt*istprj c istpl0 plot step = dsdt*istpl c igrad=0 optical ray-tracing, initial conditions as for beam c igrad=1 quasi-optical ray-tracing c igrad=-1 ray-tracing, init. condit. 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 read(2,*) dst,nstep,istpr0,istpl0,idst read(2,*) igrad,ipass,rwallm read(2,*) iox, psipol0,chipol0 c c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004 c sspl spline parameter for psi interpolation c read(2,*) filenmeqq read(2,*) ipsinorm,sspl c c factb factor for magnetic field (only for numerical equil) c scaling adopted: beta=const, qpsi=const, nustar=const c iscal ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling c factT factn factor for Te&ne scaling c read(2,*) factb, iscal,factt,factn if (factb.eq.0.0d0) factb=1.0d0 c c psbnd value of psi ( > 1 ) of density boundary c read(2,*) filenmprf read(2,*) psdbnd if(iprof.eq.0) psdbnd=1.0d0 c c signum of toroidal B and I c icocos index for equilibrium from COCOS - O. Sauter Feb 2012 read(2,*) sgnbphi,sgniphi,icocos c c read parameters for analytical density and temperature c profiles if iprof=0 c if(iprof.eq.0) then read(2,*) dens0,aln1,aln2 read(2,*) te0,dte0,alt1,alt2 else read(2,*) dummy,dummy,dummy read(2,*) dummy,dummy,dummy,dummy end if read(2,*) zeff c c read parameters for analytical simple cylindical equilibrium c if iequil=0,1 if(iequil.lt.2) then read(2,*) rr0,zr0,rpa read(2,*) b0 read(2,*) q0,qa,alq rr0m=rr0/1.0d2 zr0m=zr0/1.0d2 rpam=rpa/1.0d2 else read(2,*) dummy,dummy,dummy read(2,*) dummy read(2,*) dummy,dummy,dummy end if c close(unit=2) c if(nrayr.eq.1) igrad=0 if (nrayr.lt.5) then igrad=0 print*,' nrayr < 5 ! => OPTICAL CASE ONLY' print*,' ' end if c fhz=fghz*1.0d9 ak0=2.0d9*pi*fghz/vc akinv=1.0d0/ak0 c bresg=2.0d0*pi*fhz*me*vc/qe bres=bresg*1.0d-4 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 else zrcsi=0.5d0*ak0*w0csi**2 zreta=0.5d0*ak0*w0eta**2 if(igrad.gt.0) then wcsi=w0csi*sqrt(1.0d0+(d0csi/zrcsi)**2) weta=w0eta*sqrt(1.0d0+(d0eta/zreta)**2) rcicsi=-d0csi/(d0csi**2+zrcsi**2) rcieta=-d0eta/(d0eta**2+zreta**2) phir=phiw else d0eta=d0csi wcsi=w0csi*abs(d0csi/zrcsi) weta=w0eta*abs(d0eta/zreta) rcicsi=w0csi/zrcsi rcieta=w0eta/zreta if(d0csi.gt.0) then rcicsi=-rcicsi rcieta=-rcieta end if phir=phiw end if end if c print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0 c r00=sqrt(x00**2+y00**2) phi0=acos(x00/r00) if(y00.lt.0) phi0=-phi0 print'(a,4f8.3)','x00, y00, R00, z00 = ',x00,y00,r00,z00 print*,' ' c c anx0c=(anr0c*x00-anphi0c*y00)/r00 c any0c=(anr0c*y00+anphi0c*x00)/r00 c c anr0c=(anx0c*x00+any0c*y00)/r00 c anphi0c=(any0c*x00-anx0c*y00)/r00 c c angles alpha0, beta0 in a local reference system as proposed by Gribov et al c c anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0) c anphi0c=sin(cvdr*beta0) c anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0) anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0) anphi0c=sin(cvdr*beta0) anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0) anx0c=(anr0c*x00-anphi0c*y00)/r00 any0c=(anr0c*y00+anphi0c*x00)/r00 c c read data for Te , ne , Zeff from file if iprof>0 c if (iprof.eq.1) then nprof=98 lprf=length(filenmprf) open(file=filenmprf(1:lprf)//'.prf', . status= 'unknown',unit=nprof) call profdata close(nprof) end if c c read equilibrium data from file if iequil=2 c if (iequil.eq.2) then neqdsk=99 leqq=length(filenmeqq) open(file=filenmeqq(1:leqq)//'.eqdsk', . status= 'unknown', unit=neqdsk) call equidata close(neqdsk) c print density, temperature, safecty factor, toroidal current dens c versus psi, rhop, rhot call print_prof end if if (iequil.eq.1) call surf_anal if (iequil.ne.2.or.ipass.lt.0) then c set simple limiter as two cylindrical walls at rwallm and r00 nlim=5 rlim(1)=rwallm rlim(2)=r00*1.d-2 rlim(3)=rlim(2) rlim(4)=rlim(1) rlim(5)=rlim(1) zlim(1)=-dst*nmx*1.d-2 zlim(2)=zlim(1) zlim(3)=dst*nmx*1.d-2 zlim(4)=zlim(3) zlim(5)=zlim(1) ipass=abs(ipass) end if nfil=78 open(nfil,file='headers.txt',status= 'unknown') call date_and_time(wdat,wtim) write(nfil,916) wdat(1:4),wdat(5:6),wdat(7:8), . wtim(1:2),wtim(3:4),wtim(5:6) write(nfil,904) REVISION if (iequil.eq.2) then write(nfil,905) filenmeqq else if (iequil.eq.1) write(nfil,912) rr0,zr0,rpa,B0,q0,qa,alq if (iequil.eq.0) write(nfil,'("# VACUUM CASE")') end if if (iprof.eq.1) then write(nfil,907) filenmprf else write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeff end if write(nfil,911) fghz,p0mw,iox write(nfil,902) x00,y00,z00 write(nfil,908) alpha0,beta0 if(ibeam.ge.1) write(nfil,909) filenmbm if(ibeam.eq.0) write(nfil,903) w0csi,w0eta,d0csi,d0eta,phiw write(nfil,900) nrayr, nrayth, rwmax write(nfil,901) igrad,iwarm,ilarm,ieccd,idst write(nfil,915) dst,nstep write(nfil,914) sgnbphi,sgniphi,icocos write(nfil,906) factb,factt,factn,iscal write(nfil,910) sspl,psdbnd,nnd,ipec,ipsinorm write(nfil,*) ' ' close(nfil) return 900 format('# Nray_r Nray_th rwmax : ',2i5,1x,es12.5) 901 format('# igrad iwarm ilarm ieccd idst : ',6i5) 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) 905 format('# EQUILIBRIUM file : ',a24) 906 format('# fact_B fact_T fact_n iscal : ',(3(1x,es12.5),i5)) 907 format('# PROFILES file : ',a24) 908 format('# alpha0 beta0 launch angles (deg) CYL : ',2(1x,es12.5)) 909 format('# LAUNCHER file : ',a24) 910 format('# sspl psdbnd nd ipec ipsinorm : ',2(1x,es12.5),3i5) 911 format('# fghz P0 IOX : ',2(1x,es12.5),i5) 912 format('# AN. EQUILIBRIUM R0 z0 a B0 q0 qa alq : ',7(1x,es12.5)) 913 format('# AN. PROFILES ne0 aln1 aln2 Te0 Tea alt1 alt2 Zeff : ' . ,8(1x,es12.5)) 914 format('# sgnB_phi sgnI_phi icocos : ',2(1x,es12.5),i5) 915 format('# dst nstep : ',1x,es12.5,i5) 916 format('# Date : ',a4,2('/',a2),1x,a2,2(':',a2)) 99 format(20(1x,e16.8e3)) end c c c subroutine surf_anal implicit real*8(a-h,o-z) parameter(pi=3.14159265358979d0) common/parban/b0,rr0m,zr0m,rpam common/parbres/bres c c print circular magnetic surfaces iequil=1 c write(71,*) '#i psin R z' dal=2.0d-2*pi drn=0.1d0 do i=1,10 rni=i*drn do k=1,101 drrm=rpam*rni*cos((k-1)*dal) dzzm=rpam*rni*sin((k-1)*dal) rrm=rr0m+drrm zzm=zr0m+dzzm write(71,111) i,rni,rrm,zzm end do write(71,*) ' ' write(71,*) ' ' end do c c print resonant B iequil=1 c write(70,*)'#i Btot R z' rres=b0*rr0m/bres zmx=zr0m+rpam zmn=zr0m-rpam do i=1,21 zzres=zmn+(i-1)*(zmx-zmn)/2.0d1 write(70,111) i,bres,rres,zzres end do return 111 format(i4,20(1x,e16.8e3)) end c c subroutine read_beams implicit real*8(a-h,o-z) character*24 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 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 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=97 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 phi1=delta phi2=delta zr1=0.5d-1*ak0*waist1**2 zr2=0.5d-1*ak0*waist2**2 w1=waist1*sqrt(1.0d0+(dvir1/zr1)**2) 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 difcs(alphastv,betastv,nisteer,iopt,cbeta,ier) call difcs(alphastv,waist1v,nisteer,iopt,cwaist1,ier) call difcs(alphastv,rci1v,nisteer,iopt,crci1,ier) call difcs(alphastv,waist2v,nisteer,iopt,cwaist2,ier) call difcs(alphastv,rci2v,nisteer,iopt,crci2,ier) call difcs(alphastv,phi1v,nisteer,iopt,cphi1,ier) call difcs(alphastv,phi2v,nisteer,iopt,cphi2,ier) call difcs(alphastv,x00v,nisteer,iopt,cx0,ier) call difcs(alphastv,y00v,nisteer,iopt,cy0,ier) call difcs(alphastv,z00v,nisteer,iopt,cz0,ier) c if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then call locate(alphastv,nisteer,alpha0,k) dal=alpha0-alphastv(k) betst=spli(cbeta,nisteer,k,dal) x00=spli(cx0,nisteer,k,dal) y00=spli(cy0,nisteer,k,dal) z00=spli(cz0,nisteer,k,dal) wcsi=spli(cwaist1,nisteer,k,dal) weta=spli(cwaist2,nisteer,k,dal) rcicsi=spli(crci1,nisteer,k,dal) rcieta=spli(crci2,nisteer,k,dal) phiw=spli(cphi1,nisteer,k,dal) phir=spli(cphi2,nisteer,k,dal) else print*,' 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) end if beta0=betst c close(nfbeam) c return end c c c subroutine equidata implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) parameter(pi=3.14159265358979d0) parameter(nbb=1000) c parameter(np=100) character*48 stringa dimension fpol(nnw),pres(nnw),qpsi(nnw) dimension ffprim(nnw),pprim(nnw) dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw) dimension rbbbs(nbb),zbbbs(nbb) dimension rlim(nbb),zlim(nbb) c dimension rcon(2*np+1),zcon(2*np+1) c parameter(nrest=nnw+4,nzest=nnh+4) parameter(lwrk=4*(nnw+nnh)+11*(nrest+nzest)+nrest*nnh+nzest+54) parameter(liwrk=nnh+nnw+nrest+nzest+3,kspl=3) dimension fvpsi(nnw*nnh),cc(nnw*nnh),ffvpsi(nnw*nnh) dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk) parameter(lwrkbsp=4*(nnw+nnh),liwrkbsp=nnw+nnh) dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp) parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh) parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) parameter(lw11=nnw*3+nnh*3+nnw*nnh,ldiwrk=nnw+nnh) dimension cc01(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) dimension iwrkd(ldiwrk) parameter(lwrkf=nnw*4+nrest*16) dimension tfp(nrest),cfp(nrest),wrkf(lwrkf),iwrkf(nrest),wf(nnw) dimension fpoli(nnw) c common/pareq1/psia common/pareq1t/psiant,psinop common/cent/btrcen,rcen common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/psinr/psinr common/qpsi/qpsi common/cpsin/rv,zv,psin common/cpsi/psi common/eqnn/nr,nz,npp,nintp common/ipsn/ipsinorm common/sspl/sspl common/nfile/neqdsk,nprof common/bound/zbmin,zbmax common/sgnib/sgnbphi,sgniphi common/factb/factb common/ixp/ixp common/icocos/icocos common/coffeqt/tr,tz common/coffeqtp/tfp common/coffeq/cc common/coffeqd/cc01,cc10,cc20,cc02,cc11 common/coffeqn/nsrt,nszt,nsft common/cofffp/cfp common/fpas/fpolas common/rhotsx/rhotsx common/rrtor/rrtor common/cnt/rup,zup,rlw,zlw common/limiter/rlim,zlim,nlim c c read from file eqdsk c (see http://fusion.gat.com/efit/g_eqdsk.html) c c spline interpolation of psi and derivatives c if(icocos.gt.0) then read (neqdsk,'(a48,3i4)') stringa,idum,nr,nz else read (neqdsk,*) nr,nz end if if(ipsinorm.eq.0) then read (neqdsk,2020) drnr1,dznz1,rcen,rleft,zmid read (neqdsk,2020) rmaxis,zmaxis,psiaxis,psiedge,btrcen read (neqdsk,2020) current,xdum,xdum,xdum,xdum read (neqdsk,2020) xdum,xdum,xdum,xdum,xdum read (neqdsk,2020) (fpol(i),i=1,nr) read (neqdsk,2020) (pres(i),i=1,nr) read (neqdsk,2020) (ffprim(i),i=1,nr) read (neqdsk,2020) (pprim(i),i=1,nr) read (neqdsk,2020) ((psi(i,j),i=1,nr),j=1,nz) read (neqdsk,2020) (qpsi(i),i=1,nr) else read (neqdsk,*) drnr1,dznz1,rcen,rleft,zmid read (neqdsk,*) rmaxis,zmaxis,psiaxis,psiedge,btrcen read (neqdsk,*) current,xdum,xdum,xdum,xdum read (neqdsk,*) xdum,xdum,xdum,xdum,xdum read (neqdsk,*) (fpol(i),i=1,nr) read (neqdsk,*) (pres(i),i=1,nr) read (neqdsk,*) (ffprim(i),i=1,nr) read (neqdsk,*) (pprim(i),i=1,nr) read (neqdsk,*) ((psin(i,j),i=1,nr),j=1,nz) read (neqdsk,*) (qpsi(i),i=1,nr) end if 2020 format (5e16.9) c c compensate for different reference systems c icocosmod=mod(icocos,10) if (mod(icocos,2).eq.0.and.icocosmod.gt.0) then c icocos mod 2 = 0: toroidal angle phi CW (opposite to gray convention) btrcen=-btrcen current=-current do i=1,nr fpol(i)=-fpol(i) end do end if c if (icocosmod.eq.1 .or. icocosmod.eq.4 .or. & icocosmod.eq.5 .or. icocosmod.eq.8) then c icocos mod 10 = 1,4,5,8: psi increasing with CCW Ip c icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip psiedge=-psiedge psiaxis=-psiaxis if (ipsinorm.eq.0) then do j=1,nz do i=1,nr psi(i,j)=-psi(i,j) end do end do end if end if c c add check for Ip/psi and B0/Fpol sign consistency? c current=sign(current,psiax-psiedge) btrcen=sign(btrcen,fpol(nr)) c c length in m !!! c dr=drnr1/dble(nr-1) dz=dznz1/dble(nz-1) rv(1)=rleft zv(1)=zmid-dznz1/2.0d0 dpsinr=1.0d0/dble(nr-1) c do i=1,nr psinr(i)=(i-1)*dpsinr rv(i)=rv(1)+(i-1)*dr end do c do j=1,nz zv(j)=zv(1)+(j-1)*dz end do c rmnm=rv(1) rmxm=rv(nr) zmnm=zv(1) zmxm=zv(nz) c psi function psia0=psiedge-psiaxis c icocos=0: adapt psi to force Ip sign, otherwise maintain psi if (icocosmod.ne.0) sgniphi=sign(1.0d0,-psia0) current=sign(current,sgniphi) psia=-sgniphi*abs(psia0)*factb c icocos>10: input psi is in Wb c icocos<10: input psi is in Wb/rad (gray convention) if (icocos.ge.10) psia=psia/(2.0d0*pi) c c do j=1,nz c do i=1,nr c write(80,2021) rv(i),zv(j),psi(i,j) c enddo c write(80,*) ' ' c enddo do j=1,nz do i=1,nr if(ipsinorm.eq.0) then psin(i,j)=(psi(i,j)-psiaxis)/(psia0) psi(i,j)=psin(i,j)*psia else psi(i,j)=psin(i,j)*psia end if ij=nz*(i-1)+j fvpsi(ij)=psin(i,j) enddo enddo c c spline fitting/interpolation of psin(i,j) and derivatives c iopt=0 call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, . wrk,lwrk,iwrk,liwrk,ier) c if ier=-1 data are fitted using sspl=0 if(ier.eq.-1) then sspl=0.0d0 call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, . wrk,lwrk,iwrk,liwrk,ier) end if nsrt=nsr nszt=nsz if (sspl.gt.0.0d0) then call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi, . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier) c do j=1,nz do i=1,nr ij=nz*(i-1)+j psin(i,j)=ffvpsi(ij) psi(i,j)=psin(i,j)*psia c write(79,2021) rv(i),zv(j),psin(i,j) enddo c write(79,*) ' ' enddo end if 2021 format(5(1x,e16.9)) c nur=1 nuz=0 call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, . ffvpsi,cc10,lw10,iwrkd,ldiwrk,ier) c nur=0 nuz=1 call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, . ffvpsi,cc01,lw01,iwrkd,ldiwrk,ier) c nur=2 nuz=0 call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, . ffvpsi,cc20,lw20,iwrkd,ldiwrk,ier) c nur=0 nuz=2 call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, . ffvpsi,cc02,lw02,iwrkd,ldiwrk,ier) c nur=1 nuz=1 call parder(tr,nsr,tz,nsz,cc,kspl,kspl,nur,nuz,rv,nr,zv,nz, . ffvpsi,cc11,lw11,iwrkd,ldiwrk,ier) c c scaling of f_poloidal c c icocos=0: adapt fpol to force Ip sign, otherwise maintain fpol if (icocosmod.ne.0) sgnbphi=sign(1.0d0,fpol(nr)) btrcen=sign(btrcen,sgnbphi) do i=1,nr fpol(i)=sgnbphi*abs(fpol(i))*factb wf(i)=1.0d0 end do wf(nr)=1.0d2 c c spline interpolation of fpol(i) c iopt=0 xb=0.0d0 xe=1.0d0 ssfp=0.01d0 call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsft, . tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier) c call splev(tfp,nsft,cfp,3,psinr,fpoli,nr,ier) fpolas=fpoli(nr) c c read plasma boundaries from eqdsk file c read (neqdsk,*) nbbbs,nlim if(nbbbs.gt.0) then if(ipsinorm.eq.1) . read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs) if(ipsinorm.eq.0) . read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs) c do i=1,nbbbs c write(51,*) rbbbs(i),zbbbs(i) c end do end if if(nlim.gt.0) then if(ipsinorm.eq.1) . read (neqdsk,*) (rlim(i),zlim(i),i=1,nlim) if(ipsinorm.eq.0) . read (neqdsk,2020) (rlim(i),zlim(i),i=1,nlim) end if c c compute max and min z of last closed surface c rbmin=rax rbmax=rax if (nbbbs.gt.1) then zbmin=1.0d+30 zbmax=-1.0d+30 do i=1,nbbbs if(zbbbs(i).le.zbmin) then zbmin=zbbbs(i) rbmin=rbbbs(i) end if if(zbbbs(i).ge.zbmax) then zbmax=zbbbs(i) rbmax=rbbbs(i) end if end do else zbmin=-1.0d+30 zbmax=1.0d+30 end if if(zbmin.le.zmnm) zbmin=zbmin+dz if(rbmin.le.rmnm) rbmin=rbmin+dr if(zbmax.ge.zmxm) zbmax=zbmax-dz if(rbmax.ge.rmxm) rbmax=rbmax-dr c c start with uncorrected normalized psi c psinop=0.0d0 psinxp=1.0d0 psiant=1.0d0 c c search for O-point c call points_ox(rmaxis,zmaxis,rmop,zmop,psinoptmp,info) rmaxis=rmop zmaxis=zmop print'(a,2f8.4,es12.5)','O-point',rmop,zmop,psinoptmp c c search for X-point if ixp not = 0 c if(ixp.ne.0) then if(ixp.lt.0) then r10=rbmin z10=zbmin call points_ox(r10,z10,rxp,zxp,psinxptmp,info) if(psinxp.ne.-1.0d0) then print'(a,2f8.4,es12.5)','X-point',rxp,zxp,psinxptmp rbmin=rxp zbmin=zxp psinop=psinoptmp psinxp=psinxptmp psiant=psinxp-psinop psin1=1.0d0 r10=rmaxis z10=(zbmax+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) rbmax=r1 zbmax=z1 else ixp=0 c print'(a)','no X-point' end if else r10=rmop z10=zbmax call points_ox(r10,z10,rxp,zxp,psinxptmp,info) if(psinxp.ne.-1.0d0) then print'(a,2f8.4,e16.8)','X-point',rxp,zxp,psinxptmp zbmax=zxp psinop=psinoptmp psinxp=psinxptmp psiant=psinxp-psinop psin1=1.0d0 z10=(zbmin+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) zbmin=z1 else ixp=0 c print'(a)','no X-point' end if end if end if c if (ixp.eq.0) then psin1=1.0d0 psinop=psinoptmp psiant=psin1-psinop r10=rmaxis z10=(zbmax+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) zbmax=z1 rbmax=r1 z10=(zbmin+zmaxis)/2.0d0 call points_tgo(r10,z10,r1,z1,psin1,info) zbmin=z1 rbmin=r1 print'(a,4f8.4)','no X-point ',rbmin,zbmin,rbmax,zbmax end if print*,' ' 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) c compute normalized rho_tor from eqdsk q profile call rhotor(nr) c phitedge=abs(psia)*rhotsx*2*pi c rrtor=sqrt(phitedge/abs(btrcen)/pi) c compute flux surface averaged quantities rup=rmaxis rlw=rmaxis zup=zmaxis+(zbmax-zmaxis)/10.0d0 zlw=zmaxis-(zmaxis-zbmin)/10.0d0 call flux_average c ipr=1 c call contours_psi(1.0d0,np,rcon,zcon,ipr) c do ii=1,2*np+1 c write(52,*) rcon(ii), zcon(ii) c end do c c locate psi surface for q=1.5 and q=2 rup=rmaxis rlw=rmaxis zup=(zbmax+zmaxis)/2.0d0 zlw=(zmaxis+zbmin)/2.0d0 q2=2.0d0 q15=1.5d0 call vmaxmini(qpsi,nr,qmin,qmax,iqmn,iqmx) if (q15.gt.qmin.and.q15.lt.qmax) then call surfq(q15,psi15) rhot15=frhotor(psi15) print'(3(a,f8.5))','psi_1.5 = ',psi15,' rhop_1.5 = ' . ,sqrt(psi15),' rhot_1.5 = ',rhot15 end if if (q2.gt.qmin.and.q2.lt.qmax) then call surfq(q2,psi2) rhot2=frhotor(psi2) print'(3(a,f8.5))','psi_2 = ',psi2,' rhop_2 = ' . ,sqrt(psi2),' rhot_2 = ',rhot2 end if c c locate btot=bres c call bfield_res c return end c c c subroutine points_ox(rz,zz,rf,zf,psinvf,info) implicit real*8 (a-h,o-z) parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) external fcnox common/psival/psinv xvec(1)=rz xvec(2)=zz tol = sqrt(dpmpar(1)) call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) if(info.gt.1) then print'(a,i2,a,2f8.4)',' info subr points_ox =',info, . ' O/X coord.',xvec end if rf=xvec(1) zf=xvec(2) psinvf=psinv return end c c c subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) implicit real*8 (a-h,o-z) dimension x(n),fvec(n),fjac(ldfjac,n) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/pareq1/psia call equinum(x(1),x(2)) if (iflag.ne.2) then fvec(1) = dpsidr/psia fvec(2) = dpsidz/psia else fjac(1,1) = ddpsidrr/psia fjac(1,2) = ddpsidrz/psia fjac(2,1) = ddpsidrz/psia fjac(2,2) = ddpsidzz/psia end if return end c c c subroutine points_tgo(rz,zz,rf,zf,psin,info) implicit real*8 (a-h,o-z) parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) external fcntgo common/cnpsi/h h=psin xvec(1)=rz xvec(2)=zz tol = sqrt(dpmpar(1)) call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) if(info.gt.1) then end if rf=xvec(1) zf=xvec(2) return end c c c subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag) implicit real*8 (a-h,o-z) dimension x(n),fvec(n),fjac(ldfjac,n) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/psival/psinv common/cnpsi/h common/pareq1/psia call equinum(x(1),x(2)) if (iflag.ne.2) then fvec(1) = psinv-h fvec(2) = dpsidr/psia else fjac(1,1) = dpsidr/psia fjac(1,2) = dpsidz/psia fjac(2,1) = ddpsidrr/psia fjac(2,2) = ddpsidrz/psia end if return end c c c subroutine print_prof implicit real*8 (a-h,o-z) parameter(nnw=501,eps=1.d-4) dimension psinr(nnw),qpsi(nnw) c common/psinr/psinr common/qpsi/qpsi common/eqnn/nr,nz,npp,nintp common/dens/dens,ddens c write(55,*) ' #psi rhot ne Te q Jphi' psin=0.0d0 rhop=0.0d0 rhot=0.0d0 call density(psin) call tor_curr_psi(eps,ajphi) te=temp(psin) qq=qpsi(1) c write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6 c nst=nr do i=2,nst psin=dble(i-1)/dble(nst-1) rhop=sqrt(psin) c call density(psin) te=temp(psin) c ips=int((nr-1)*psin+1) if(i.lt.nst) then call intlin(psinr(ips),qpsi(ips),psinr(ips+1),qpsi(ips+1), . psin,qq) else qq=qpsi(nr) end if rhot=frhotor(psin) call tor_curr_psi(psin,ajphi) write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6 end do c return 111 format(12(1x,e12.5)) end c c c subroutine surfq(qval,psival) implicit real*8 (a-h,o-z) parameter(nnw=501) parameter(ncnt=100,ncntt=2*ncnt+1) dimension psinr(nnw),qpsi(nnw) dimension rcon(ncntt),zcon(ncntt) c common/psinr/psinr common/qpsi/qpsi common/eqnn/nr,nz,npp,nintp c c locate psi surface for q=qval c call locate(qpsi,nr,qval,i1) call intlin(qpsi(i1),psinr(i1),qpsi(i1+1),psinr(i1+1), . qval,psival) ipr=1 call contours_psi(psival,ncnt,rcon,zcon,ipr) return end c c c subroutine bfield_res implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) dimension rv(nnw),zv(nnh),psin(nnw,nnh) dimension btotal(nnw,nnh) parameter(icmx=2002) dimension rrcb(icmx),zzcb(icmx),ncpts(10) c common/cpsin/rv,zv,psin common/eqnn/nr,nz,npp,nintp common/parbres/bres common/btt/btotal c c Btotal on psi grid c btmx=-1.0d30 btmn=1.0d30 do j=1,nr rrj=rv(j) do k=1,nz zzk=zv(k) call bfield(rrj,zzk,bbphi,bbr,bbz) btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2) if(btotal(j,k).ge.btmx) btmx=btotal(j,k) if(btotal(j,k).le.btmn) btmn=btotal(j,k) c write(90,113) j,rrj,zzk,btotal(j,k) enddo c write(90,*) ' ' enddo c c compute Btot=Bres and Btot=Bres/2 c write(70,*)'#i Btot R z' do n=1,5 bbb=bres/dble(n) if (bbb.ge.btmn.and.bbb.le.btmx) then call cniteq(bbb,nconts,ncpts,nctot,rrcb,zzcb,1) do inc=1,nctot write(70,113) inc,bbb,rrcb(inc),zzcb(inc) end do end if write(70,*) ' ' end do c return 113 format(i6,12(1x,e12.5)) end c c subroutine profdata implicit real*8 (a-h,o-z) parameter(npmx=250,npest=npmx+4) dimension psrad(npmx),terad(npmx),derad(npmx),zfc(npmx) dimension ct(npmx,4),cz(npmx,4) parameter(lwrkf=npmx*4+npest*16) dimension tfn(npest),cfn(npest),wrkf(lwrkf),iwrkf(npest),wf(npmx) dimension densi(npest),ddensi(npest),d2densi(npest),wrkfd(npest) c common/densbnd/psdbnd common/denspp/psnpp,aa,bb,cc common/eqnn/nr,nz,npp,nintp common/iipr/iprof common/nfile/neqdsk,nprof common/crad/psrad,derad,terad,zfc common/coffte/ct common/coffz/cz common/factb/factb common/coffdt/tfn common/coffdnst/nsfd common/cofffn/cfn common/scal/iscal common/facttn/factt,factn c c read plasma profiles from file if iprof>0 c if (iscal.eq.0) then aat=2.0d0/3.0d0 aan=4.0d0/3.0d0 else aat=1.0d0 aan=1.0d0 end if ffact=factb if(iscal.eq.2) ffact=1.0d0 c if (iprof.gt.0) then read(nprof,*) npp read(nprof,*) psrad0,terad0,derad0,zfc0 if(psrad0.ne.0.0d0) psrad0=0.0d0 psrad(1)=psrad0 terad(1)=terad0*ffact**aat*factt derad(1)=derad0*ffact**aan*factn zfc(1)=zfc0 wf(1)=1.0d0 do i=2,npp read(nprof,*) psradi,teradi,deradi,zfci psrad(i)=psradi terad(i)=teradi*ffact**aat*factt derad(i)=deradi*ffact**aan*factn zfc(i)=zfci wf(i)=1.0d0 end do c c spline approximation of temperature and Zeff c iopt=0 call difcsn(psrad,terad,npmx,npp,iopt,ct,ier) c iopt=0 call difcsn(psrad,zfc,npmx,npp,iopt,cz,ier) c c spline approximation of density c iopt=0 xb=0.0d0 xe=psrad(npp) kspl=3 sspl=.001d0 c call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd, . tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier) c call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier) nu=1 call splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier) dnpp=densi(npp) ddnpp=ddensi(npp) c nu=2 call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier) d2dnpp=d2densi(npp) if(derad(npp).eq.0.0d0) then psdbnd=psrad(npp) else psnpp=psrad(npp) dpsb=-psnpp+psdbnd nn=3 nn1=nn+1 nn2=nn+2 aa=(nn1*nn2*dnpp+2*nn1*ddnpp*dpsb+d2dnpp*dpsb**2) aa=aa/(-dpsb)**nn/2.0d0 bb=-(nn*nn2*dnpp+(2*nn+1)*ddnpp*dpsb+d2dnpp*dpsb**2) bb=bb/(-dpsb)**nn1 cc=(nn1*nn*dnpp+2*nn*ddnpp*dpsb+d2dnpp*dpsb**2) cc=cc/(-dpsb)**nn2/2.0d0 end if c end if c return end c c c subroutine rhotor(nr) implicit real*8(a-h,o-z) parameter(nnw=501) dimension psinr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),crhot(nnw,4) common/psinr/psinr common/qpsi/qpsi common/rhotsx/rhotsx common/crhot/crhot c c normalized toroidal rho : ~ Integral q(psi) dpsi c iopt=0 call difcsn(psinr,qpsi,nnw,nr,iopt,cq,ier) c rhotnr(1)=0.0d0 do k=1,nr-1 dx=psinr(k+1)-psinr(k) drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0d0+dx*(cq(k,3)/3.0d0 . +dx*cq(k,4)/4.0d0))) rhotnr(k+1)=rhotnr(k)+drhot end do rhotsx=rhotnr(nr) do k=1,nr rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr)) end do c c spline interpolation of rhotor c iopt=0 call difcs(psinr,rhotnr,nr,iopt,crhot,ier) return end function frhotor_eq(psi) implicit real*8(a-h,o-z) parameter(nnw=501) dimension psinr(nnw),crhot(nnw,4) common/psinr/psinr common/eqnn/nr,nz,npp,nintp common/crhot/crhot c irt=int((nr-1)*psi+1) if(irt.eq.0) irt=1 if(irt.eq.nr) irt=nr-1 dps=psi-psinr(irt) frhotor_eq=spli(crhot,nr,irt,dps) return end function frhotor(psi) implicit real*8(a-h,o-z) frhotor=frhotor_eq(psi) return end function frhotor_av(psi) implicit real*8(a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp),crhotq(nnintp,4) common/pstab/rpstab common/eqnn/nr,nz,npp,nintp 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 dps=rpsi-rpstab(ip) frhotor_av=spli(crhotq,nintp,ip,dps) return end subroutine cniteq(h, ncon, npts, icount, rcon, zcon,ichoi) implicit real*8 (a-h,o-z) c c v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. c (based on an older code) c parameter(nnw=501,nnh=501,icmx=2002,nna=nnw*nnh) dimension a(nna),ja(3,2),lx(1000),npts(10) dimension rcon(icmx),zcon(icmx) dimension rqgrid(nnw),zqgrid(nnh),psin(nnw,nnh),btotal(nnw,nnh) c common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/cpsin/rqgrid,zqgrid,psin common/btt/btotal common/eqnn/nr,nz,npp,nintp c data px/0.5d0/ c if(ichoi.eq.0) then do j=1,nz do i=1,nr a(nr*(j-1)+i)=psin(i,j) enddo enddo endif c if(ichoi.eq.1) then do j=1,nz do i=1,nr a(nr*(j-1)+i)=btotal(i,j) enddo enddo endif c do ico=1,icmx rcon(ico)=0.0d0 zcon(ico)=0.0d0 enddo c nrqmax=nr nreq=nr nzeq=nz drgrd=dr dzgrd=dz c ncon = 0 do i=1,10 npts(i) = 0 end do iclast = 0 icount = 0 mpl=0 ix=0 mxr = nrqmax * (nzeq - 1) n1 = nreq - 1 c do 3 jx=2,n1 do 3 jm=jx,mxr,nrqmax j = jm + nrqmax ah=a(j)-h if(ah) 60,60,61 60 if(a(jm)-h) 62,62,1 61 if(a(jm)-h) 1,1,63 1 ix=ix+1 lx(ix)=-j if(ah) 62,62,63 62 if(a(j-1)-h) 3,3,4 63 if(a(j-1)-h) 4,4,3 4 ix=ix+1 lx(ix)=j 3 continue c do 79 jm=nreq,mxr,nrqmax j = jm + nrqmax ah=a(j)-h if(ah) 64,64,65 64 if(a(j-1)-h) 66,66,76 65 if(a(j-1)-h) 76,76,67 76 ix=ix+1 lx(ix)=j if(ah) 66,66,67 66 if(a(jm)-h) 79,79,78 67 if(a(jm)-h) 78,78,79 78 ix=ix+1 lx(ix)=-j 79 continue c do 75 jm=1,mxr,nrqmax j = jm + nrqmax if(a(j)-h) 68,68,69 68 if(a(jm)-h)75,75,71 69 if(a(jm)-h)71,71,75 71 ix=ix+1 lx(ix) =-j 75 continue c do 73 j=2,nreq if(a(j)-h) 58,58,59 58 if(a(j-1)-h) 73,73,72 59 if(a(j-1)-h) 72,72,73 72 ix=ix+1 lx(ix)=j 73 continue c if(ix) 50,50,8 108 if(mpl.lt.4) go to 8 ncon = ncon + 1 npts(ncon) = icount - iclast iclast = icount 8 in=ix jx=lx(in) jfor=0 lda=1 ldb=2 30 if(jx) 21,22,22 21 jabs=-jx jnb = jabs - nrqmax go to 23 22 jabs=jx jnb=jabs-1 23 adn=a(jabs)-a(jnb) if(adn) 24,9,24 24 px=(a(jabs)-h)/adn 9 kx = (jabs - 1) / nrqmax ikx = jabs - nrqmax * kx - 1 if(jx) 25,26,26 25 x = drgrd * ikx y = dzgrd * (kx - px) go to 27 26 x = drgrd * (ikx - px) y = dzgrd * kx 27 continue icount = icount + 1 rcon(icount) = x + rqgrid(1) zcon(icount) = y + zqgrid(1) mpl= icount itm = 1 ja(1,1) = jabs + nrqmax j=1 if(jx) 10,10,11 10 ja(1,1) = -jabs-1 j=2 11 ja(2,1)=-ja(1,1) ja(3,1) = -jx + 1 - nrqmax ja(3,2) = -jx ja(j,2) = jabs - nrqmax k= 3-j ja(k,2) = 1-jabs if (kx) 14,14,39 39 if(ikx) 14,14,36 36 if(ikx + 1 - nreq) 35, 37, 37 37 if(jx) 38,38,35 35 if(jfor) 28,29,28 28 do 13 i=1,3 if(jfor-ja(i,2)) 13,14,13 13 continue 38 lda=2 go to 15 14 lda=1 15 ldb=lda 29 do 18 k=1,3 do 18 l=lda,ldb do 16 i=1,ix if(lx(i)-ja(k,l)) 16,17,16 16 continue go to 18 17 itm=itm+1 inext= i if(jfor) 19,33,19 33 if(itm .gt. 3) goto 20 18 continue 19 lx(in)=0 if(itm .eq. 1) goto 6 20 jfor=jx jx=lx(inext) in = inext go to 30 6 if(lx(ix)) 108,7,108 7 ix= ix-1 if(ix) 51,51,6 51 if(mpl.lt.4) go to 50 ncon = ncon + 1 npts(ncon) = icount - iclast iclast = icount 50 continue c return end c c c subroutine contours_psi(h,np,rcn,zcn,ipr) implicit real*8 (a-h,o-z) parameter(pi=3.14159265358979d0) parameter(mest=4,kspl=3) parameter(nnw=501,nnh=501) parameter(nrest=nnw+4,nzest=nnh+4) dimension rcn(2*np+1),zcn(2*np+1) dimension cc(nnw*nnh),tr(nrest),tz(nzest) dimension czc(nrest),zeroc(mest) c common/pareq1/psia common/pareq1t/psiant,psinop common/coffeqn/nsr,nsz,nsft common/coffeq/cc common/coffeqt/tr,tz common/cnt/rup,zup,rlw,zlw common/rwallm/rwallm c ra=rup rb=rlw za=zup zb=zlw call points_tgo(ra,za,rup,zup,h,info) call points_tgo(rb,zb,rlw,zlw,h,info) th=pi/dble(np) rcn(1)=rlw zcn(1)=zlw rcn(2*np+1)=rlw zcn(2*np+1)=zlw rcn(np+1)=rup zcn(np+1)=zup do ic=2,np zc=zlw+(zup-zlw)*(1.0d0-cos(th*(ic-1)))/2.0d0 iopt=1 call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) print*,' profil =',ier val=h*psiant+psinop call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) if (zeroc(1).gt.rwallm) then rcn(ic)=zeroc(1) zcn(ic)=zc rcn(2*np+2-ic)=zeroc(2) zcn(2*np+2-ic)=zc else rcn(ic)=zeroc(2) zcn(ic)=zc rcn(2*np+2-ic)=zeroc(3) zcn(2*np+2-ic)=zc end if end do if (ipr.gt.0) then do ii=1,2*np+1 write(71,111) ii,h,rcn(ii),zcn(ii) end do write(71,*) write(71,*) end if return 111 format(i6,12(1x,e12.5)) 99 format(12(1x,e12.5)) end c c c subroutine flux_average implicit real*8 (a-h,o-z) real*8 lam parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41) parameter(zero=0.0d0,one=1.0d0) parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi)) parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ . njest*nnintp+nlest+54) parameter(kwrk=nnintp+nlam+njest+nlest+3) parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) dimension pstab(nnintp),varea(nnintp),vvol(nnintp),rpstab(nnintp) dimension rri(nnintp),rbav(nnintp),bav(nnintp),rhotqv(nnintp) dimension bmxpsi(nnintp),bmnpsi(nnintp),ffc(nnintp) dimension vcurrp(nnintp),vajphiav(nnintp),qqv(nnintp) dimension rcon(2*ncnt+1),zcon(2*ncnt+1) dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1) dimension cbmx(nnintp,4),cbmn(nnintp,4),crbav(nnintp,4) dimension cvol(nnintp,4),crri(nnintp,4),carea(nnintp,4) dimension cratjpl(nnintp,4),cratja(nnintp,4),cratjb(nnintp,4) dimension cfc(nnintp,4),crhotq(nnintp,4) dimension vratjpl(nnintp),vratja(nnintp),vratjb(nnintp) dimension alam(nlam),fhlam(nnintp,nlam) dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam) dimension tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)) dimension iwrk(kwrk),wrk(lwrk) dimension ch01(lw01),weights(nlam) c common/cent/btrcen,rcen common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/eqnn/nr,nz,npp,nintp common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv c common/pstab/rpstab common/flav/vvol,rri,rbav,bmxpsi,bmnpsi common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc common/cratj/cratja,cratjb,cratjpl common/crhotq/crhotq common/cnt/rup,zup,rlw,zlw common/bound/zbmin,zbmax common/rarea/rarea c common/coffvl/ch common/coffdvl/ch01 common/coffvlt/tjp,tlm common/coffvln/njpt,nlmt c computation of flux surface averaged quantities write(71,*)' #i psin R z' nintp=nnintp ninpr=(nintp-1)/10 dlam=1.0d0/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam fhlam(1,l)=sqrt(1.0d0-alam(l)) ffhlam(l)=fhlam(1,l) dffhlam(l)=-0.5d0/sqrt(1.0d0-alam(l)) weights(l)=1.0d0 end do weights(1)=0.5d0 weights(nlam)=0.5d0 alam(nlam)=1.0d0 fhlam(1,nlam)=0.0d0 ffhlam(nlam)=0.0d0 dffhlam(nlam)=-99999.0d0 jp=1 anorm=2.0d0*pi*rmaxis/abs(btaxis) b2av=btaxis**2 dvdpsi=2.0d0*pi*anorm dadpsi=2.0d0*pi/abs(btaxis) ratio_cdator=abs(btaxis/btrcen) ratio_cdbtor=1.0d0 ratio_pltor=1.0d0 qq=1.0d0 fc=1.0d0 pstab(1)=0.0d0 rpstab(1)=0.0d0 rhotqv(1)=0.0d0 vcurrp(1)=0.0d0 vajphiav(1)=0.0d0 bmxpsi(1)=abs(btaxis) bmnpsi(1)=abs(btaxis) bav(1)=abs(btaxis) rbav(1)=1.0d0 rri(1)=rmaxis varea(1)=0.0d0 vvol(1)=0.0d0 vratjpl(1)=ratio_pltor vratja(1)=ratio_cdator vratjb(1)=ratio_cdbtor ffc(1)=fc rhot2q=0.0d0 C rup=rmaxis C rlw=rmaxis C zup=(zbmax+zmaxis)/2.0d0 C zlw=(zmaxis+zbmin)/2.0d0 do jp=2,nintp height=dble(jp-1)/dble(nintp-1) if(jp.eq.nintp) height=0.9999d0 ipr=0 jpr=mod(jp,ninpr) if(jpr.eq.1) ipr=1 height2=height*height call contours_psi(height2,ncnt,rcon,zcon,ipr) c r2iav=0.0d0 anorm=0.0d0 dadpsi=0.0d0 currp=0.0d0 b2av=0.0d0 area=0.0d0 volume=0.0d0 ajphiav=0.0d0 bbav=0.0d0 bmmx=-1.0d+30 bmmn=1.0d+30 call bfield(rcon(1),zcon(1),bphi,brr,bzz) call tor_curr(rcon(1),zcon(1),ajphi0) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) bv(1)=btot0 bpv(1)=bpoloid0 rpsim0=rcon(1) do inc=1,ncntt-1 inc1=inc+1 dla=sqrt((rcon(inc)-rmaxis)**2+(zcon(inc)-zmaxis)**2) dlb=sqrt((rcon(inc1)-rmaxis)**2+(zcon(inc1)-zmaxis)**2) dlp=sqrt((rcon(inc1)-rcon(inc))**2+ . (zcon(inc1)-zcon(inc))**2) drc=(rcon(inc1)-rcon(inc)) c compute length, area and volume defined by psi=height^2 ph=0.5d0*(dla+dlb+dlp) area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) area=area+sqrt(area2) rzp=rcon(inc1)*zcon(inc1) rz=rcon(inc)*zcon(inc) volume=pi*(rzp+rz)*drc+volume c compute line integral on the contour psi=height^2 rpsim=rcon(inc1) zpsim=zcon(inc1) call bfield(rpsim,zpsim,bphi,brr,bzz) call tor_curr(rpsim,zpsim,ajphi) btot=sqrt(bphi**2+brr**2+bzz**2) bpoloid=sqrt(brr**2+bzz**2) dlpv(inc)=dlp bv(inc1)=btot bpv(inc1)=bpoloid dlph=0.5d0*dlp anorm=anorm+dlph*(1.0d0/bpoloid+1.0d0/bpoloid0) dadpsi=dadpsi+dlph* . (1.0d0/(bpoloid*rpsim)+1.0d0/(bpoloid0*rpsim0)) currp=currp+dlph*(bpoloid+bpoloid0) b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) r2iav=r2iav+dlph* . (1.0d0/(bpoloid*rpsim**2)+1.0d0/(bpoloid0*rpsim0**2)) ajphiav=ajphiav+dlph* . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) ajphi0=ajphi rpsim0=rpsim bpoloid0=bpoloid btot0=btot c computation maximum/minimum B values on given flux surface if(btot.le.bmmn) bmmn=btot if(btot.ge.bmmx) bmmx=btot end do c bav= [T] , b2av= [T^2] , rbav=/b_min c anorm = int d l_p/B_p = dV/dpsi/(2pi) c r2iav=<1/R^2> [m^-2] , c riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)), c rri = /(|R B_tor|<1/R^2>) , used to compute I_tor [m^-1] c currp = plasma current within psi=const bbav=bbav/anorm r2iav=r2iav/anorm dvdpsi=2.0d0*pi*anorm riav=dadpsi/anorm b2av=b2av/anorm vcurrp(jp)=ccj*currp vajphiav(jp)=ajphiav/dadpsi c area == varea, volume == vvol c flux surface minor radius == (area/pi)^1/2 c ratio_cdator = Jcd_astra/J_phi Jcd_astra = /B0 c ratio_cdbtor = Jcd_jintrac/J_phi Jcd_jintrac = / c ratio_pltor = Jcd_||/J_phi Jcd_|| = pstab(jp)=height2 rpstab(jp)=height vvol(jp)=abs(volume) varea(jp)=area bav(jp)=bbav rbav(jp)=bbav/bmmn bmxpsi(jp)=bmmx bmnpsi(jp)=bmmn rri(jp)=bav(jp)/abs(fpolv*r2iav) ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen)) ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav)) ratio_pltor=abs(bbav*riav/(fpolv*r2iav)) vratjpl(jp)=ratio_pltor vratja(jp)=ratio_cdator vratjb(jp)=ratio_cdbtor qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) qqv(jp)=qq c c write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp) c computation of rhot from calculated q profile rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) . /dble(nintp-1) rhotqv(jp)=sqrt(rhot2q) c computation of fraction of circulating/trapped fraction fc, ft c and of function H(lambda,rhop) c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ fc=0.0d0 shlam=0.0d0 do l=nlam,1,-1 lam=alam(l) srl=0.0d0 rl2=1.0d0-lam*bv(1)/bmmx rl0=0.d0 if(rl2.gt.0) rl0=sqrt(rl2) do inc=1,ncntt-1 rl2=1.0d0-lam*bv(inc+1)/bmmx rl=0.0d0 if(rl2.gt.0) rl=sqrt(rl2) srl=srl+0.5d0*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) rl0=rl end do srl=srl/anorm dhlam=0.5d0/srl fc=fc+lam/srl*weights(l) if(l.eq.nlam) then fhlam(jp,l)=0.0d0 ffhlam(nlam*(jp-1)+l)=0.0d0 dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam else shlam=shlam+0.5d0*(dhlam+dhlam0)*dlam fhlam(jp,l)=shlam dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam end if end do fc=0.75d0*b2av/bmmx**2*fc*dlam ffc(jp)=fc ccfh=bmmn/bmmx/fc do l=1,nlam ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) end do end do write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// .' Area Vol |I_pl| qq fc ratio_cdbtor' qqv(1)=qqv(2) vajphiav(1)=vajphiav(2) do jp=1,nintp rhotqv(jp)=rhotqv(jp)/rhotqv(nintp) if(jp.eq.nintp) then rhotqv(jp)=1.0d0 rpstab(jp)=1.0d0 pstab(jp)=1.0d0 end if rhot_eq=frhotor_eq(pstab(jp)) write(56,99) pstab(jp),rhot_eq,rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp) end do rarea=sqrt(varea(nintp)/pi) c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi c used for computations of dP/dV and J_cd c spline coefficients of rhot iopt=0 call difcs(rpstab,vvol,nintp,iopt,cvol,ier) iopt=0 call difcs(rpstab,rbav,nintp,iopt,crbav,ier) iopt=0 call difcs(rpstab,rri,nintp,iopt,crri,ier) iopt=0 call difcs(rpstab,bmxpsi,nintp,iopt,cbmx,ier) iopt=0 call difcs(rpstab,bmnpsi,nintp,iopt,cbmn,ier) iopt=0 call difcs(rpstab,vratja,nintp,iopt,cratja,ier) iopt=0 call difcs(rpstab,vratjb,nintp,iopt,cratjb,ier) iopt=0 call difcs(rpstab,vratjpl,nintp,iopt,cratjpl,ier) iopt=0 call difcs(rpstab,varea,nintp,iopt,carea,ier) iopt=0 call difcs(rpstab,ffc,nintp,iopt,cfc,ier) iopt=0 call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier) c spline interpolation of H(lambda,rhop) and dH/dlambda iopt=0 s=0.0d0 call regrid(iopt,nintp,rpstab,nlam,alam,ffhlam, . zero,one,zero,one,ksp,ksp,s, . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) njpt=njp nlmt=nlm call coeff_parder(tjp,njp,tlm,nlm,ch,ksp,ksp,0,1, . ch01,lw01,ier) return 99 format(20(1x,e12.5)) end c c c subroutine vectinit implicit real*8 (a-h,o-z) 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 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) 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/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj common/pcjki/ppabs,ccci common/dcjki/didst common/nray/nrayr,nrayth common/nstep/nstep common/tau1v/tau1v c if(nstep.gt.nmx) nstep=nmx if(nrayr.gt.jmx) nrayr=jmx if(nrayth.gt.kmx) nrayth=kmx c do i=1,nstep do k=1,nrayth do j=1,nrayr psjki(j,k,i)=0.0d0 tauv(j,k,i)=0.0d0 alphav(j,k,i)=0.0d0 pdjki(j,k,i)=0.0d0 ppabs(j,k,i)=0.0d0 didst(j,k,i)=0.0d0 ccci(j,k,i)=0.0d0 currj(j,k,i)=0.0d0 iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 tau1v(j,k)=0.0d0 end do end do end do c if(iwarm.gt.1) then dt=2.0d0*tmax/dble(npts) do i = 1, npts+1 ttv(i) = -tmax+dble(i-1)*dt extv(i)=exp(-ttv(i)*ttv(i)) end do end if c return end subroutine vectinit2 implicit real*8 (a-h,o-z) 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 currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx) common/iiv/iiv common/iov/iop,iow common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj common/pcjki/ppabs,ccci common/dcjki/didst common/nray/nrayr,nrayth common/nstep/nstep do i=1,nstep do k=1,nrayth do j=1,nrayr psjki(j,k,i)=0.0d0 tauv(j,k,i)=0.0d0 alphav(j,k,i)=0.0d0 pdjki(j,k,i)=0.0d0 ppabs(j,k,i)=0.0d0 didst(j,k,i)=0.0d0 ccci(j,k,i)=0.0d0 currj(j,k,i)=0.0d0 iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 end do end do end do return end c c c subroutine paraminit implicit real*8(a-h,o-z) c common/istep/istep 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 c c subroutine updatepos implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36) dimension xc(3,jmx,kmx),xco(3,jmx,kmx) dimension du1(3,jmx,kmx),du1o(3,jmx,kmx) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) c common/nray/nrayr,nrayth common/grco/xco,du1o common/grc/xc,du1 common/wrk/ywrk,ypwrk c do j=1,nrayr do k=1,nrayth if(j.eq.1.and.k.gt.1) then xco(1,j,k)=xco(1,j,1) xco(2,j,k)=xco(2,j,1) xco(3,j,k)=xco(3,j,1) xc(1,j,k)=xc(1,j,1) xc(2,j,k)=xc(2,j,1) xc(3,j,k)=xc(3,j,1) else xco(1,j,k)=xc(1,j,k) xco(2,j,k)=xc(2,j,k) xco(3,j,k)=xc(3,j,k) xc(1,j,k)=ywrk(1,j,k) xc(2,j,k)=ywrk(2,j,k) xc(3,j,k)=ywrk(3,j,k) end if du1o(1,j,k)=du1(1,j,k) du1o(2,j,k)=du1(2,j,k) du1o(3,j,k)=du1(3,j,k) end do end do c return end c c c subroutine gradi implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36) dimension dffiu(jmx),ddffiu(jmx) dimension xc(3,jmx,kmx),xco(3,jmx,kmx) dimension du1(3,jmx,kmx),du1o(3,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) dimension dgrad2v(3,jmx,kmx) dimension grad2(jmx,kmx) dimension dxv1(3),dxv2(3),dxv3(3),dgu(3) dimension dgg1(3),dgg2(3),dgg3(3) dimension df1(3),df2(3),df3(3) c common/nray/nrayr,nrayth common/fi/dffiu,ddffiu common/grco/xco,du1o common/grc/xc,du1 common/grad2jk/grad2 common/dgrad2jk/dgrad2v common/gradjk/gri common/ggradjk/ggri c c compute grad u and grad(S_I) c do k=1,nrayth do j=1,nrayr if(j.eq.1) then gri(1,j,k)=0.0d0 gri(2,j,k)=0.0d0 gri(3,j,k)=0.0d0 jp=j+1 km=k-1 if(k.eq.1) km=nrayth kp=k+1 if(k.eq.nrayth) kp=1 do iv=1,3 dxv1(iv)=xc(iv,jp,k)-xc(iv,j,k) dxv2(iv)=xc(iv,jp,kp)-xc(iv,jp,km) dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) end do call solg0(dxv1,dxv2,dxv3,dgu) else jm=j-1 km=k-1 if(k.eq.1) km=nrayth kp=k+1 if(k.eq.nrayth) kp=1 do iv=1,3 dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k) dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km) dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) end do call solg0(dxv1,dxv2,dxv3,dgu) end if du1(1,j,k)=dgu(1) du1(2,j,k)=dgu(2) du1(3,j,k)=dgu(3) gri(1,j,k)=dgu(1)*dffiu(j) gri(2,j,k)=dgu(2)*dffiu(j) gri(3,j,k)=dgu(3)*dffiu(j) grad2(j,k)=gri(1,j,k)**2+gri(2,j,k)**2+gri(3,j,k)**2 end do end do c c compute derivatives of grad u and grad(S_I) c do k=1,nrayth do j=1,nrayr if(j.eq.1) then jp=j+1 km=k-1 if(k.eq.1) km=nrayth kp=k+1 if(k.eq.nrayth) kp=1 do iv=1,3 dxv1(iv)=xc(iv,jp,kp)-xc(iv,jp,km) dxv2(iv)=xc(iv,jp,k)-xc(iv,j,k) dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) end do df1(1)=du1(1,jp,kp)-du1(1,jp,km) df1(2)=du1(1,jp,k)-du1(1,j,k) df1(3)=du1(1,j,k)-du1o(1,j,k) df2(1)=du1(2,jp,kp)-du1(2,jp,km) df2(2)=du1(2,jp,k)-du1(2,j,k) df2(3)=du1(2,j,k)-du1o(2,j,k) df3(1)=du1(3,jp,kp)-du1(3,jp,km) df3(2)=du1(3,jp,k)-du1(3,j,k) df3(3)=du1(3,j,k)-du1o(3,j,k) call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3) else jm=j-1 km=k-1 if(k.eq.1) km=nrayth kp=k+1 if(k.eq.nrayth) kp=1 do iv=1,3 dxv1(iv)=xc(iv,j,k)-xc(iv,jm,k) dxv2(iv)=xc(iv,j,kp)-xc(iv,j,km) dxv3(iv)=xc(iv,j,k)-xco(iv,j,k) end do df1(1)=du1(1,j,k)-du1(1,jm,k) df1(2)=du1(1,j,kp)-du1(1,j,km) df1(3)=du1(1,j,k)-du1o(1,j,k) df2(1)=du1(2,j,k)-du1(2,jm,k) df2(2)=du1(2,j,kp)-du1(2,j,km) df2(3)=du1(2,j,k)-du1o(2,j,k) df3(1)=du1(3,j,k)-du1(3,jm,k) df3(2)=du1(3,j,kp)-du1(3,j,km) df3(3)=du1(3,j,k)-du1o(3,j,k) call solg3(dxv1,dxv2,dxv3,df1,df2,df3,dgg1,dgg2,dgg3) end if c c derivatives of u c ux=du1(1,j,k) uy=du1(2,j,k) uz=du1(3,j,k) uxx=dgg1(1) uyy=dgg2(2) uzz=dgg3(3) uxy=(dgg1(2)+dgg2(1))/2.0d0 uxz=(dgg1(3)+dgg3(1))/2.0d0 uyz=(dgg2(3)+dgg3(2))/2.0d0 c c derivatives of S_I and Grad(S_I) c dfu=dffiu(j) dfuu=ddffiu(j) gx=ux*dfu gy=uy*dfu gz=uz*dfu gxx=dfuu*ux*ux+dfu*uxx gyy=dfuu*uy*uy+dfu*uyy gzz=dfuu*uz*uz+dfu*uzz gxy=dfuu*ux*uy+dfu*uxy gxz=dfuu*ux*uz+dfu*uxz gyz=dfuu*uy*uz+dfu*uyz c ggri(1,1,j,k)=gxx ggri(2,2,j,k)=gyy ggri(3,3,j,k)=gzz ggri(1,2,j,k)=gxy ggri(2,1,j,k)=gxy ggri(1,3,j,k)=gxz ggri(3,1,j,k)=gxz ggri(2,3,j,k)=gyz ggri(3,2,j,k)=gyz c c derivatives of |Grad(S_I)|^2 c dgrad2v(1,j,k)=2.0d0*(gx*gxx+gy*gxy+gz*gxz) dgrad2v(2,j,k)=2.0d0*(gx*gxy+gy*gyy+gz*gyz) dgrad2v(3,j,k)=2.0d0*(gx*gxz+gy*gyz+gz*gzz) end do end do c return end c c solution of the linear system of 3 eqs : dgg . dxv = dff c input vectors : dxv1, dxv2, dxv3, dff c output vector : dgg c c dff=(1,0,0) c subroutine solg0(dxv1,dxv2,dxv3,dgg) double precision denom,aa1,aa2,aa3 double precision dxv1(3),dxv2(3),dxv3(3),dgg(3) aa1=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3)) aa2=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2)) aa3=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2)) denom = dxv1(1)*aa1-dxv2(1)*aa2+dxv3(1)*aa3 dgg(1) = aa1/denom dgg(2) = -(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3))/denom dgg(3) = (dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2))/denom return end c c three rhs vectors df1, df2, df3 c subroutine solg3(dxv1,dxv2,dxv3,df1,df2,df3,dg1,dg2,dg3) double precision denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 double precision dxv1(3),dxv2(3),dxv3(3) double precision df1(3),df2(3),df3(3) double precision dg1(3),dg2(3),dg3(3) a11=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3)) a21=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2)) a31=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2)) a12=(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3)) a22=(dxv1(1)*dxv3(3)-dxv1(3)*dxv3(1)) a32=(dxv1(1)*dxv2(3)-dxv1(3)*dxv2(1)) a13=(dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2)) a23=(dxv1(1)*dxv3(2)-dxv1(2)*dxv3(1)) a33=(dxv1(1)*dxv2(2)-dxv1(2)*dxv2(1)) denom = dxv1(1)*a11-dxv2(1)*a21+dxv3(1)*a31 dg1(1) = (df1(1)*a11-df1(2)*a21+df1(3)*a31)/denom dg1(2) = -(df1(1)*a12-df1(2)*a22+df1(3)*a32)/denom dg1(3) = (df1(1)*a13-df1(2)*a23+df1(3)*a33)/denom dg2(1) = (df2(1)*a11-df2(2)*a21+df2(3)*a31)/denom dg2(2) = -(df2(1)*a12-df2(2)*a22+df2(3)*a32)/denom dg2(3) = (df2(1)*a13-df2(2)*a23+df2(3)*a33)/denom dg3(1) = (df3(1)*a11-df3(2)*a21+df3(3)*a31)/denom dg3(2) = -(df3(1)*a12-df3(2)*a22+df3(3)*a32)/denom dg3(3) = (df3(1)*a13-df3(2)*a23+df3(3)*a33)/denom return end c c Runge-Kutta integrator c subroutine rkint4 implicit real*8 (a-h,o-z) parameter(ndim=6,jmx=31,kmx=36) dimension y(ndim),yy(ndim) dimension fk1(ndim),fk2(ndim),fk3(ndim),fk4(ndim) dimension dgr2(3),dgr(3),ddgr(3,3) dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) c common/nray/nrayr,nrayth common/dsds/dst c common/wrk/ywrk,ypwrk common/grad2jk/grad2 common/dgrad2jk/dgrad2v common/gradjk/gri common/ggradjk/ggri common/gr/gr2 common/dgr/dgr2,dgr,ddgr common/igrad/igrad c h=dst hh=h*0.5d0 h6=h/6.0d0 c do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk do ieq=1,ndim y(ieq)=ywrk(ieq,j,k) fk1(ieq)=ypwrk(ieq,j,k) yy(ieq)=y(ieq)+fk1(ieq)*hh end do gr2=grad2(j,k) do iv=1,3 dgr2(iv)=dgrad2v(iv,j,k) dgr(iv)=gri(iv,j,k) do jv=1,3 ddgr(iv,jv)=ggri(iv,jv,j,k) end do end do call fwork(yy,fk2) c do ieq=1,ndim yy(ieq)=y(ieq)+fk2(ieq)*hh end do call fwork(yy,fk3) c do ieq=1,ndim yy(ieq)=y(ieq)+fk3(ieq)*h end do call fwork(yy,fk4) c do ieq=1,ndim ywrk(ieq,j,k)=y(ieq) . +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq)) end do end do end do c call updatepos c if(igrad.eq.1) call gradi c return end c c c subroutine gwork(j,k) implicit real*8 (a-h,o-z) parameter(ndim=6,jmx=31,kmx=36) dimension yy(ndim),yyp(ndim) dimension dgr2(3),dgr(3),ddgr(3,3) dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) c common/igrad/igrad c common/wrk/ywrk,ypwrk common/grad2jk/grad2 common/dgrad2jk/dgrad2v common/gradjk/gri common/ggradjk/ggri common/gr/gr2 common/dgr/dgr2,dgr,ddgr c c begin --- update vector yy c do ieq=1,ndim yy(ieq)=ywrk(ieq,j,k) end do c if(igrad.eq.1) then gr2=grad2(j,k) do iv=1,3 dgr2(iv)=dgrad2v(iv,j,k) dgr(iv)=gri(iv,j,k) do jv=1,3 ddgr(iv,jv)=ggri(iv,jv,j,k) end do end do end if c call fwork(yy,yyp) c do ieq=1,ndim ypwrk(ieq,j,k)=yyp(ieq) end do c c end --- update vector yy c return end c c c subroutine fwork(y,dery) implicit real*8 (a-h,o-z) parameter(ndim=6) dimension y(ndim),dery(ndim) dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3) dimension derdxv(3),danpldxv(3),derdnv(3) dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3) c common/gr/gr2 common/dgr/dgr2,dgr,ddgr common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/mode/sox common/nplr/anpl,anpr common/bb/bv common/dbb/derbv common/igrad/igrad common/xgxg/xg common/ygyg/yg common/dxgyg/derxg,deryg common/vgv/vgm,derdnm common/dersdst/dersdst common/ierr/ierr common/anv/anv common/xv/xv common/idst/idst c xx=y(1) yy=y(2) zz=y(3) xv(1)=y(1) xv(2)=y(2) xv(3)=y(3) c anv(1) = y(4) anv(2) = y(5) anv(3) = y(6) c C rr=sqrt(xx**2+yy**2) C phi=acos(xx/rr) C if (yy.lt.0.0d0) then C phi=-phi C end if c call plas_deriv(xx,yy,zz) c an2=anv(1)*anv(1)+anv(2)*anv(2)+anv(3)*anv(3) anpl=anv(1)*bv(1)+anv(2)*bv(2)+anv(3)*bv(3) c if(abs(anpl).gt.0.99d0) then if(abs(anpl).le.1.05d0) then ierr=97 else ierr=98 end if end if c anpl2=anpl*anpl dnl=1.0d0-anpl2 anpr2=an2-anpl2 anpr=0.0d0 if(anpr2.gt.0.0d0) anpr=sqrt(anpr2) yg2=yg**2 c an2s=1.0d0 dan2sdxg=0.0d0 dan2sdyg=0.0d0 dan2sdnpl=0.0d0 del=0.0d0 fdia=0.0d0 dfdiadnpl=0.0d0 dfdiadxg=0.0d0 dfdiadyg=0.0d0 c duh=1.0d0-xg-yg2 if(xg.gt.0.0d0) then del=sqrt(dnl*dnl+4.0d0*anpl*anpl*(1.0d0-xg)/yg2) an2s=1.0d0-xg -xg*yg2*(1.0d0+anpl2+sox*del)/duh/2.0d0 c dan2sdxg=-yg2*(1.0d0-yg2)*(1.0d0+anpl2+sox*del)/duh**2/2.0d0 . +sox*xg*anpl2/(del*duh)-1.0d0 dan2sdyg=-xg*yg*(1.0d0-xg)*(1.0d0+anpl2+sox*del)/duh**2 . +2.0d0*sox*xg*(1.0d0-xg)*anpl2/(yg*del*duh) dan2sdnpl=-xg*yg2*anpl/duh . -sox*xg*anpl*(2.0d0*(1.0d0-xg)-yg2*dnl)/(del*duh) c if(igrad.gt.0) then ddelnpl2=2.0d0*(2.0d0*(1.0d0-xg)*(1.0d0+3.0d0*anpl2**2) . -yg2*dnl**3)/yg2/del**3 fdia=-xg*yg2*(1.0d0+sox*ddelnpl2/2.0d0)/duh derdel=2.0d0*(1.0d0-xg)*anpl2*(1.0d0+3.0d0*anpl2**2) . -dnl**2*(1.0d0+3.0d0*anpl2)*yg2 derdel=4.0d0*derdel/(yg**5*del**5) ddelnpl2y=2.0d0*(1.0d0-xg)*derdel ddelnpl2x=yg*derdel dfdiadnpl=24.0d0*sox*xg*(1.0d0-xg)*anpl*(1.0d0-anpl**4) . /(yg2*del**5) dfdiadxg=-yg2*(1.0d0-yg2)/duh**2-sox*yg2*((1.0d0-yg2)*ddelnpl2 . +xg*duh*ddelnpl2x)/(2.0d0*duh**2) dfdiadyg=-2.0d0*yg*xg*(1.0d0-xg)/duh**2 . -sox*xg*yg*(2.0d0*(1.0d0-xg)*ddelnpl2 . +yg*duh*ddelnpl2y)/(2.0d0*duh**2) c end if end if c bdotgr=0.0d0 do iv=1,3 bdotgr=bdotgr+bv(iv)*dgr(iv) dbgr(iv)=0.0d0 do jv=1,3 dbgr(iv)=dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv) end do end do c derdnm=0.0d0 c do iv=1,3 danpldxv(iv)=anv(1)*derbv(1,iv)+anv(2)*derbv(2,iv) . +anv(3)*derbv(3,iv) derdxv(iv)=-(derxg(iv)*dan2sdxg . +deryg(iv)*dan2sdyg+danpldxv(iv)*dan2sdnpl) derdxv(iv)=derdxv(iv)-igrad*dgr2(iv) c derdxv(iv)=derdxv(iv)+fdia*bdotgr*dbgr(iv)+0.5d0*bdotgr**2* . (derxg(iv)*dfdiadxg+deryg(iv)*dfdiadyg+danpldxv(iv)*dfdiadnpl) c derdnv(iv)=2.0d0*anv(iv)-bv(iv)*dan2sdnpl . +0.5d0*bdotgr**2*bv(iv)*dfdiadnpl c derdnm=derdnm+derdnv(iv)**2 end do c derdnm=sqrt(derdnm) c derdom=-2.0d0*an2+2.0d0*xg*dan2sdxg+yg*dan2sdyg+anpl*dan2sdnpl . +2.0d0*igrad*gr2 . -bdotgr**2*(fdia+xg*dfdiadxg+yg*dfdiadyg/2.0d0 . +anpl*dfdiadnpl/2.0d0) c if (idst.eq.0) then c integration variable: s denom=-derdnm else if (idst.eq.1) then c integration variable: c*t denom=derdom else c integration variable: Sr denom=-(anv(1)*derdnv(1)+anv(2)*derdnv(2)+anv(3)*derdnv(3)) end if c c coefficient for integration in s c ds/dst, where st is the integration variable dersdst=-derdnm/denom c dery(1) = -derdnv(1)/denom dery(2) = -derdnv(2)/denom dery(3) = -derdnv(3)/denom dery(4) = derdxv(1)/denom dery(5) = derdxv(2)/denom dery(6) = derdxv(3)/denom c c vgv : ~ group velocity c vgm=0 do iv=1,3 vgv(iv)=-derdnv(iv)/derdom vgm=vgm+vgv(iv)**2 end do vgm=sqrt(vgm) c c dd : dispersion relation (real part) c ddi : dispersion relation (imaginary part) c dd=an2-an2s-igrad*(gr2-0.5d0*bdotgr**2*fdia) ddi=derdnv(1)*dgr(1)+derdnv(2)*dgr(2)+derdnv(3)*dgr(3) c return end c c c subroutine plas_deriv(xx,yy,zz) implicit real*8 (a-h,o-z) parameter(pi=3.14159265358979d0) dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3) dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3) c common/parbres/bres common/parpl/brr,bphi,bzz,ajphi common/btot/btot common/bb/bv common/dbb/derbv common/xgxg/xg common/dxgdps/dxgdpsi common/ygyg/yg common/dxgyg/derxg,deryg common/iieq/iequil common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/psival/psinv common/sgnib/sgnbphi,sgniphi c xg=0.0d0 yg=9.9d1 c do iv=1,3 derxg(iv)=0.0d0 deryg(iv)=0.0d0 bv(iv)=0.0d0 dbtot(iv)=0.0d0 do jv=1,3 dbv(iv,jv)=0.0d0 derbv(iv,jv)=0.0d0 dbvcdc(iv,jv)=0.0d0 dbvcdc(iv,jv)=0.0d0 dbvdc(iv,jv)=0.0d0 end do end do c if(iequil.eq.0) return c c cylindrical coordinates c rr2=xx**2+yy**2 rr=sqrt(rr2) csphi=xx/rr snphi=yy/rr c bv(1)=-snphi*sgnbphi bv(2)=csphi*sgnbphi c c convert from cm to meters c zzm=1.0d-2*zz rrm=1.0d-2*rr c if(iequil.eq.1) then call equian(rrm,zzm) yg=fpolv/(rrm*bres) end if c if(iequil.eq.2) then call equinum(rrm,zzm) call sub_xg_derxg yg=fpolv/(rrm*bres) bphi=fpolv/rrm btot=abs(bphi) if(psinv.lt.0.0d0) return end if c c B = f(psi)/R e_phi+ grad psi x e_phi/R c c bvc(i) = B_i in cylindrical coordinates c bv(i) = B_i in cartesian coordinates c bphi=fpolv/rrm brr=-dpsidz/rrm bzz= dpsidr/rrm c dfpolv=ffpv/fpolv c bvc(1)=brr bvc(2)=bphi bvc(3)=bzz c bv(1)=bvc(1)*csphi-bvc(2)*snphi bv(2)=bvc(1)*snphi+bvc(2)*csphi bv(3)=bvc(3) c b2tot=bv(1)**2+bv(2)**2+bv(3)**2 btot=sqrt(b2tot) c c dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) c dbvcdc(1,1)=-ddpsidrz/rrm-brr/rrm dbvcdc(1,3)=-ddpsidzz/rrm dbvcdc(2,1)= dfpolv*dpsidr/rrm-bphi/rrm dbvcdc(2,3)= dfpolv*dpsidz/rrm dbvcdc(3,1)= ddpsidrr/rrm-bzz/rrm dbvcdc(3,3)= ddpsidrz/rrm c c dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) c dbvdc(1,1) = dbvcdc(1,1)*csphi-dbvcdc(2,1)*snphi dbvdc(1,2) = -bv(2) dbvdc(1,3) = dbvcdc(1,3)*csphi-dbvcdc(2,3)*snphi dbvdc(2,1) = dbvcdc(1,1)*snphi+dbvcdc(2,1)*csphi dbvdc(2,2) = bv(1) dbvdc(2,3) = dbvcdc(1,3)*snphi+dbvcdc(2,3)*csphi dbvdc(3,1) = dbvcdc(3,1) dbvdc(3,2) = dbvcdc(3,2) dbvdc(3,3) = dbvcdc(3,3) c drrdx=csphi drrdy=snphi dphidx=-snphi/rrm dphidy=csphi/rrm c c dbv(iv,jv) = d Bcart(iv) / dxvcart(jv) c dbv(1,1)=drrdx*dbvdc(1,1)+dphidx*dbvdc(1,2) dbv(1,2)=drrdy*dbvdc(1,1)+dphidy*dbvdc(1,2) dbv(1,3)=dbvdc(1,3) dbv(2,1)=drrdx*dbvdc(2,1)+dphidx*dbvdc(2,2) dbv(2,2)=drrdy*dbvdc(2,1)+dphidy*dbvdc(2,2) dbv(2,3)=dbvdc(2,3) dbv(3,1)=drrdx*dbvdc(3,1)+dphidx*dbvdc(3,2) dbv(3,2)=drrdy*dbvdc(3,1)+dphidy*dbvdc(3,2) dbv(3,3)=dbvdc(3,3) c dbtot(1)=(bv(1)*dbv(1,1)+bv(2)*dbv(2,1)+bv(3)*dbv(3,1))/btot dbtot(2)=(bv(1)*dbv(1,2)+bv(2)*dbv(2,2)+bv(3)*dbv(3,2))/btot dbtot(3)=(bv(1)*dbv(1,3)+bv(2)*dbv(2,3)+bv(3)*dbv(3,3))/btot c yg=btot/Bres c c convert spatial derivatives from dummy/m -> dummy/cm c to be used in fwork c c bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z c do iv=1,3 deryg(iv)=1.0d-2*dbtot(iv)/Bres bv(iv)=bv(iv)/btot do jv=1,3 derbv(iv,jv)=1.0d-2*(dbv(iv,jv)-bv(iv)*dbtot(jv))/btot end do end do c derxg(1)=1.0d-2*drrdx*dpsidr*dxgdpsi derxg(2)=1.0d-2*drrdy*dpsidr*dxgdpsi derxg(3)=1.0d-2*dpsidz*dxgdpsi c c current density computation in Ampere/m^2 c ccj==1/mu_0 c ccj=1.0d+7/(4.0d0*pi) ajphi=ccj*(dbvcdc(1,3)-dbvcdc(3,1)) c ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3)) c ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2)) c return end c c c subroutine equian(rrm,zzm) implicit real*8 (a-h,o-z) c common/parqq/q0,qa,alq common/parban/b0,rr0m,zr0m,rpam common/psival/psinv common/pareq1/psia common/densbnd/psdbnd common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/xgxg/xg common/dxgdps/dxgdpsi common/xgcn/xgcn common/dens/dens,ddens common/sgnib/sgnbphi,sgniphi common/bmxmn/bmxi,bmni common/fc/fci common/iipr/iprof c if(iprof.eq.0) psdbnd=1.0d0 c dpsidrp=0.0d0 d2psidrp=0.0d0 c c simple model for equilibrium c rpm=sqrt((rrm-rr0m)**2+(zzm-zr0m)**2) rn=rpm/rpam c snt=0.0d0 cst=1.0d0 if (rpm.gt.0.0d0) then snt=(zzm-zr0m)/rpm cst=(rrm-rr0m)/rpm end if c c valore fittizio di psinv e di psia: psinv=rn**2 psia=sgniphi c sgn=sgniphi*sgnbphi if(rn.le.1.0d0) then qq=q0+(qa-q0)*rn**alq dpsidrp=B0*rpam*rn/qq*sgn dqq=alq*(qa-q0)*rn**(alq-1.0d0) d2psidrp=B0*(1.0d0-rn*dqq/qq)/qq*sgn else dpsidrp=B0*rpam/qa/rn*sgn d2psidrp=-B0/qa/rn**2 *sgn end if c fpolv=sgnbphi*b0*rr0m dfpolv=0.0d0 ffpv=sgniphi*fpolv*dfpolv c dpsidr=dpsidrp*cst dpsidz=dpsidrp*snt ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm) ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp c dens=0.0d0 xg=0.0d0 dxgdpsi=0.0d0 if(psinv.lt.psdbnd) then call density(psinv) end if xg=xgcn*dens ddensdrp=2.0d0*rn*ddens if(dpsidrp.ne.0.0d0) dxgdpsi=xgcn*ddensdrp/(dpsidrp*rpam) c bmax=sqrt(fpolv**2+dpsidrp**2)/(rr0m-rpam*rn) bmin=sqrt(fpolv**2+dpsidrp**2)/(rr0m+rpam*rn) c return end c c c subroutine equinum(rpsim,zpsim) implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) parameter(nrest=nnw+4,nzest=nnh+4) parameter(lwrk=8,liwrk=2) dimension ccspl(nnw*nnh) dimension tr(nrest),tz(nzest),wrk(lwrk),iwrk(liwrk) dimension rrs(1),zzs(1),ffspl(1) parameter(lw10=nnw*3+nnh*4+nnw*nnh,lw01=nnw*4+nnh*3+nnw*nnh) parameter(lw20=nnw*2+nnh*4+nnw*nnh,lw02=nnw*4+nnh*2+nnw*nnh) parameter(lw11=nnw*3+nnh*3+nnw*nnh) parameter(nrs=1,nzs=1) dimension cc01(lw10),cc10(lw01),cc02(lw02),cc20(lw20),cc11(lw11) dimension tfp(nrest),cfp(nrest),wrkfd(nrest) c common/eqnn/nr,nz,npp,nintp common/psival/psinv common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv common/pareq1/psia common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/pareq1t/psiant,psinop c common/coffeqt/tr,tz common/coffeqtp/tfp common/coffeq/ccspl common/coffeqd/cc01,cc10,cc20,cc02,cc11 common/coffeqn/nsrt,nszt,nsft common/cofffp/cfp common/fpas/fpolas c psinv=-1.0d0 c dpsidr=0.0d0 dpsidz=0.0d0 ddpsidrr=0.0d0 ddpsidzz=0.0d0 ddpsidrz=0.0d0 c c here lengths are measured in meters c fpolv=fpolas ffpv=0.0d0 c if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return c rrs(1)=rpsim zzs(1)=zpsim nsr=nsrt nsz=nszt 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 nur=1 nuz=0 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc10,kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc10(iwr),cc10(iwz),iwrk(1),iwrk(2)) dpsidr= ffspl(1)*psia c nur=0 nuz=1 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc01,kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc01(iwr),cc01(iwz),iwrk(1),iwrk(2)) dpsidz= ffspl(1)*psia c nur=2 nuz=0 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc20,kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc20(iwr),cc20(iwz),iwrk(1),iwrk(2)) ddpsidrr= ffspl(1)*psia c nur=0 nuz=2 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc02,kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc02(iwr),cc02(iwz),iwrk(1),iwrk(2)) ddpsidzz= ffspl(1)*psia c nur=1 nuz=1 kkr=3-nur kkz=3-nuz iwr=1+(nr-nur-4)*(nz-nuz-4) iwz=iwr+4-nur call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc11,kkr,kkz, . rrs,nrs,zzs,nzs,ffspl,cc11(iwr),cc11(iwz),iwrk(1),iwrk(2)) ddpsidrz= ffspl(1)*psia c if(psinv.le.1.0d0.and.psinv.gt.0.0d0) then rrs(1)=psinv call splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier) fpolv=ffspl(1) nu=1 call splder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier) dfpolv=ffspl(1) ffpv=fpolv*dfpolv/psia end if c return end c c c subroutine bfield(rpsim,zpsim,bphi,brr,bzz) implicit real*8 (a-h,o-z) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv call equinum(rpsim,zpsim) bphi=fpolv/rpsim brr=-dpsidz/rpsim bzz= dpsidr/rpsim return end c subroutine tor_curr_psi(h,ajphi) implicit real*8 (a-h,o-z) common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz call psi_raxis(h,r1,r2) call tor_curr(r2,zmaxis,ajphi) return end c subroutine tor_curr(rpsim,zpsim,ajphi) implicit real*8 (a-h,o-z) parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi)) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv call equinum(rpsim,zpsim) bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim ajphi=ccj*(dbvcdc13-dbvcdc31) return end c subroutine psi_raxis(h,r1,r2) implicit real*8 (a-h,o-z) parameter(mest=4,kspl=3) parameter(nnw=501,nnh=501) parameter(nrest=nnw+4,nzest=nnh+4) dimension cc(nnw*nnh),tr(nrest),tz(nzest) dimension czc(nrest),zeroc(mest) c common/pareq1t/psiant,psinop common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz common/coffeqn/nsr,nsz,nsft common/coffeq/cc common/coffeqt/tr,tz c iopt=1 zc=zmaxis call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) print*,' profil =',ier val=h*psiant+psinop call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) r1=zeroc(1) r2=zeroc(2) return end c c c subroutine sub_xg_derxg implicit real*8 (a-h,o-z) common/psival/psinv common/pareq1/psia common/xgxg/xg common/dxgdps/dxgdpsi common/xgcn/xgcn common/dens/dens,ddenspsin xg=0.0d0 dxgdpsi=0.0d0 c if(psinv.le.psdbnd.and.psinv.ge.0) then call density(psinv) xg=xgcn*dens dxgdpsi=xgcn*ddenspsin/psia c end if return end c c c subroutine density(arg) implicit real*8 (a-h,o-z) parameter(npmx=250,npest=npmx+4) dimension xxs(1),ffs(1) dimension tfn(npest),cfn(npest),wrkfd(npest) c common/densbnd/psdbnd common/denspp/psnpp,aad,bbd,ccd common/iipr/iprof common/pardens/dens0,aln1,aln2 common/dens/dens,ddens common/coffdt/tfn common/coffdnst/nsfd common/cofffn/cfn c c computation of density [10^19 m^-3] and derivative wrt psi c dens=0.0d0 ddens=0.0d0 if(arg.gt.psdbnd.or.arg.lt.0.0d0) return c if(iprof.eq.0) then if(arg.gt.1.0d0) return profd=(1.0d0-arg**aln1)**aln2 dens=dens0*profd dprofd=-aln1*aln2*arg**(aln1-1.0d0) . *(1.0d0-arg**aln1)**(aln2-1.0d0) ddens=dens0*dprofd else if(arg.le.psdbnd.and.arg.gt.psnpp) then c c cubic interpolation for 1 < psi < psdbnd c nn=3 nn1=nn+1 nn2=nn+2 dpsib=arg-psdbnd dens=aad*dpsib**nn+bbd*dpsib**nn1+ccd*dpsib**nn2 ddens=nn*aad*dpsib**(nn-1)+nn1*bbd*dpsib**nn . +nn2*ccd*dpsib**nn1 else xxs(1)=arg ier=0 call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier) dens=ffs(1) nu=1 ier=0 call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier) ddens=ffs(1) if(ier.gt.0) print*,ier if(abs(dens).lt.1.0d-10) dens=0.0d0 end if if(dens.lt.0.0d0) print*,' DENSITY NEGATIVE',dens c end if return end c c c function temp(arg) implicit real*8 (a-h,o-z) parameter(npmx=250) dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),ct(npmx,4) c common/parqte/te0,dte0,alt1,alt2 common/iipr/iprof common/crad/psrad,derad,terad,zfc common/coffte/ct common/eqnn/nr,nz,npp,nintp c temp=0.0d0 if(arg.ge.1.0d0.or.arg.lt.0.0d0) return if(iprof.eq.0) then proft=(1.0d0-arg**alt1)**alt2 temp=(te0-dte0)*proft+dte0 else call locate(psrad,npp,arg,k) if(k.eq.0) k=1 if(k.eq.npp) k=npp-1 dps=arg-psrad(k) temp=spli(ct,npmx,k,dps) endif return end c c c function fzeff(arg) implicit real*8 (a-h,o-z) parameter(npmx=250) dimension psrad(npmx),derad(npmx),terad(npmx),zfc(npmx),cz(npmx,4) c common/iipr/iprof common/crad/psrad,derad,terad,zfc common/coffz/cz common/eqnn/nr,nz,npp,nintp common/zz/Zeff c fzeff=1 ps=arg if(ps.gt.1.0d0.and.ps.lt.0.0d0) return if(iprof.eq.0) then fzeff=zeff else call locate(psrad,npp,ps,k) if(k.eq.0) k=1 if(k.eq.npp) k=npp-1 dps=ps-psrad(k) fzeff=spli(cz,npmx,k,dps) endif return end c c beam tracing initial conditions igrad=1 c subroutine ic_gb implicit real*8 (a-h,o-z) 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) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) 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 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 complex*16 catand external catand c parameter(ui=(0.0d0,1.0d0)) c common/nray/nrayr,nrayth common/rwmax/rwmax common/parwv/ak0,akinv,fhz common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/wrk/ywrk0,ypwrk0 common/grc/xc0,du10 common/fi/dffiu,ddffiu common/mirr/x00,y00,z00 common/grad2jk/grad2 common/dgrad2jk/dgrad2v common/gradjk/gri common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 c ui=(0.0d0,1.0d0) csth=anz0c snth=sqrt(1.0d0-csth**2) csps=1.0d0 snps=0.0d0 if(snth.gt.0.0d0) then csps=any0c/snth snps=anx0c/snth end if c phiwrad=phiw*cvdr phirrad=phir*cvdr csphiw=cos(phiwrad) snphiw=sin(phiwrad) c csphir=cos(phirrad) c snphir=sin(phirrad) c wwcsi=2.0d0*akinv/wcsi**2 wweta=2.0d0*akinv/weta**2 c if(phir.ne.phiw) then sk=(rcicsi+rcieta) sw=(wwcsi+wweta) dk=(rcicsi-rcieta) dw=(wwcsi-wweta) ts=-(dk*sin(2.0d0*phirrad)-ui*dw*sin(2.0d0*phiwrad)) tc=(dk*cos(2.0d0*phirrad)-ui*dw*cos(2.0d0*phiwrad)) phic=0.5d0*catand(ts/tc) ddd=dk*cos(2*(phirrad+phic))-ui*dw*cos(2*(phiwrad+phic)) sss=sk-ui*sw qi1=0.5d0*(sss+ddd) qi2=0.5d0*(sss-ddd) rci1=dble(qi1) ww1=-dimag(qi1) rci2=dble(qi2) ww2=-dimag(qi2) else rci1=rcicsi rci2=rcieta ww1=wwcsi ww2=wweta phic=-phiwrad qi1=rci1-ui*ww1 qi2=rci2-ui*ww2 end if c w01=sqrt(2.0d0*akinv/ww1) c z01=-rci1/(rci1**2+ww1**2) c w02=sqrt(2.0d0*akinv/ww2) c z02=-rci2/(rci2**2+ww2**2) qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2 qqyy=qi1*sin(phic)**2+qi2*cos(phic)**2 qqxy=-(qi1-qi2)*sin(2.0d0*phic) wwxx=-dimag(qqxx) wwyy=-dimag(qqyy) wwxy=-dimag(qqxy)/2.0d0 rcixx=dble(qqxx) rciyy=dble(qqyy) rcixy=dble(qqxy)/2.0d0 dqi1=-qi1**2 dqi2=-qi2**2 d2qi1=2*qi1**3 d2qi2=2*qi2**3 dqqxx=dqi1*cos(phic)**2+dqi2*sin(phic)**2 dqqyy=dqi1*sin(phic)**2+dqi2*cos(phic)**2 dqqxy=-(dqi1-dqi2)*sin(2.0d0*phic) d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2 d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2 d2qqxy=-(d2qi1-d2qi2)*sin(2.0d0*phic) dwwxx=-dimag(dqqxx) dwwyy=-dimag(dqqyy) dwwxy=-dimag(dqqxy)/2.0d0 d2wwxx=-dimag(d2qqxx) d2wwyy=-dimag(d2qqyy) d2wwxy=-dimag(d2qqxy)/2.0d0 drcixx=dble(dqqxx) drciyy=dble(dqqyy) drcixy=dble(dqqxy)/2.0d0 dr=1.0d0 if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) da=2.0d0*pi/dble(nrayth) c ddfu=2.0d0*dr**2*akinv do j=1,nrayr u=dble(j-1) c ffi=u**2*ddfu/2.0d0 dffiu(j)=u*ddfu ddffiu(j)=ddfu do k=1,nrayth alfak=(k-1)*da dcsiw=dr*cos(alfak)*wcsi detaw=dr*sin(alfak)*weta dx0t=dcsiw*csphiw-detaw*snphiw dy0t=dcsiw*snphiw+detaw*csphiw x0t=u*dx0t y0t=u*dy0t z0t=-0.5d0*(rcixx*x0t**2+rciyy*y0t**2+2*rcixy*x0t*y0t) c c csiw=u*dcsiw c etaw=u*detaw c csir=x0t*csphir+y0t*snphir c etar=-x0t*snphir+y0t*csphir c dx0= x0t*csps+snps*(y0t*csth+z0t*snth) dy0=-x0t*snps+csps*(y0t*csth+z0t*snth) dz0= z0t*csth-y0t*snth x0=x00+dx0 y0=y00+dy0 z0=z00+dz0 gxt=x0t*wwxx+y0t*wwxy gyt=x0t*wwxy+y0t*wwyy gzt=0.5d0*(x0t**2*dwwxx+y0t**2*dwwyy)+x0t*y0t*dwwxy gr2=gxt*gxt+gyt*gyt+gzt*gzt gxxt=wwxx gyyt=wwyy gzzt=0.5d0*(x0t**2*d2wwxx+y0t**2*d2wwyy)+x0t*y0t*d2wwxy gxyt=wwxy gxzt=x0t*dwwxx+y0t*dwwxy gyzt=x0t*dwwxy+y0t*dwwyy dgr2xt=2.0d0*(gxt*gxxt+gyt*gxyt+gzt*gxzt) dgr2yt=2.0d0*(gxt*gxyt+gyt*gyyt+gzt*gyzt) dgr2zt=2.0d0*(gxt*gxzt+gyt*gyzt+gzt*gzzt) dgr2x= dgr2xt*csps+snps*(dgr2yt*csth+dgr2zt*snth) dgr2y=-dgr2xt*snps+csps*(dgr2yt*csth+dgr2zt*snth) dgr2z= dgr2zt*csth-dgr2yt*snth gri(1,j,k)=gxt*csps+snps*(gyt*csth+gzt*snth) gri(2,j,k)=-gxt*snps+csps*(gyt*csth+gzt*snth) gri(3,j,k)=gzt*csth-gyt*snth ggri(1,1,j,k)=gxxt*csps**2+ . snps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)+ . 2.0d0*snps*csps*(gxyt*csth+gxzt*snth) ggri(2,2,j,k)=gxxt*snps**2+ . csps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)- . 2.0d0*snps*csps*(gxyt*csth+gxzt*snth) ggri(3,3,j,k)=gzzt*csth**2+gyyt*snth**2-2.0d0*csth*snth*gyzt ggri(1,2,j,k)=csps*snps*(-gxxt+csth**2*gyyt+snth**2*gzzt . +2.0d0*csth*snth*gyzt) . +(csps**2-snps**2)*(snth*gxzt+csth*gxyt) ggri(1,3,j,k)=csth*snth*snps*(gzzt-gyyt) . +(csth**2-snth**2)*snps*gyzt+csps*(csth*gxzt-snth*gxyt) ggri(2,3,j,k)=csth*snth*csps*(gzzt-gyyt) . +(csth**2-snth**2)*csps*gyzt+snps*(snth*gxyt-csth*gxzt) ggri(2,1,j,k)=ggri(1,2,j,k) ggri(3,1,j,k)=ggri(1,3,j,k) ggri(3,2,j,k)=ggri(2,3,j,k) c du1tx=(dx0t*wwxx+dy0t*wwxy)/ddfu du1ty=(dx0t*wwxy+dy0t*wwyy)/ddfu du1tz=0.5d0*u*(dx0t**2*dwwxx+dy0t**2*dwwyy+ . 2.0d0*dx0t*dy0t*dwwxy)/ddfu c pppx=x0t*rcixx+y0t*rcixy pppy=x0t*rcixy+y0t*rciyy denpp=pppx*gxt+pppy*gyt if (denpp.ne.0.0d0) then ppx=-pppx*gzt/denpp ppy=-pppy*gzt/denpp else ppx=0.0d0 ppy=0.0d0 end if c anzt=sqrt((1.0d0+gr2)/(1.0d0+ppx**2+ppy**2)) anxt=ppx*anzt anyt=ppy*anzt c anx= anxt*csps+snps*(anyt*csth+anzt*snth) any=-anxt*snps+csps*(anyt*csth+anzt*snth) anz= anzt*csth-anyt*snth c an20=1.0d0+gr2 an0=sqrt(an20) anx0=anx any0=any anz0=anz c xc0(1,j,k)=x0 xc0(2,j,k)=y0 xc0(3,j,k)=z0 c ywrk0(1,j,k)=x0 ywrk0(2,j,k)=y0 ywrk0(3,j,k)=z0 ywrk0(4,j,k)=anx0 ywrk0(5,j,k)=any0 ywrk0(6,j,k)=anz0 c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 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 grad2(j,k)=gr2 dgrad2v(1,j,k)=dgr2x dgrad2v(2,j,k)=dgr2y dgrad2v(3,j,k)=dgr2z c du10(1,j,k)= du1tx*csps+snps*(du1ty*csth+du1tz*snth) du10(2,j,k)=-du1tx*snps+csps*(du1ty*csth+du1tz*snth) du10(3,j,k)= du1tz*csth-du1ty*snth c dd=anx0**2+any0**2+anz0**2-an20 vgradi=anxt*gxt+anyt*gyt+anzt*gzt ddi=2.0d0*vgradi c 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 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 ddr110=dd end if if(j.eq.nrayr.and.k.eq.1) then write(17,99) zero,ddr110,dd,ddi end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end c c ray tracing initial conditions igrad=0 c subroutine ic_rt implicit real*8 (a-h,o-z) 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) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) c common/nray/nrayr,nrayth common/rwmax/rwmax common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/wrk/ywrk0,ypwrk0 common/grc/xc0,du10 common/fi/dffiu,ddffiu common/mirr/x00,y00,z00 common/grad2jk/grad2 common/dgrad2jk/dgrad2v common/gradjk/gri common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 c csth=anz0c snth=sqrt(1.0d0-csth**2) csps=1.0d0 snps=0.0d0 if(snth.gt.0.0d0) then csps=any0c/snth snps=anx0c/snth end if c phiwrad=phiw*cvdr csphiw=cos(phiwrad) snphiw=sin(phiwrad) c dr=1.0d0 if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) da=2.0d0*pi/dble(nrayth) z0t=0.0d0 c do j=1,nrayr u=dble(j-1) dffiu(j)=0.0d0 ddffiu(j)=0.0d0 do k=1,nrayth alfak=(k-1)*da dcsiw=dr*cos(alfak)*wcsi detaw=dr*sin(alfak)*weta dx0t=dcsiw*csphiw-detaw*snphiw dy0t=dcsiw*snphiw+detaw*csphiw x0t=u*dx0t y0t=u*dy0t c c csiw=u*dcsiw c etaw=u*detaw c csir=csiw c etar=etaw c dx0= x0t*csps+snps*(y0t*csth+z0t*snth) dy0=-x0t*snps+csps*(y0t*csth+z0t*snth) dz0= z0t*csth-y0t*snth c x0=x00+dx0 y0=y00+dy0 z0=z00+dz0 c ppcsi=u*dr*cos(alfak)*rcicsi ppeta=u*dr*sin(alfak)*rcieta c anzt=1.0d0/sqrt(1.0d0+ppcsi**2+ppeta**2) ancsi=ppcsi*anzt aneta=ppeta*anzt c anxt=ancsi*csphiw-aneta*snphiw anyt=ancsi*snphiw+aneta*csphiw c anx= anxt*csps+snps*(anyt*csth+anzt*snth) any=-anxt*snps+csps*(anyt*csth+anzt*snth) anz= anzt*csth-anyt*snth c an20=1.0d0 an0=sqrt(an20) anx0=anx any0=any anz0=anz c xc0(1,j,k)=x0 xc0(2,j,k)=y0 xc0(3,j,k)=z0 c ywrk0(1,j,k)=x0 ywrk0(2,j,k)=y0 ywrk0(3,j,k)=z0 ywrk0(4,j,k)=anx0 ywrk0(5,j,k)=any0 ywrk0(6,j,k)=anz0 c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 ypwrk0(4,j,k) = 0.0d0 ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 c do iv=1,3 gri(iv,j,k)=0.0d0 dgrad2v(iv,j,k)=0.0d0 du10(iv,j,k)=0.0d0 do jv=1,3 ggri(iv,jv,j,k)=0.0d0 end do end do grad2(j,k)=0.0d0 c dd=anx0**2+any0**2+anz0**2-an20 vgradi=0.0d0 ddi=2.0d0*vgradi c 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 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 end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end subroutine ic_rt2 implicit real*8 (a-h,o-z) 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) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension yyrfl(jmx,kmx,ndim) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) dimension dffiu(jmx),ddffiu(jmx) dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) c common/nray/nrayr,nrayth common/wrk/ywrk0,ypwrk0 common/grc/xc0,du10 common/grad2jk/grad2 common/dgrad2jk/dgrad2v common/gradjk/gri common/ggradjk/ggri common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/yyrfl/yyrfl do j=1,nrayr do k=1,nrayth x0=yyrfl(j,k,1) y0=yyrfl(j,k,2) z0=yyrfl(j,k,3) anx0=yyrfl(j,k,4) any0=yyrfl(j,k,5) anz0=yyrfl(j,k,6) an20=anx0*anx0+any0*any0+anz0*anz0 an0=sqrt(an20) c xc0(1,j,k)=x0 xc0(2,j,k)=y0 xc0(3,j,k)=z0 c ywrk0(1,j,k)=x0 ywrk0(2,j,k)=y0 ywrk0(3,j,k)=z0 ywrk0(4,j,k)=anx0/an0 ywrk0(5,j,k)=any0/an0 ywrk0(6,j,k)=anz0/an0 c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 ypwrk0(4,j,k) = 0.0d0 ypwrk0(5,j,k) = 0.0d0 ypwrk0(6,j,k) = 0.0d0 c do iv=1,3 gri(iv,j,k)=0.0d0 dgrad2v(iv,j,k)=0.0d0 du10(iv,j,k)=0.0d0 do jv=1,3 ggri(iv,jv,j,k)=0.0d0 end do end do grad2(j,k)=0.0d0 c dd=anx0**2+any0**2+anz0**2-an20 vgradi=0.0d0 ddi=2.0d0*vgradi c 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 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 end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end c c ray power weigth coefficient q(j) c subroutine pweigth implicit real*8(a-h,o-z) parameter(jmx=31) dimension q(jmx) common/qw/q common/nray/nrayr,nrayth common/rwmax/rwmax c dr=1.0d0 if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) r1=0.0d0 summ=0.0d0 q(1)=1.0d0 if(nrayr.gt.1) then do j=1,nrayr r2=(dble(j)-0.5d0)*dr q(j)=(exp(-2.0d0*r1**2)-exp(-2.0d0*r2**2)) r1=r2 summ=summ+q(j) end do c q(1)=q(1)/summ sm=q(1) j=1 k=1 do j=2,nrayr q(j)=q(j)/nrayth/summ do k=1,nrayth sm=sm+q(j) end do end do end if c return end c c c subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi, . bmxi,bmni,fci,intp) implicit real*8 (a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp),cbmx(nnintp,4),cbmn(nnintp,4) dimension cvol(nnintp,4),crri(nnintp,4),crbav(nnintp,4) dimension carea(nnintp,4),cfc(nnintp,4) c common/cflav/cvol,crri,crbav,cbmx,cbmn,carea,cfc 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 dps=rpsi-rpstab(ip) c areai=spli(carea,nintp,ip,dps) voli=spli(cvol,nintp,ip,dps) dervoli=splid(cvol,nintp,ip,dps) rrii=spli(crri,nintp,ip,dps) c if(intp.eq.0) return c rbavi=spli(crbav,nintp,ip,dps) bmxi=spli(cbmx,nintp,ip,dps) bmni=spli(cbmn,nintp,ip,dps) fci=spli(cfc,nintp,ip,dps) c return end c c c subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) implicit real*8 (a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp) dimension cratja(nnintp,4),cratjb(nnintp,4),cratjpl(nnintp,4) 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 dps=rpsi-rpstab(ip) ratjai=spli(cratja,nintp,ip,dps) ratjbi=spli(cratjb,nintp,ip,dps) ratjpli=spli(cratjpl,nintp,ip,dps) return end c c c subroutine pabs_curr(i,j,k) implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36,nmx=8000) parameter(pi=3.14159265358979d0) dimension psjki(jmx,kmx,nmx) 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 q(jmx),tau1v(jmx,kmx) c common/psjki/psjki common/atjki/tauv,alphav common/dpjjki/pdjki,currj common/pcjki/ppabs,ccci common/dcjki/didst common/tau1v/tau1v c common/qw/q common/p0/p0mw c common/dsds/dst common/dersdst/dersdst common/iieq/iequil c common/parban/b0,rr0m,zr0m,rpam common/absor/alpha,effjcd,akim,tau0 c common/psival/psinv common/sgnib/sgnbphi,sgniphi common/bmxmn/bmxi,bmni common/fc/fci c dvvl=1.0d0 rbavi=1.0d0 rrii=rr0m rhop0=sqrt(psjki(j,k,i-1)) if(iequil.eq.2) then intp=1 psinv=psjki(j,k,i) rhop=sqrt(psinv) call valpsispl(rhop,voli,dervoli,area,rrii, . rbavi,bmxi,bmni,fci,intp) dvvl=dervoli*abs(rhop-rhop0) else dvvl=4.0d0*pi*rr0m*pi*abs(rhop*rhop-rhop0*rhop0)*rpam**2 rbavi=b0/bmni fci=1.0d0-1.469d0*sqrt(rpam/rr0m) end if if(dvvl.le.0.0d0) dvvl=1.0d0 c adnm=2.0d0*pi*rrii c c calcolo della corrente cohen: currtot [MA] c calcolo della densita' corrente cohen: currj [MA/m^2] c calcolo della densita' potenza: pdjki [MW/m^3] c calcolo della efficienza : effjcdav [A m/W ] c tau0=tauv(j,k,i-1) c call ecrh_cd c alphav(j,k,i)=alpha dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0d0 tauv(j,k,i)=tauv(j,k,i-1)+dtau dpdst=p0mw*q(j)*exp(-tauv(j,k,i)-tau1v(j,k))* . alphav(j,k,i)*dersdst pdjki(j,k,i)=dpdst*dst/dvvl ppabs(j,k,i)=p0mw*q(j)*exp(-tau1v(j,k))* . (1.0d0-exp(-tauv(j,k,i))) effjcdav=rbavi*effjcd currj(j,k,i)=sgnbphi*effjcdav*pdjki(j,k,i) didst(j,k,i)=sgnbphi*effjcdav*dpdst/adnm dcidst=(didst(j,k,i)+didst(j,k,i-1))/2.0d0 ccci(j,k,i)=ccci(j,k,i-1)+dcidst*dst return end c c c subroutine ecrh_cd implicit real*8 (a-h,o-z) real*8 mc2,mesi parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8) parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) parameter(vcsi=2.99792458d+8) parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3) c common/ithn/ithn common/nharm/nharm,nhf common/warm/iwarm,ilarm common/ieccd/ieccd c common/ygyg/yg common/nplr/anpl,anpr common/vgv/vgm,derdnm common/parwv/ak0,akinv,fhz c common/nprw/anprre,anprim c common/absor/alpha,effjcd,akim,tau c common/psival/psinv common/tete/tekev common/amut/amu common/zz/Zeff c c absorption computation c alpha=0.0d0 akim=0.0d0 effjcd=0.0d0 c tekev=temp(psinv) amu=mc2/tekev zeff=fzeff(psinv) c if(tekev.le.0.0d0.or.tau.gt.taucr) return c dnl=1.0d0-anpl*anpl if(iwarm.eq.1) then ygc=1.0d0-anpl**2/2.0d0 else ygc=sqrt(1.0d0-anpl**2) end if c nharm=int(ygc/yg-eps)+1 c if(nharm.eq.0) return c nhf=0 do nn=nharm,nharm+4 ygn=nn*yg if(iwarm.eq.1) then rdu2=2.0d0*(ygn-ygc) u1=anpl+sqrt(rdu2) u2=anpl-sqrt(rdu2) uu2=min(u1*u1,u2*u2) argex=amu*uu2/2.0d0 else rdu2=ygn**2-ygc**2 g1=(ygn+anpl*sqrt(rdu2))/dnl g2=(ygn-anpl*sqrt(rdu2))/dnl gg=min(g1,g2) argex=amu*(gg-1.0d0) u1=(ygn*anpl+sqrt(rdu2))/dnl u2=(ygn*anpl-sqrt(rdu2))/dnl end if if(argex.le.xxcr) nhf=nn end do if(nhf.eq.0) return c lrm=ilarm if(lrm.lt.nhf) lrm=nhf call dispersion(lrm) c akim=ak0*anprim ratiovgr=2.0d0*anpr/derdnm c ratiovgr=2.0d0*anpr/derdnm*vgm alpha=2.0d0*akim*ratiovgr if(alpha.lt.0.0d0) then ierr=94 print*,' IERR = ', ierr,' alpha negative' end if c ithn=1 if(lrm.gt.nharm) ithn=2 if(ieccd.gt.0) call eccd(effjcd) return 999 format(12(1x,e12.5)) end c 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 complex*16 epsl(3,3,lrm),sepsl(3,3),e330 complex*16 e11,e22,e12,e13,e23 c complex*16 e33,e21,e31,e32 complex*16 a13,a31,a23,a32,a33 c common/ygyg/yg 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 c if (iwarm.eq.1) then call diel_tens_wr(e330,epsl,lrm) else call diel_tens_fr(e330,epsl,lrm) end if c do i=1,imx c do j=1,3 do k=1,3 sepsl(k,j)=dcmplx(0.0d0,0.0d0) do ilrm=1,lrm sepsl(k,j)=sepsl(k,j)+epsl(k,j,ilrm)*anpr2a**(ilrm-1) end do end do end do c anpra=sqrt(anpr2a) c e11=sepsl(1,1) e22=sepsl(2,2) e12=sepsl(1,2) a33=sepsl(3,3) a13=sepsl(1,3) a23=sepsl(2,3) a31=a13 a32=-a23 c e33=e330+anpr2a*a33 e13=anpra*a13 e23=anpra*a23 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) . -a32*e12*(a13+anpl)+a23*e12*(a31+anpl) . -(a23*a32+e330+(e22-anpl2)*(1.0d0-a33))*(e11-anpl2) . -(a13+anpl)*(a31+anpl)*(e22-anpl2) cc0=e330*((e11-anpl2)*(e22-anpl2)+e12*e12) c rr=cc2*cc2-4.0d0*cc0*cc4 c if(yg.gt.1.0d0) then s=sox if(dimag(rr).LE.0.0d0) s=-s else s=-sox if(dble(rr).le.0.d0.and.dimag(rr).ge.0.d0) s=-s 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)) anpr2a=anpr2 end do c 999 continue if(i.gt.imx) print*,' i>imx ',yg,errnpr,i 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) ez=dcmplx(0.0d0,0.0d0) c if (abs(anpl).gt.1.0d-6) then den=e12*e23-(e13+anpr*anpl)*(e22-anpr2-anpl2) ey=-(e12*(e13+anpr*anpl)+(e11-anpl2)*e23)/den ez=(e12*e12+(e22-anpr2-anpl2)*(e11-anpl2))/den ez2=abs(ez)**2 ey2=abs(ey)**2 enx2=1.0d0/(1.0d0+ey2+ez2) ex=dcmplx(sqrt(enx2),0.0d0) ez2=ez2*enx2 ey2=ey2*enx2 ez=ez*ex ey=ey*ex else if(sox.lt.0.0d0) then ez=dcmplx(1.0d0,0.0d0) ez2=abs(ez)**2 else ex2=1.0d0/(1.0d0+abs(-e11/e12)**2) ex=sqrt(ex2) ey=-ex*e11/e12 ey2=abs(ey)**2 ez2=0.0d0 end if end if c return end c c Fully relativistic case c computation of dielectric tensor elements c up to third order in Larmor radius for hermitian part c subroutine diel_tens_fr(e330,epsl,lrm) implicit real*8(a-h,o-z) complex*16 ui complex*16 e330,epsl(3,3,lrm) complex*16 ca11,ca12,ca22,ca13,ca23,ca33 complex*16 cq0p,cq0m,cq1p,cq1m,cq2p parameter(pi=3.14159265358979d0,rpi=1.77245385090552d0) parameter(ui=(0.0d0,1.0d0)) dimension rr(-lrm:lrm,0:2,0:lrm),ri(lrm,0:2,lrm) c common/xgxg/xg common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/resah/anpl2,dnl c common/cri/cr,ci common/warm/iwarm,ilarm c anpl2=anpl**2 dnl=1.0d0-anpl2 c cmxw=1.0d0+15.0d0/(8.0d0*amu)+105.0d0/(128.0d0*amu**2) cr=-amu*amu/(rpi*cmxw) ci=sqrt(2.0d0*pi*amu)*amu**2/cmxw c do l=1,lrm do j=1,3 do i=1,3 epsl(i,j,l)=dcmplx(0.0d0,0.0d0) end do end do end do c if(iwarm.eq.2) call hermitian(rr,lrm) if(iwarm.eq.4) call hermitian_2(rr,lrm) c call antihermitian(ri,lrm) c do l=1,lrm lm=l-1 fal=-0.25d0**l*fact(2*l)/(fact(l)**2*yg**(2*lm)) ca11=(0.0d0,0.0d0) ca12=(0.0d0,0.0d0) ca13=(0.0d0,0.0d0) ca22=(0.0d0,0.0d0) ca23=(0.0d0,0.0d0) ca33=(0.0d0,0.0d0) do is=0,l k=l-is w=(-1.0d0)**k c asl=w/(fact(is+l)*fact(l-is)) bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) c if(is.gt.0) then cq0p=rr(is,0,l)+rr(-is,0,l)+ui*ri(is,0,l) cq0m=rr(is,0,l)-rr(-is,0,l)+ui*ri(is,0,l) cq1p=rr(is,1,l)+rr(-is,1,l)+ui*ri(is,1,l) cq1m=rr(is,1,l)-rr(-is,1,l)+ui*ri(is,1,l) cq2p=rr(is,2,l)+rr(-is,2,l)+ui*ri(is,2,l) else cq0p=rr(is,0,l) cq0m=rr(is,0,l) cq1p=rr(is,1,l) cq1m=rr(is,1,l) cq2p=rr(is,2,l) end if c ca11=ca11+is**2*asl*cq0p ca12=ca12+is*l*asl*cq0m ca22=ca22+bsl*cq0p ca13=ca13+is*asl*cq1m/yg ca23=ca23+l*asl*cq1p/yg ca33=ca33+asl*cq2p/yg**2 end do epsl(1,1,l) = - xg*ca11*fal epsl(1,2,l) = + ui*xg*ca12*fal epsl(2,2,l) = - xg*ca22*fal epsl(1,3,l) = - xg*ca13*fal epsl(2,3,l) = - ui*xg*ca23*fal epsl(3,3,l) = - xg*ca33*fal end do c cq2p=rr(0,2,0) e330=1.0d0+xg*cq2p c epsl(1,1,1) = 1.d0 + epsl(1,1,1) epsl(2,2,1) = 1.d0 + epsl(2,2,1) c do l=1,lrm epsl(2,1,l) = - epsl(1,2,l) epsl(3,1,l) = epsl(1,3,l) epsl(3,2,l) = - epsl(2,3,l) end do c return end c c c subroutine hermitian(rr,lrm) implicit real*8(a-h,o-z) parameter(tmax=5,npts=500) dimension rr(-lrm:lrm,0:2,0:lrm) real*8 ttv(npts+1),extv(npts+1) c common/ttex/ttv,extv c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/warm/iwarm,ilarm common/cri/cr,ci c do n=-lrm,lrm do k=0,2 do m=0,lrm rr(n,k,m)=0.0d0 end do end do end do c llm=min(3,lrm) c dt=2.0d0*tmax/dble(npts) bth2=2.0d0/amu bth=sqrt(bth2) amu2=amu*amu amu4=amu2*amu2 amu6=amu4*amu2 c do i = 1, npts+1 t = ttv(i) rxt=sqrt(1.0d0+t*t/(2.0d0*amu)) x = t*rxt upl2=bth2*x**2 upl=bth*x gx=1.0d0+t*t/amu exdx=cr*extv(i)*gx/rxt*dt c n1=1 if(iwarm.gt.2) n1=-llm c do n=n1,llm nn=abs(n) gr=anpl*upl+n*yg zm=-amu*(gx-gr) s=amu*(gx+gr) zm2=zm*zm zm3=zm2*zm call calcei3(zm,fe0m) c do m=nn,llm if(n.eq.0.and.m.eq.0) then rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2 else ffe=0.0d0 if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))/amu2 if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s . +s*s*(1.0d0+zm-zm2*fe0m))/amu4 if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0* . (20.0d0-8.0d0*zm+zm2)+s**3*(2.0d0+zm+zm2-zm3*fe0m))/amu6 c rr(n,0,m) = rr(n,0,m) + exdx*ffe rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2 end if c end do end do end do c if(iwarm.gt.2) return c sy1=1.0d0+yg sy2=1.0d0+yg*2.0d0 sy3=1.0d0+yg*3.0d0 c bth4=bth2*bth2 bth6=bth4*bth2 c anpl2=anpl*anpl c rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2) . +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2)) rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2)) rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2)) rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2 . /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/ . sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3)) rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ . 1.5d0*anpl2/sy1**2)) rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0* . anpl2/sy1**2)) c if(llm.gt.1) then c rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+ . bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2)) rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2)) rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2)) rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* . (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4* . (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2 . -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2* . (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2)) rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* . (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2)) rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* . (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4* . (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2 . -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2* . (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2)) rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* . (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2)) c if(llm.gt.2) then c rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4* . (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2)) rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2)) rr(-1,0,3) = -12.0d0*bth4/sy1* . (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+ . bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2) . /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2)) rr(-1,2,3) = -6.0d0*bth6/sy1* . (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2)) c rr(-2,0,3) = -12.0d0*bth4/sy2* . (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+ . bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2) . /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2)) rr(-2,2,3) = -6.0d0*bth6/sy2* . (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2)) c rr(-3,0,3) = -12.0d0*bth4/sy3* . (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+ . bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2) . /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4)) rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2)) rr(-3,2,3) = -6.0d0*bth6/sy3* . (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2)) c end if end if c return end subroutine hermitian_2(rr,lrm) implicit real*8(a-h,o-z) parameter(tmax=5,npts=500) dimension rr(-lrm:lrm,0:2,0:lrm) parameter(epsa=0.0d0,epsr=1.0d-4) parameter (lw=5000,liw=lw/4) dimension w(lw),iw(liw) external fhermit c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/warm/iwarm,ilarm common/cri/cr,ci common/nmhermit/n,m,ih c do n=-lrm,lrm do k=0,2 do m=0,lrm rr(n,k,m)=0.0d0 end do end do end do c llm=min(3,lrm) c bth2=2.0d0/amu bth=sqrt(bth2) amu2=amu*amu amu4=amu2*amu2 amu6=amu4*amu2 c n1=1 if(iwarm.gt.10) n1=-llm c do n=n1,llm nn=abs(n) do m=nn,llm if(n.eq.0.and.m.eq.0) then ih=2 c call dqagi(fhermit,bound,2,epsa,epsr,resfh, call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh, . epp,neval,ier,liw,lw,last,iw,w) rr(n,2,m) = resfh else do ih=0,2 c call dqagi(fhermit,bound,2,epsa,epsr,resfh, call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh, . epp,neval,ier,liw,lw,last,iw,w) rr(n,ih,m) = resfh end do end if end do end do c write(83,'(12(1x,e12.5))') c . yg,rr(1,0,1),rr(1,1,1),rr(1,2,1) if(iwarm.gt.10) return c sy1=1.0d0+yg sy2=1.0d0+yg*2.0d0 sy3=1.0d0+yg*3.0d0 c bth4=bth2*bth2 bth6=bth4*bth2 c anpl2=anpl*anpl c rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2) . +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2)) rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2)) rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2)) rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2 . /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/ . sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3)) rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ . 1.5d0*anpl2/sy1**2)) rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0* . anpl2/sy1**2)) c if(llm.gt.1) then c rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+ . bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2)) rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2)) rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2)) rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* . (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4* . (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2 . -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2* . (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2)) rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* . (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2)) rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* . (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4* . (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2 . -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2* . (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2)) rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* . (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2)) c if(llm.gt.2) then c rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4* . (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2)) rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2)) rr(-1,0,3) = -12.0d0*bth4/sy1* . (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+ . bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2) . /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2)) rr(-1,2,3) = -6.0d0*bth6/sy1* . (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2)) c rr(-2,0,3) = -12.0d0*bth4/sy2* . (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+ . bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2) . /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2)) rr(-2,2,3) = -6.0d0*bth6/sy2* . (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2)) c rr(-3,0,3) = -12.0d0*bth4/sy3* . (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+ . bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2) . /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4)) rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2* . (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2)) rr(-3,2,3) = -6.0d0*bth6/sy3* . (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2)) c end if end if c return end function fhermit(t) implicit real*8 (a-h,o-z) c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/cri/cr,ci common/nmhermit/n,m,ih bth2=2.0d0/amu bth=sqrt(bth2) amu2=amu*amu amu4=amu2*amu2 amu6=amu4*amu2 rxt=sqrt(1.0d0+t*t/(2.0d0*amu)) x = t*rxt upl2=bth2*x**2 upl=bth*x gx=1.0d0+t*t/amu exdxdt=cr*exp(-t*t)*gx/rxt nn=abs(n) gr=anpl*upl+n*yg zm=-amu*(gx-gr) s=amu*(gx+gr) zm2=zm*zm zm3=zm2*zm call calcei3(zm,fe0m) ffe=0.0d0 uplh=upl**ih if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))*uplh/amu2 if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s . +s*s*(1.0d0+zm-zm2*fe0m))*uplh/amu4 if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*(20.0d0-8.0d0*zm+zm2) . +s**3*(2.0d0+zm+zm2-zm3*fe0m))*uplh/amu6 fhermit= exdxdt*ffe return end c c c subroutine antihermitian(ri,lrm) implicit none integer lmx,nmx,lrm,n,k,m,mm real*8 rpi parameter(rpi=1.77245385090552d0) parameter(lmx=20,nmx=lmx+2) real*8 fsbi(nmx) real*8 ri(lrm,0:2,lrm) real*8 yg,amu,anpl,cmu,anpl2,dnl,cr,ci,ygn,rdu2,rdu,anpr real*8 du,ub,aa,up,um,gp,gm,xp,xm,eep,eem,ee,cm,cim real*8 fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0 real*8 fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m real*8 fact c common/ygyg/yg common/amut/amu common/nplr/anpl,anpr common/uu/ub,cmu common/resah/anpl2,dnl c common/cri/cr,ci c do n=1,lrm do k=0,2 do m=1,lrm ri(n,k,m)=0.0d0 end do end do end do c dnl=1.0d0-anpl2 cmu=anpl*amu c do n=1,lrm ygn=n*yg rdu2=ygn**2-dnl if(rdu2.gt.0.0d0) then rdu=sqrt(rdu2) du=rdu/dnl ub=anpl*ygn/dnl aa=amu*anpl*du if (abs(aa).gt.5.0d0) then up=ub+du um=ub-du gp=anpl*up+ygn gm=anpl*um+ygn xp=up+1.0d0/cmu xm=um+1.0d0/cmu eem=exp(-amu*(gm-1.0d0)) eep=exp(-amu*(gp-1.0d0)) fi0p0=-1.0d0/cmu fi1p0=-xp/cmu fi2p0=-(1.0d0/cmu**2+xp*xp)/cmu fi0m0=-1.0d0/cmu fi1m0=-xm/cmu fi2m0=-(1.0d0/cmu**2+xm*xm)/cmu do m=1,lrm fi0p1=-2.0d0*m*(fi1p0-ub*fi0p0)/cmu fi0m1=-2.0d0*m*(fi1m0-ub*fi0m0)/cmu fi1p1=-((1.0d0+2.0d0*m)*fi2p0-2.0d0*(m+1.0d0)*ub*fi1p0 . +up*um*fi0p0)/cmu fi1m1=-((1.0d0+2.0d0*m)*fi2m0-2.0d0*(m+1.0d0)*ub*fi1m0 . +up*um*fi0m0)/cmu fi2p1=(2.0d0*(1.0d0+m)*fi1p1-2.0d0*m* . (ub*fi2p0-up*um*fi1p0))/cmu fi2m1=(2.0d0*(1.0d0+m)*fi1m1-2.0d0*m* . (ub*fi2m0-up*um*fi1m0))/cmu if(m.ge.n) then ri(n,0,m)=0.5d0*ci*dnl**m*(fi0p1*eep-fi0m1*eem) ri(n,1,m)=0.5d0*ci*dnl**m*(fi1p1*eep-fi1m1*eem) ri(n,2,m)=0.5d0*ci*dnl**m*(fi2p1*eep-fi2m1*eem) end if fi0p0=fi0p1 fi1p0=fi1p1 fi2p0=fi2p1 fi0m0=fi0m1 fi1m0=fi1m1 fi2m0=fi2m1 end do else ee=exp(-amu*(ygn-1.0d0+anpl*ub)) call ssbi(aa,n,lrm,fsbi) do m=n,lrm cm=rpi*fact(m)*du**(2*m+1) cim=0.5d0*ci*dnl**m mm=m-n+1 fi0m=cm*fsbi(mm) fi1m=-0.5d0*aa*cm*fsbi(mm+1) fi2m=0.5d0*cm*(fsbi(mm+1)+0.5d0*aa*aa*fsbi(mm+2)) ri(n,0,m)=cim*ee*fi0m ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m) ri(n,2,m)=cim*ee*(du*du*fi2m+2.0d0*du*ub*fi1m+ub*ub*fi0m) end do end if end if end do c return end c subroutine ssbi(zz,n,l,fsbi) implicit none integer n,l,nmx,lmx,k,m,mm real*8 eps,c0,c1,sbi,zz real*8 gamm parameter(eps=1.0d-10,lmx=20,nmx=lmx+2) real*8 fsbi(nmx) do m=n,l+2 c0=1.0d0/gamm(dble(m)+1.5d0) sbi=c0 do k=1,50 c1=c0*0.25d0*zz**2/(dble(m+k)+0.5d0)/dble(k) sbi=sbi+c1 if(c1/sbi.lt.eps) exit c0=c1 end do mm=m-n+1 fsbi(mm)=sbi end do return end c c Weakly relativistic dielectric tensor c computation of dielectric tensor elements c Krivenki and Orefice, JPP 30,125 (1983) c subroutine diel_tens_wr(ce330,cepsl,lrm) implicit real*8 (a-b,d-h,o-z) implicit complex*16 (c) dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2),cepsl(3,3,lrm) parameter(cui=(0.0d0,1.0d0)) c common/xgxg/xg common/ygyg/yg common/nplr/anpl,anprf common/amut/amu c anpl2=anpl*anpl c call fsup(cefp,cefm,lrm) c do l=1,lrm lm=l-1 fcl=0.5d0**l*((1.0d0/yg)**2/amu)**lm*fact(2*l)/fact(l) ca11=(0.d0,0.d0) ca12=(0.d0,0.d0) ca13=(0.d0,0.d0) ca22=(0.d0,0.d0) ca23=(0.d0,0.d0) ca33=(0.d0,0.d0) do is=0,l k=l-is w=(-1.0d0)**k c asl=w/(fact(is+l)*fact(l-is)) bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) c cq0p=amu*cefp(is,0) cq0m=amu*cefm(is,0) cq1p=amu*anpl*(cefp(is,0)-cefp(is,1)) cq1m=amu*anpl*(cefm(is,0)-cefm(is,1)) cq2p=cefp(is,1)+amu*anpl2*(cefp(is,2) . +cefp(is,0)-2.0d0*cefp(is,1)) c ca11=ca11+is**2*asl*cq0p ca12=ca12+is*l*asl*cq0m ca22=ca22+bsl*cq0p ca13=ca13+is*asl*cq1m/yg ca23=ca23+l*asl*cq1p/yg ca33=ca33+asl*cq2p/yg**2 end do cepsl(1,1,l) = - xg*ca11*fcl cepsl(1,2,l) = + cui*xg*ca12*fcl cepsl(2,2,l) = - xg*ca22*fcl cepsl(1,3,l) = - xg*ca13*fcl cepsl(2,3,l) = - cui*xg*ca23*fcl cepsl(3,3,l) = - xg*ca33*fcl end do c cq2p=cefp(0,1)+amu*anpl2*(cefp(0,2)+cefp(0,0)-2.0d0*cefp(0,1)) ce330=1.0d0-xg*amu*cq2p c cepsl(1,1,1) = 1.d0 + cepsl(1,1,1) cepsl(2,2,1) = 1.d0 + cepsl(2,2,1) c do l=1,lrm cepsl(2,1,l) = - cepsl(1,2,l) cepsl(3,1,l) = cepsl(1,3,l) cepsl(3,2,l) = - cepsl(2,3,l) end do c return end c c c subroutine fsup(cefp,cefm,lrm) implicit real*8 (a-b,d-h,o-z) implicit complex*16 (c) parameter(apsicr=0.7d0) parameter(cui=(0.0d0,1.0d0)) dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2) c common/ygyg/yg common/nplr/anpl,anprf common/amut/amu c psi=sqrt(0.5d0*amu)*anpl apsi=abs(psi) c do is=0,lrm alpha=anpl*anpl/2.0d0+is*yg-1.0d0 phi2=amu*alpha phim=sqrt(abs(phi2)) if (alpha.ge.0) then xp=psi-phim yp=0.0d0 xm=-psi-phim ym=0.0d0 x0=-phim y0=0.0d0 else xp=psi yp=phim xm=-psi ym=phim x0=0.0d0 y0=phim end if call zetac (xp,yp,zrp,zip,iflag) call zetac (xm,ym,zrm,zim,iflag) c czp=dcmplx(zrp,zip) czm=dcmplx(zrm,zim) cf12=(0.0d0,0.0d0) if (alpha.ge.0) then if (alpha.ne.0) cf12=-(czp+czm)/(2.0d0*phim) else cf12=-cui*(czp+czm)/(2.0d0*phim) end if c if(apsi.gt.apsicr) then cf32=-(czp-czm)/(2.0d0*psi) else cphi=phim if(alpha.lt.0) cphi=-cui*phim call zetac (x0,y0,zr0,zi0,iflag) cz0=dcmplx(zr0,zi0) cdz0=2.0d0*(1.0d0-cphi*cz0) cf32=cdz0 end if c cf0=cf12 cf1=cf32 cefp(is,0)=cf32 cefm(is,0)=cf32 do l=1,is+2 iq=l-1 if(apsi.gt.apsicr) then cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 else cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) end if ir=l-is if(ir.ge.0) then cefp(is,ir)=cf2 cefm(is,ir)=cf2 end if cf0=cf1 cf1=cf2 end do c if(is.ne.0) then c alpha=anpl*anpl/2.0d0-is*yg-1.0d0 phi2=amu*alpha phim=sqrt(abs(phi2)) if (alpha.ge.0.0d0) then xp=psi-phim yp=0.0d0 xm=-psi-phim ym=0.0d0 x0=-phim y0=0.0d0 else xp=psi yp=phim xm=-psi ym=phim x0=0.0d0 y0=phim end if call zetac (xp,yp,zrp,zip,iflag) call zetac (xm,ym,zrm,zim,iflag) c czp=dcmplx(zrp,zip) czm=dcmplx(zrm,zim) c cf12=(0.0d0,0.0d0) if (alpha.ge.0) then if (alpha.ne.0.0d0) cf12=-(czp+czm)/(2.0d0*phim) else cf12=-cui*(czp+czm)/(2.0d0*phim) end if if(apsi.gt.apsicr) then cf32=-(czp-czm)/(2.0d0*psi) else cphi=phim if(alpha.lt.0) cphi=-cui*phim call zetac (x0,y0,zr0,zi0,iflag) cz0=dcmplx(zr0,zi0) cdz0=2.0d0*(1.0d0-cphi*cz0) cf32=cdz0 end if c cf0=cf12 cf1=cf32 do l=1,is+2 iq=l-1 if(apsi.gt.apsicr) then cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 else cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) end if ir=l-is if(ir.ge.0) then cefp(is,ir)=cefp(is,ir)+cf2 cefm(is,ir)=cefm(is,ir)-cf2 end if cf0=cf1 cf1=cf2 end do c end if c end do c return end subroutine eccd(effjcd) use green_func_p implicit none real*8 anum,denom,resp,resj,effjcd,coullog,dens,tekev real*8 anucc,canucc,ddens real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff real*8 rbn,rbx,btot,bmin,bmax,alams,fp0s,pa,fc real*8 fjch,fjncl,fjch0,fconic real*8 cst,cst2 integer ieccd,nharm,nhf,nhn external fjch,fjncl,fjch0 parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) parameter(vcsi=2.99792458d+8) parameter(qe=4.8032068d-10,me=mesi*1.d3,vc=vcsi*1.d2) parameter(canucc=4.0d13*pi*qe**4/(me**2*vc**3)) parameter(ceff=qesi/(mesi*vcsi)) common/nharm/nharm,nhf common/nhn/nhn common/ieccd/ieccd common/tete/tekev common/dens/dens,ddens common/zz/Zeff common/btot/btot common/bmxmn/bmax,bmin common/fc/fc common/ncl/rbx common/cohen/rbn,alams,pa,fp0s common/cst/cst,cst2 anum=0.0d0 denom=0.0d0 effjcd=0.0d0 coullog=24.0d0-log(1.0d4*sqrt(0.1d0*dens)/tekev) anucc=canucc*dens*coullog c nhf=nharm select case(ieccd) case(1) c cohen model c rbn=B/B_min c rbx=B/B_max c cst2=1.0d0-B/B_max c alams=sqrt(1-B_min/B_max) c Zeff < 31 !!! c fp0s= P_a (alams) rbn=btot/bmin rbx=btot/bmax cst2=1.0d0-rbx if(cst2.lt.1d-6) cst2=0.0d0 cst=sqrt(cst2) alams=sqrt(1.0d0-bmin/bmax) pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0 fp0s=fconic(alams,pa,0) do nhn=nharm,nhf call curr_int(fjch,resj,resp) anum=anum+resj denom=denom+resp end do case(2:9) cst=0.0d0 cst2=0.0d0 do nhn=nharm,nhf call curr_int(fjch0,resj,resp) anum=anum+resj denom=denom+resp end do case(10:11) c neoclassical model: c ft=1-fc trapped particle fraction c rzfc=(1+Zeff)/fc rbn=btot/bmin rbx=btot/bmax cst2=1.0d0-rbx if(cst2.lt.1d-6) cst2=0.0d0 cst=sqrt(cst2) call SpitzFuncCoeff(Tekev,Zeff,fc) do nhn=nharm,nhf call curr_int(fjncl,resj,resp) anum=anum+resj denom=denom+resp end do nhn=nhn-1 CASE DEFAULT print*,'ieccd undefined' end select c c effjpl = / /(B_min/) [A m /W] c if(denom.gt.0.0d0) effjcd=-ceff*anum/(anucc*denom) return end c c c subroutine curr_int(fcur,resj,resp) implicit real*8(a-h,o-z) parameter(epsa=0.0d0,epsr=1.0d-2) parameter(xxcr=16.0d0) parameter (lw=5000,liw=lw/4) dimension w(lw),iw(liw) external fcur,fpp common/nhn/nhn common/ygyg/yg common/nplr/anpl,anpr common/amut/amu common/gg/uplp,uplm,ygn common/ierr/ierr common/iokh/iokhawa common/cst/cst,cst2 c EC power and current densities anpl2=anpl*anpl dnl=1.0d0-anpl2 ygn=nhn*yg ygn2=ygn*ygn resj=0.0d0 resj1=0.0d0 resj2=0.0d0 resp=0.0d0 c rdu2=anpl2+ygn2-1.0d0 uplp=0.0d0 uplm=0.0d0 upltp=0.0d0 upltm=0.0d0 c if (rdu2.ge.0.0d0) then rdu=sqrt(rdu2) uplp=(anpl*ygn+rdu)/dnl uplm=(anpl*ygn-rdu)/dnl c uu1=uplm uu2=uplp xx1=amu*(anpl*uu1+ygn-1.0d0) xx2=amu*(anpl*uu2+ygn-1.0d0) c if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl duu=abs(uu1-uu2) c if(duu.gt.1.d-6) then call dqags(fpp,uu1,uu2,epsa,epsr,resp,epp,neval,ier, . liw,lw,last,iw,w) if (ier.gt.0) ierr=90 end if c rdu2t=cst2*anpl2+ygn2-1.0d0 c if (rdu2t.lt.0.0d0.or.cst2.eq.0.0d0) then c c resonance curve does not cross the trapping region c iokhawa=0 if(duu.gt.1.d-4) then call dqags(fcur,uu1,uu2,epsa,epsr, . resj,ej,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) ierr=91 end if else c c resonance curve crosses the trapping region c iokhawa=1 rdut=sqrt(rdu2t) upltm=(cst2*anpl*ygn-cst*rdut)/(1.0d0-cst2*anpl2) upltp=(cst2*anpl*ygn+cst*rdut)/(1.0d0-cst2*anpl2) c uu1=uplm uu2=upltm xx1=amu*(anpl*uu1+ygn-1.0d0) xx2=amu*(anpl*uu2+ygn-1.0d0) if(xx1.lt.xxcr.or.xx2.lt.xxcr) then if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl duu=abs(uu1-uu2) if(duu.gt.1.d-6) then call dqags(fcur,uu1,uu2,epsa,epsr, . resj1,ej1,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) then if (abs(resj1).lt.1.0d-10) then resj1=0.0d0 else ierr=92 end if end if end if end if c uu1=upltp uu2=uplp xx1=amu*(anpl*uu1+ygn-1.0d0) xx2=amu*(anpl*uu2+ygn-1.0d0) if(xx1.lt.xxcr.or.xx2.lt.xxcr) then if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl duu=abs(uu1-uu2) if(duu.gt.1.d-6) then call dqags(fcur,uu1,uu2,epsa,epsr, . resj2,ej2,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) then if(ier.ne.2) ierr=93 end if end if end if resj=resj1+resj2 end if end if return end c c computation of integral for power density, integrand function fpp c c ith=0 : polarization term = const c ith=1 : polarization term Larmor radius expansion to lowest order c ith=2 : full polarization term (J Bessel) c function fpp(upl) implicit real*8 (a-h,o-z) complex*16 ex,ey,ez,ui,emxy,epxy parameter(ui=(0.0d0,1.0d0)) common/ygyg/yg common/nplr/anpl,anpr common/amut/amu common/gg/uplp,uplm,ygn common/epolar/ex,ey,ez common/nprw/anprre,anprim common/ithn/ithn common/nhn/nhn upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm) gam=anpl*upl+ygn ee=exp(-amu*(gam-1)) thn2=1.0d0 thn2u=upr2*thn2 if(ithn.gt.0) then emxy=ex-ui*ey epxy=ex+ui*ey if(upr2.gt.0.0d0) then upr=sqrt(upr2) bb=anprre*upr/yg if(ithn.eq.1) then c Larmor radius expansion polarization term at lowest order cth=1.0d0 if(nhn.gt.1) cth=(0.5d0*bb)**(nhn-1)*nhn/fact(nhn) thn2=(0.5d0*cth*abs(emxy+ez*anprre*upl/ygn))**2 thn2u=upr2*thn2 else c Full polarization term nm=nhn-1 np=nhn+1 ajbnm=dbesjn(nm, bb) ajbnp=dbesjn(np, bb) ajbn=dbesjn(nhn, bb) thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0d0))**2 end if end if end if fpp=ee*thn2u return end c computation of integral for current density c fjch integrand for Cohen model with trapping c fjch0 integrand for Cohen model without trapping c fjncl integrand for momentum conserv. model K(u) from Maruschenko c gg=F(u)/u with F(u) as in Cohen paper function fjch(upl) implicit real*8 (a-h,o-z) common/gg/uplp,uplm,ygn common/nplr/anpl,anpr common/zz/Zeff common/cohen/rb,alams,pa,fp0s upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm) gam=anpl*upl+ygn u2=gam*gam-1.0d0 u=sqrt(u2) z5=Zeff+5.0d0 xi=1.0d0/z5**2 xib=1.0d0-xi xibi=1.0d0/xib fu2b=1.0d0+xib*u2 fu2=1.0d0+xi*u2 gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2) gg=u*gu/z5 dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5 alam=sqrt(1.0d0-upr2/u2/rb) fp0=fconic(alam,pa,0) dfp0=-(pa*pa/2.0d0+0.125d0) if (alam.lt.1.0d0) then dfp0=-fconic(alam,pa,1)/sqrt(1.0d0-alam**2) end if fh=alam*(1.0d0-alams*fp0/(alam*fp0s)) dfhl=1.0d0-alams*dfp0/fp0s eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam) if(upl.lt.0.0d0) eta=-eta fjch=eta*fpp(upl) return end function fjch0(upl) implicit real*8 (a-h,o-z) common/gg/uplp,uplm,ygn common/nplr/anpl,anpr common/zz/Zeff gam=anpl*upl+ygn u2=gam*gam-1.0d0 u=sqrt(u2) z5=Zeff+5.0d0 xi=1.0d0/z5**2 xib=1.0d0-xi xibi=1.0d0/xib fu2b=1.0d0+xib*u2 fu2=1.0d0+xi*u2 gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2) gg=u*gu/z5 dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5 eta=anpl*gg+gam*upl*dgg/u fjch0=eta*fpp(upl) return end function fjncl(upl) use green_func_p implicit real*8 (a-h,o-z) common/gg/uplp,uplm,ygn common/nplr/anpl,anpr common/fc/fc common/ncl/hb common/psival/psinv common/amut/amu common/tete/tekev common/zz/Zeff gam=anpl*upl+ygn u2=gam*gam-1.0d0 u=sqrt(u2) upr2=u2-upl*upl bth=sqrt(2.0d0/amu) uth=u/bth call GenSpitzFunc(Tekev,Zeff,fc,uth,u,gam,fk,dfk) fk=fk*(4.0d0/amu**2) dfk=dfk*(2.0d0/amu)*bth alam=upr2/u2/hb psi=psinv call vlambda(alam,psi,fu,dfu) eta=gam*fu*dfk/u-2.0d0*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb if(upl.lt.0) eta=-eta fjncl=eta*fpp(upl) return end subroutine vlambda(alam,psi,fv,dfv) implicit real*8 (a-h,o-z) parameter (nnintp=101,nlam=41) parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ . njest*nnintp+nlest+54) parameter(kwrk=nnintp+nlam+njest+nlest+3) parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) external fpbisp dimension xxs(1),yys(1),ffs(1) dimension ch01(lw01),ch((njest-4)*(nlest-4)) dimension tjp(njest),tlm(nlest) dimension iwrk(kwrk),wrk(lwrk) common/coffvl/ch common/coffdvl/ch01 common/coffvlt/tjp,tlm common/coffvln/njpt,nlmt njp=njpt nlm=nlmt xxs(1)=sqrt(psi) yys(1)=alam call fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs, . wrk(1),wrk(5),iwrk(1),iwrk(2)) fv=ffs(1) iwp=1+(njp-4)*(nlm-5) iwl=iwp+4 call fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2, . xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2)) dfv=ffs(1) return end subroutine projxyzt(iproj,nfile) implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) c common/nray/nrayr,nrayth common/wrk/ywrk,ypwrk common/psinv11/psinv11 common/istep/istep common/ss/st c rtimn=1.d+30 rtimx=-1.d-30 c jd=1 if(iproj.eq.0) jd=nrayr-1 do j=1,nrayr,jd kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk c dx=ywrk(1,j,k)-ywrk(1,1,1) dy=ywrk(2,j,k)-ywrk(2,1,1) dz=ywrk(3,j,k)-ywrk(3,1,1) c dirx=ywrk(4,j,k) diry=ywrk(5,j,k) dirz=ywrk(6,j,k) dir=sqrt(dirx*dirx+diry*diry+dirz*dirz) c if(j.eq.1.and.k.eq.1) then csth1=dirz/dir snth1=sqrt(1.0d0-csth1**2) csps1=1.0d0 snps1=0.0d0 if(snth1.gt.0.0d0) then csps1=diry/(dir*snth1) snps1=dirx/(dir*snth1) end if end if xti= dx*csps1-dy*snps1 yti=(dx*snps1+dy*csps1)*csth1-dz*snth1 zti=(dx*snps1+dy*csps1)*snth1+dz*csth1 rti=sqrt(xti**2+yti**2) c dirxt= (dirx*csps1-diry*snps1)/dir diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir c if(k.eq.1) then xti1=xti yti1=yti zti1=zti rti1=rti end if if(.not.(iproj.eq.0.and.j.eq.1)) . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 c if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti c end do c c if(.not.(iproj.eq.0.and.j.eq.1)) c . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 if(iproj.eq.1) write(nfile,*) ' ' end do c write(nfile,*) ' ' c write(12,99) istep,st,psinv11,rtimn,rtimx return 99 format(i5,12(1x,e16.8e3)) 111 format(3i5,12(1x,e16.8e3)) end c c c subroutine pec(pabs,currt) implicit real*8(a-h,o-z) parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000) parameter(rtbc=1.0d0) parameter(pi=3.14159265358979d0) dimension psjki(jmx,kmx,nmx),iiv(jmx,kmx) dimension ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx) dimension pdjki(jmx,kmx,nmx),currj(jmx,kmx,nmx) dimension xxi(nmx),ypt(nmx),yamp(nmx) dimension rtab(nndmx),rhotv(nndmx) dimension rtab1(0:nndmx) dimension ajphiv(nndmx),dpdv(nndmx) dimension dvolt(nndmx),darea(nndmx) dimension ratjav(nndmx),ratjbv(nndmx),ratjplv(nndmx) dimension ajplv(nndmx),ajcdav(nndmx),ajcdbv(nndmx) dimension pins(nndmx),currins(nndmx),fi(nndmx) parameter(llmx=21) dimension isev(llmx) c common/nray/nrayr,nrayth common/istep/istep common/dsds/dst common/ipec/ipec,nnd common/ieccd/ieccd common/index_rt/index_rt c common/iiv/iiv common/psjki/psjki common/pcjki/ppabs,ccci common/dpjjki/pdjki,currj c common/cent/btrcen,rcen common/angles/alpha0,beta0 common/iieq/iequil common/parban/b0,rr0m,zr0m,rpam common/taumnx/taumn,taumx,pabstot,currtot common/jmxmn/rhot1,rhot2,aj1,aj2 c common/polcof/psipol,chipol c stf=istep*dst nd=nnd if(pabstot.gt.0.0d0) then do ll=1,llmx isev(ll)=0 end do intp=0 voli0=0.0d0 areai0=0.0d0 rtab1(0)=0.0d0 do it=1,nd drt=1.0d0/dble(nd-1) rt=dble(it-1)*drt if(it.lt.nd) then rt1=rt+drt/2.0d0 else rt1=rt end if rtab(it)=rt rtab1(it)=rt1 dpdv(it)=0.0d0 ajphiv(it)=0.0d0 if (ipec.eq.0) then psit=rt psit1=rt1 else psit=rt**2 psit1=rt1**2 end if if(iequil.eq.2) then rhotv(it)=frhotor(psit) call valpsispl(sqrt(psit1),voli1,dervoli,areai1,rrii, . rbavi,bmxi,bmni,fci,intp) dvolt(it)=abs(voli1-voli0) darea(it)=abs(areai1-areai0) voli0=voli1 areai0=areai1 call ratioj(sqrt(psit),ratjai,ratjbi,ratjpli) ratjav(it)=ratjai ratjbv(it)=ratjbi ratjplv(it)=ratjpli else rhotv(it)=sqrt(psit) area1=pi*psit1*rpam**2 voli1=2.0d0*pi*rr0m*area1 dvolt(it)=abs(voli1-voli0) darea(it)=abs(areai1-areai0) voli0=voli1 areai0=areai1 end if end do kkk=1 do j=1,nrayr if(j.gt.1) kkk=nrayth 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 idecr=-1 is=1 do i=1,ii if(ipec.eq.0) xxi(i)=abs(psjki(j,k,i)) if(ipec.eq.1) xxi(i)=sqrt(abs(psjki(j,k,i))) if(psjki(j,k,i).ge.0.and.psjki(j,k,i).le.rtbc) then ypt(i)=ppabs(j,k,i) yamp(i)=ccci(j,k,i) else ypt(i)=0.0d0 yamp(i)=0.0d0 end if if(ise0.eq.0) then if(xxi(i).lt.rtbc) then ise0=i isev(is)=i-1 is=is+1 end if else if (idecr.eq.-1) then if(xxi(i).gt.xxi(i-1)) then isev(is)=i-1 is=is+1 idecr=1 end if else if(xxi(i).gt.rtbc) exit if(xxi(i).lt.xxi(i-1)) then isev(is)=i is=is+1 idecr=-1 end if end if end if end do c isev(is)=i-1 ppa1=0.0d0 cci1=0.0d0 do iis=1,is-1 iis1=iis+1 iise0=isev(iis) iise=isev(iis1) if (mod(iis,2).ne.0) then idecr=-1 ind1=nd ind2=2 iind=-1 else idecr=1 ind1=1 ind2=nd iind=1 end if do ind=ind1,ind2,iind iind=ind if (idecr.eq.-1) iind=ind-1 rt1=rtab1(iind) call locatex(xxi,iise,iise0,iise,rt1,itb1) if(itb1.gt.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), . yamp(itb1+1),rt1,cci2) dppa=ppa2-ppa1 dpdv(ind)=dpdv(ind)+dppa didst=(cci2-cci1) ajphiv(ind)=ajphiv(ind)+didst ppa1=ppa2 cci1=cci2 end if end do end do end do end do h=1.0d0/dble(nd-1) rhotpav=0.0d0 drhotpav=0.0d0 rhotjav=0.0d0 rhotjava=0.0d0 rhot2java=0.0d0 fi=dpdv/h call simpson (nd,h,fi,spds) fi=rhotv*dpdv/h call simpson (nd,h,fi,rhotpav) fi=rhotv*rhotv*dpdv/h call simpson (nd,h,fi,rhot2pav) rhotpav=rhotpav/spds rhot2pav=rhot2pav/spds if (ieccd.ne.0) then fi=ajphiv/h call simpson (nd,h,fi,sccs) fi=rhotv*ajphiv/h call simpson (nd,h,fi,rhotjav) fi=abs(ajphiv)/h call simpson (nd,h,fi,sccsa) fi=rhotv*abs(ajphiv)/h call simpson (nd,h,fi,rhotjava) fi=rhotv*rhotv*abs(ajphiv)/h call simpson (nd,h,fi,rhot2java) rhotjav=rhotjav/sccs rhotjava=rhotjava/sccsa rhot2java=rhot2java/sccsa end if c factor sqrt(8)=2 sqrt(2) to match with full width c of gaussian profile drhot2pav=rhot2pav-rhotpav**2 drhotpav=sqrt(8.d0*drhot2pav) drhot2java=rhot2java-rhotjava**2 drhotjava=sqrt(8.d0*drhot2java) spds=0.0d0 sccs=0.0d0 do i=1,nd spds=spds+dpdv(i) sccs=sccs+ajphiv(i) pins(i)=spds currins(i)=sccs dpdv(i)=dpdv(i)/dvolt(i) ajphiv(i)=ajphiv(i)/darea(i) end do facpds=1.0d0 facjs=1.0d0 if(spds.gt.0.0d0) facpds=pabs/spds if(sccs.ne.0.0d0) facjs=currt/sccs do i=1,nd dpdv(i)=facpds*dpdv(i) ajphiv(i)=facjs*ajphiv(i) ajcdav(i)=ajphiv(i)*ratjav(i) ajcdbv(i)=ajphiv(i)*ratjbv(i) ajplv(i)=ajphiv(i)*ratjplv(i) end do call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx, . drhotp,drhopp) if(ieccd.ne.0) then call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi, . drhotjfi,drhopfi) xps=rhopfi call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx) else rhotjfi=0.0d0 rhopfi=0.0d0 ajmxfi=0.0d0 drhotjfi=0.0d0 drhopfi=0.0d0 xps=rhopp ratjamx=1.0d0 ratjbmx=1.0d0 ratjplmx=1.0d0 end if iif1=iiv(1,1) istmx=1 do i=2,iif1 if(psjki(1,1,i).ge.0.0d0) then if(pdjki(1,1,i).gt.pdjki(1,1,i-1)) istmx=i end if end do stmx=istmx*dst pins_02=0.0d0 pins_05=0.0d0 pins_085=0.0d0 xrhot=0.2d0 call locate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_02) xrhot=0.5d0 call locate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_05) xrhot=0.85d0 call locate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_085) else ajmxfi=0.0d0 dpdvmx=0.0d0 rhotjfi=1.0d0 rhotjav=1.0d0 rhotjava=1.0d0 rhotp=1.0d0 rhotpav=1.0d0 drhotjfi=0.0d0 drhotjava=0.0d0 drhotp=0.0d0 drhotpav=0.0d0 taumn=0 taumx=0 stmx=stf pins_02=0.0d0 pins_05=0.0d0 pins_085=0.0d0 c end of pabstot > 0 end if c dPdV [MW/m^3], Jcd [MA/m^2] if(ieccd.eq.0) currt=0.0d0 currtka=currt*1.0d3 write(6,*)' ' write(6,*)'#beta0 alpha0 Icd Pa Jphimx dPdVmx '// .'rhotj rhotjava rhotp rhotpav '// .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' write(6,99) beta0,alpha0,currtka,pabstot,ajmxfi,dpdvmx, . rhotjfi,rhotjava,rhotp,rhotpav, . drhotjava,drhotpav, . stmx,psipol,chipol,real(index_rt) write(7,99) currtka,pabstot,ajmxfi,dpdvmx, . rhotjfi,rhotjava,rhotp,rhotpav, . drhotjava,drhotpav,ratjbmx, . stmx,psipol,chipol,real(index_rt) do i=1,nd if (ipec.eq.0) then psin=rtab(i) rhop=sqrt(rtab(i)) else psin=rtab(i)**2 rhop=rtab(i) end if pinsr=0.0d0 if(pabstot.gt.0) pinsr=pins(i)/pabstot write(48,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i), . currins(i),pins(i),pinsr,real(index_rt) end do return 49 format(i5,20(1x,e12.5)) 99 format(30(1x,e12.5)) end c c c subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) implicit real*8(a-h,o-z) parameter(emn1=0.367879441171442d0) dimension xx(nd),yy(nd) common/jmxmn/rhotp,rhotm,ypkp,ypkm common/ipec/ipec,nnd common/iieq/iequil c call vmaxmini(yy,nd,ymn,ymx,imn,imx) ypk=0.0d0 xmx=xx(imx) xmn=xx(imn) if (abs(ymx).gt.abs(ymn)) then ipk=imx ypkp=ymx xpkp=xmx if(abs(ymn/ymx).lt.1d-2) ymn=0.0d0 ypkm=ymn xpkm=xmn else ipk=imn ypkp=ymn xpkp=xmn if(abs(ymx/ymn).lt.1d-2) ymx=0.0d0 ypkm=ymx xpkm=xmx end if if(xpkp.gt.0.0d0) then xpk=xpkp ypk=ypkp yye=ypk*emn1 call locatex(yy,nd,1,ipk,yye,ie1) if(ie1.gt.0.and.ie1.lt.nd) then call intlin(yy(ie1),xx(ie1),yy(ie1+1),xx(ie1+1),yye,rte1) else 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) else ipk=2 xpk=xx(2) ypk=yy(2) rte1=0.0d0 yye=ypk*emn1 call locate(yy,nd,yye,ie2) call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2) end if c ipk=1 if(ymx.ne.0.and.ymn.ne.0) ipk=2 c drhop=0.0d0 drhot=0.0d0 psie1=0.0d0 psie2=1.0d0 rhopmx=1.0d0 rhopmn=0.0d0 if (ie1.gt.0.or.ie2.gt.0) then if(ipec.eq.0) then rhopmx=sqrt(xpkp) rhopmn=sqrt(xpkm) psie2=rte2 psie1=rte1 drhop=sqrt(rte2)-sqrt(rte1) else rhopmx=xpkp rhopmn=xpkm drhop=rte2-rte1 psie2=rte2**2 psie1=rte1**2 end if end if c if(iequil.eq.2) then rhotmx=frhotor(rhopmx**2) rhotmn=frhotor(rhopmn**2) rhote2=frhotor(psie2) rhote1=frhotor(psie1) drhot=rhote2-rhote1 else rhotmx=rhopmx rhotmn=rhopmn drhot=drhop end if c if(ipk.eq.2) then drhop=-drhop drhot=-drhot end if ypk=ypkp rhotp=rhotmx rhotm=rhotmn c return end c subroutine polarcold(exf,eyif,ezf,elf,etf) implicit real*8(a-h,o-z) common/nplr/anpl,anpr common/xgxg/xg common/ygyg/yg common/mode/sox c c dcold dispersion c dielectric tensor (transposed) c c exf=0.0d0 c eyif=0.0d0 c ezf=0.0d0 c if(xg.le.0) return c anpl2=anpl*anpl anpr2=anpr*anpr an2=anpl2+anpr2 yg2=yg**2 dy2=1.0d0-yg2 aa=1.0d0-xg-yg2 e3=1.0d0-xg c if(xg.gt.0.0d0) then if (anpl.ne.0.0d0) then qq=xg*yg/(an2*dy2-aa) p=(anpr2-e3)/(anpl*anpr) ezf=1.0d0/sqrt(1.0d0+p*p*(1.0d0+qq*qq)) exf=p*ezf eyif=qq*exf else if(sox.lt.0.d0) then ezf=1 exf=0 eyif=0 else ezf=0 qq=-aa/(xg*yg) exf=1.0d0/sqrt(1.0d0+qq*qq) eyif=qq*exf end if end if elf=(anpl*ezf+anpr*exf)/sqrt(an2) etf=sqrt(1.0d0-elf*elf) else if(sox.lt.0.0d0) then ezf=1 exf=0 eyif=0 else ezf=0 exf=0.0d0 eyif=1.0d0 end if elf=0 etf=1 end if c return end subroutine pol_limit(qq,uu,vv) 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 complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) c common/anv/anv 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) anz=anv(3) anpl2=anpl*anpl dnl=1.0d0-anpl2 del0=sqrt(dnl**2+4.0d0*anpl2/yg**2) ffo=0.5d0*yg*(dnl+del0) ffx=0.5d0*yg*(dnl-del0) an2=anx*anx+any*any+anz*anz an=sqrt(an2) anxy=sqrt(anx*anx+any*any) sngam=(anz*anpl-an2*bv(3))/(an*anxy*anpr) csgam=-(any*bv(1)-anx*bv(2))/anxy/anpr csg2=csgam**2 sng2=sngam**2 ffo2=ffo*ffo ffx2=ffx*ffx deno=ffo2+anpl2 denx=ffx2+anpl2 xe2om=(ffo2*csg2+anpl2*sng2)/deno ye2om=(ffo2*sng2+anpl2*csg2)/deno xe2xm=(ffx2*csg2+anpl2*sng2)/denx ye2xm=(ffx2*sng2+anpl2*csg2)/denx 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) 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 real*8 pi real*8 psipol,chipol,psitr,chitr complex*16 ui,extr,eytr,eztr,ext,eyt complex*16 evin(3),evrfl(3) parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) common/xv/xv common/anv/anv common/polcof/psipol,chipol common/evt/ext,eyt common/wrefl/walln anv=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) c computation of reflection coordinates and normal to the wall 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 c rotation matrix from local to lab frame vv1(1)=anv(2) vv1(2)=-anv(1) vv1(3)=0.0d0 vv2(1)=anv(1)*anv(3) vv2(2)=anv(2)*anv(3) vv2(3)=-anv(1)*anv(1)-anv(2)*anv(2) vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2) vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2) vv3=anv evin=ext*vv1+eyt*vv2 c wave vector and electric field after reflection in lab frame anvrfl=anv-2.0d0* . (anv(1)*walln(1)+anv(2)*walln(2)+anv(3)*walln(3))*walln evrfl=-evin+2.0d0* . (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln vv1(1)=anvrfl(2) vv1(2)=-anvrfl(1) vv1(3)=0.0d0 vv2(1)=anvrfl(1)*anvrfl(3) vv2(2)=anvrfl(2)*anvrfl(3) vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2) vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2) vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2) vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2) extr=dot_product(vv1,evrfl) eytr=dot_product(vv2,evrfl) eztr=dot_product(vv3,evrfl) qqtr=abs(extr)**2-abs(eytr)**2 uutr=2.0d0*dble(extr*dconjg(eytr)) vvtr=2.0d0*dimag(extr*dconjg(eytr)) psitr=0.5d0*atan2(uutr,qqtr)*180.d0/pi chitr=0.5d0*asin(vvtr)*180.d0/pi ivac=2 anv=anvrfl xvrfl=xvout xv=xvout call vacuum_rt(xv,xvout,ivac) xv=xvout return 111 format(20(1x,e12.5)) end subroutine vacuum_rt(xvstart,xvend,ivac) use reflections, only : inters_linewall,inside implicit none integer*4 ivac integer nbb,nlim,i,imax parameter(nbb=1000) real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6) real*8 xv0(3) real*8 rlim(nbb),zlim(nbb) common/wrefl/walln common/limiter/rlim,zlim,nlim common/anv/anv common/dsds/dst common/psival/psinv common/densbnd/psdbnd common/dstvac/dstvac c ivac=1 first interface plasma-vacuum c ivac=2 second interface vacuum-plasma after wall reflection c ivac=-1 second interface vacuum-plasma WITHOUT wall reflection ! the real argument associated to xvstart is in a common block ! used by fwork!!! xv0=xvstart call inters_linewall(xv0/1.d2,anv,rlim(1:nlim),zlim(1:nlim), . nlim,smax,walln) smax=smax*1.d2 rrm=1.d-2*sqrt(xv0(1)**2+xv0(2)**2) zzm=1.d-2*xv0(3) if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then ! first wall interface is outside-inside if (dot_product(walln,walln)