subroutine gray_integration implicit none integer istep,nstep,istop,index_rt real*8 st,dst,strfl11 integer i real*8 st0 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(rhopin,nrho,dpdvout,ajcdout) implicit real*8 (a-h,o-z) parameter(zero=0.0d0) character*255 filenmeqq,filenmprf,filenmbm dimension rhopin(nrho),dpdvout(nrho),ajcdout(nrho) 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=609 call projxyzt(iproj,nfilp) end if c c print final results on screen c write(*,*) write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st if(iwarm.gt.0) then write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx else write(*,'(a,2f9.4)') 'taumn, taumx = ', zero,zero end if c write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot currtka =currtot*1.0d3 write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka c if (index_rt.eq.1) then if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq if(iequil.eq.1) write(*,*) 'ANALTYCAL EQUILIBRIUM' if(iprof.eq.1) write(*,*) 'PROFILE file : ',filenmprf if(iprof.eq.0) write(*,*) 'ANALTYCAL PROFILES' if(ibeam.ge.1) write(*,*) 'LAUNCHER CASE : ',filenmbm end if c c compute power and current density profiles for all rays c pabs=pabstot currt=currtot call pec(rhopin,nrho,pabs,currt,dpdvout,ajcdout) c return 99 format(20(1x,e16.8e3)) end c c c subroutine after_onestep(i,istop) use itm_constants, only : pi=>itm_pi implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36,nmx=8000) parameter(taucr=12.0d0,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 write(*,*) ' IERR = ', ierr if(ierr.eq.97) then c igrad=0 c ierr=0 else istop=1 c 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 write(*,*) ' ' c write(*,*) 'Coupled power fraction =',powa c write(*,*) ' ' c write(*,*) 'Input coupled power (MW) =',p0mw c write(*,*) ' ' 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(*,*) '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 write(*,*) 'istep = ',i if(nrayr.gt.1) then iproj=0 nfilp=608 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) use itm_constants, only : pi=>itm_pi implicit real*8 (a-h,o-z) parameter(ndim=6,jmx=31,kmx=36,nmx=8000) parameter(taucr=12.0d0) 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(604,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(617,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) then c if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then write(633,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(633,*) ' ' 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(604,*)' #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(608,*) ' #istep j k xt yt zt rt psin' write(609,*) ' #istep j k xt yt zt rt psin' write(617,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(633,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' write(612,*) ' #i sst psi w1 w2' write(607,*)'#Icd Pa Jphip dPdVp '// .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp' write(648,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' else c If(index_rt.eq.3) then write(604,*) ' ' write(608,*) ' ' write(609,*) ' ' write(617,*) ' ' write(612,*) ' ' write(648,*) ' ' c end if end if return end c c c subroutine read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, . dne, zeff, qsf, powin) use green_func_p, only:Setup_SpitzFunc use itm_constants, only : pi=>itm_pi implicit real*8 (a-h,o-z) integer ijetto, mr, mz, nrho, nbnd real*8 r(mr), z(mz), psin(mr,mz) real*8 psiax, psibnd, rax, zax, powin real*8 rbnd(nbnd), zbnd(nbnd) real*8 psijet(nrho), f(nrho), te(nrho), dne(nrho), . zeff(nrho), qsf(nrho) real*8 me character*8 wdat character*10 wtim character*255 filenmeqq,filenmprf,filenmbm parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) parameter(cvdr=pi/180.0d0) parameter(nmx=8000,nbb=5000) 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/Zeffan 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(602,file='gray.data',status= 'unknown') c c alpha0, beta0 (cartesian) launching angles c fghz wave frequency (GHz) c p0mw injected power (MW) c read(602,*) alpha0,beta0 read(602,*) fghz read(602,*) p0mw c power value overwritten with value passed as argument (in W) p0mw=powin*1.d-6 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(602,*) nrayr,nrayth,rwmax if(nrayr.eq.1) nrayth=1 c c x00,y00,z00 coordinates of launching point c read(602,*) 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(602,*) 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(602,*) ibeam read(602,*) 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(602,*) iequil,ixp c c iprof=0 :analytical density and temp. profiles c iprof>0 :numerical density and temp. profiles c read(602,*) 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(602,*) 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(602,*) ieccd if (ieccd.ge.10) call Setup_SpitzFunc c c ipec=0/1/-1 :pec profiles grid equispaced in psi/rhop/as input c nnd :number of grid steps for pec profiles +1 c read(602,*) ipec,nnd c input radial grid forced on output. nnd is ignored in this case ipec=-1 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(602,*) dst,nstep,istpr0,istpl0,idst read(602,*) igrad,ipass,rwallm read(602,*) iox,psipol0,chipol0 c c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004 c sspl spline parameter for psi interpolation c read(602,*) filenmeqq read(602,*) ipsinorm,sspl c psin grid is normalized ipsinorm=1 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(602,*) factb, iscal,factt,factn if (factb.eq.0.0d0) factb=1.0d0 c c psbnd value of psi ( > 1 ) of density boundary c read(602,*) filenmprf read(602,*) 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(602,*) sgnbphi,sgniphi,icocos c data are already in the coordinate system used by gray icocos=3 c c read parameters for analytical density and temperature c profiles if iprof=0 c if(iprof.eq.0) then read(602,*) dens0,aln1,aln2 read(602,*) te0,dte0,alt1,alt2 else read(602,*) dummy,dummy,dummy read(602,*) dummy,dummy,dummy,dummy end if read(602,*) zeffan c c read parameters for analytical simple cylindical equilibrium c if iequil=0,1 if(iequil.lt.2) then read(602,*) rr0,zr0,rpa read(602,*) b0 read(602,*) q0,qa,alq rr0m=rr0/1.0d2 zr0m=zr0/1.0d2 rpam=rpa/1.0d2 else read(602,*) dummy,dummy,dummy read(602,*) dummy read(602,*) dummy,dummy,dummy end if c close(unit=602) c if(nrayr.eq.1) igrad=0 if (nrayr.lt.5) then igrad=0 write(*,*) ' nrayr < 5 ! => OPTICAL CASE ONLY' write(*,*) 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 write(*,'(a,2f8.3)') 'alpha0, beta0 = ',alpha0,beta0 c r00=sqrt(x00**2+y00**2) phi0=acos(x00/r00) if(y00.lt.0) phi0=-phi0 write(*,'(a,4f8.3)') 'x00, y00, R00, z00 = ',x00,y00,r00,z00 write(*,*) 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 c nprof=98 c lprf=length(filenmprf) c open(file=filenmprf(1:lprf)//'.prf', c . status= 'unknown',unit=nprof) c read profiles from input arguments filenmprf='JETTO' call profdata(nrho, psijet, te, dne, zeff) c close(nprof) end if c c read equilibrium data from file if iequil=2 c if (iequil.eq.2) then c neqdsk=99 c leqq=length(filenmeqq) c open(file=filenmeqq(1:leqq)//'.eqdsk', c . status= 'unknown', unit=neqdsk) filenmeqq='JETTO' call equidata(ijetto, mr, mz, r, z, psin, psiax, psibnd, . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, qsf) c 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=638 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) trim(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) trim(filenmprf) else write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeffan 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) trim(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 : ',a) 906 format('# fact_B fact_T fact_n iscal : ',(3(1x,es12.5),i5)) 907 format('# PROFILES file : ',a) 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 use itm_constants, only : pi=>itm_pi implicit real*8(a-h,o-z) common/parban/b0,rr0m,zr0m,rpam common/parbres/bres c c print circular magnetic surfaces iequil=1 c write(631,*) '#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(631,111) i,rni,rrm,zzm end do write(631,*) ' ' write(631,*) ' ' end do c c print resonant B iequil=1 c write(630,*)'#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(630,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*255 filenmeqq,filenmprf,filenmbm parameter(nstrmx=50) c dimension alphastv(nstrmx),betastv(nstrmx),cbeta(nstrmx,4) dimension x00v(nstrmx),y00v(nstrmx),z00v(nstrmx) dimension cx0(nstrmx,4),cy0(nstrmx,4),cz0(nstrmx,4) dimension waist1v(nstrmx),waist2v(nstrmx) dimension rci1v(nstrmx),rci2v(nstrmx) dimension cwaist1(nstrmx,4),cwaist2(nstrmx,4) dimension crci1(nstrmx,4),crci2(nstrmx,4) dimension phi1v(nstrmx),phi2v(nstrmx) dimension cphi1(nstrmx,4),cphi2(nstrmx,4) c 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=603 lbm=length(filenmbm) open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam) c read(nfbeam,*) nisteer do i=1,nisteer if(ibeam.eq.1) then read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm, . waist1,dvir1,waist2,dvir2,delta 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 difcsg(alphastv,betastv,nisteer,iopt,cbeta,ier) call difcsg(alphastv,waist1v,nisteer,iopt,cwaist1,ier) call difcsg(alphastv,rci1v,nisteer,iopt,crci1,ier) call difcsg(alphastv,waist2v,nisteer,iopt,cwaist2,ier) call difcsg(alphastv,rci2v,nisteer,iopt,crci2,ier) call difcsg(alphastv,phi1v,nisteer,iopt,cphi1,ier) call difcsg(alphastv,phi2v,nisteer,iopt,cphi2,ier) call difcsg(alphastv,x00v,nisteer,iopt,cx0,ier) call difcsg(alphastv,y00v,nisteer,iopt,cy0,ier) call difcsg(alphastv,z00v,nisteer,iopt,cz0,ier) c if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then call vlocate(alphastv,nisteer,alpha0,k) dal=alpha0-alphastv(k) betst=fspli(cbeta,nisteer,k,dal) x00=fspli(cx0,nisteer,k,dal) y00=fspli(cy0,nisteer,k,dal) z00=fspli(cz0,nisteer,k,dal) wcsi=fspli(cwaist1,nisteer,k,dal) weta=fspli(cwaist2,nisteer,k,dal) rcicsi=fspli(crci1,nisteer,k,dal) rcieta=fspli(crci2,nisteer,k,dal) phiw=fspli(cphi1,nisteer,k,dal) phir=fspli(cphi2,nisteer,k,dal) else write(*,*) ' alpha0 outside table range !!!' if(alpha0.ge.alphastv(nisteer)) ii=nisteer if(alpha0.le.alphastv(1)) ii=1 betst=betastv(ii) x00=x00v(ii) y00=y00v(ii) z00=z00v(ii) wcsi=waist1v(ii) weta=waist2v(ii) rcicsi=rci1v(ii) rcieta=rci2v(ii) phiw=phi1v(ii) phir=phi2v(ii) end if beta0=betst c close(nfbeam) c return end c c c subroutine equidata(ijetto,mr,mz,r,z,psin2d,psiaxis,psiedge, . rax,zax,nbnd,rbnd,zbnd,mrho,psijet,fpjet,qjet) use itm_constants, only : pi=>itm_pi use reflections, only : inside implicit real*8 (a-h,o-z) integer ijetto,mr,mz,nbnd,mrho real*8 r(mr),z(mz),psin2d(mr,mz) real*8 psiaxis,psiedge,rax,zax real*8 rbnd(nbnd),zbnd(nbnd) real*8 psijet(mrho),fpjet(mrho),qjet(mrho) parameter(nnw=501,nnh=501) parameter(nbb=5000) dimension fpol(nnw),qpsi(nnw) dimension psi(nnw,nnh),rv(nnw),zv(nnh),psin(nnw,nnh),psinr(nnw) dimension rlim(nbb),zlim(nbb) 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) dimension rv1d(nnw*nnh),zv1d(nnw*nnh),wpsi(nnw*nnh) 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,nrho,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/phitedge/phitedge 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 nr=mr nz=mz nrho=mrho do i=1,nr rv(i)=r(i) enddo do j=1,nz zv(j)=z(j) do i=1,nr psin(i,j)=psin2d(i,j) enddo enddo do i=1,nrho psinr(i)=psijet(i) fpol(i)=fpjet(i) qpsi(i)=qjet(i) end do rmaxis=rax zmaxis=zax c assume last point of fpol array corresponds to vacuum value c and compute vacuum field at magnetic axis as reference rcen=rax btrcen=fpol(nrho)/rcen sgnbphi=sign(1.0d0,btrcen) c c length in m !!! c dr=rv(2)-rv(1) dz=zv(2)-zv(1) rmnm=rv(1) rmxm=rv(nr) zmnm=zv(1) zmxm=zv(nz) c psi function psia0=psiedge-psiaxis psia=psia0*factb sgniphi=sign(1.0d0,-psia0) c c do j=1,nz c do i=1,nr c write(620,2021) rv(i),zv(j),psin(i,j) c enddo c write(620,*) ' ' c enddo c c spline fitting/interpolation of psin(i,j) and derivatives c if (ijetto.eq.1) then c valid data on the whole grid do j=1,nz do i=1,nr c psi(i,j)=psin(i,j)*psia ij=nz*(i-1)+j fvpsi(ij)=psin(i,j) enddo enddo 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 else c valid data only inside boundary ij=0 do j=1,nz do i=1,nr if (.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle ij=ij+1 rv1d(ij)=rv(i) zv1d(ij)=zv(j) fvpsi(ij)=psin(i,j) wpsi(ij)=1.0d0 end do end do c add boundary to the fit if (ij+nbnd.le.nnw*nnh) then do i=1,nbnd ij=ij+1 rv1d(ij)=rbnd(i) zv1d(ij)=zbnd(i) fvpsi(ij)=1.0d0 wpsi(ij)=1.0d0 end do end if nrz=ij c c fit as a scattered set of points c nsr=nr+4 nsz=nz+4 call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier) !!!!!!!!!!!! Sistemare anche dipendenze Makefile scatterspl:... c if ier=-1 data are fitted using sspl=0 if(ier.eq.-1) then sspl=0.0d0 nsr=nr+4 nsz=nz+4 call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl, . rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier) end if nsrt=nsr nszt=nsz end if c c re-evaluate psi on the grid using the spline (only for debug) c c call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi, c . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier) c c do j=1,nz c do i=1,nr c ij=nz*(i-1)+j c psin(i,j)=ffvpsi(ij) cc psi(i,j)=psin(i,j)*psia c write(619,2021) rv(i),zv(j),psin(i,j) c enddo c write(619,*) ' ' c enddo c c2021 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 do i=1,nrho fpol(i)=fpol(i)*factb wf(i)=1.0d0 end do wf(nrho)=1.0d2 c c spline interpolation of fpol(i) c iopt=0 xb=0.0d0 xe=1.0d0 ssfp=0.01d0 call curfit(iopt,nrho,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,nrho,ier) fpolas=fpoli(nrho) c c no limiter shape provided yet nlim=0 c c compute max and min z of last closed surface c rbmin=rmaxis rbmax=rmaxis if (nbnd.gt.1) then zbmin=1.0d+30 zbmax=-1.0d+30 do i=1,nbnd if(zbnd(i).le.zbmin) then zbmin=zbnd(i) rbmin=rbnd(i) end if if(zbnd(i).ge.zbmax) then zbmax=zbnd(i) rbmax=rbnd(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 write(*,'(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 write(*,'(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 write(*,'(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 write(*,'(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 write(*,'(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 write(*,'(a,4f8.4)') 'no X-point ',rbmin,zbmin,rbmax,zbmax end if write(*,*) ' ' c compute B_toroidal on axis btaxis=fpol(1)/rmaxis btrcen=abs(btrcen)*factb write(*,'(a,f8.4)') 'factb = ',factb write(*,'(a,f8.4)') '|BT_centr|= ',btrcen write(*,'(a,f8.4)') '|BT_axis| = ',abs(btaxis) c compute normalized rho_tor from eqdsk q profile call rhotor(nrho) phitedge=abs(psia)*rhotsx*2*pi rrtor=sqrt(phitedge/abs(btrcen)/pi) call rhopol c write(*,*) rhotsx,phitedge,rrtor,abs(psia) 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(632,*) 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,nrho,qmin,qmax,iqmn,iqmx) if (q15.gt.qmin.and.q15.lt.qmax) then call surfq(q15,psi15) rhot15=frhotor(psi15) write(*,'(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) write(*,'(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 write(*,'(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,nrho,npp,nintp common/dens/dens,ddens c write(645,*) ' #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(645,111) psin,rhot,dens,te,qq,ajphi*1.d-6 c nst=nrho do i=2,nst psin=dble(i-1)/dble(nst-1) rhop=sqrt(psin) c call density(psin) te=temp(psin) c call vlocate(psinr,nrho,psin,ips) ips=min(max(ips,1),nrho-1) call intlin(psinr(ips),qpsi(ips),psinr(ips+1),qpsi(ips+1), . psin,qq) rhot=frhotor(psin) call tor_curr_psi(psin,ajphi) write(645,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,nrho,npp,nintp c c locate psi surface for q=qval c call vlocate(qpsi,nrho,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,nrho,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(640,113) j,rrj,zzk,btotal(j,k) enddo c write(640,*) ' ' enddo c c compute Btot=Bres and Btot=Bres/2 c write(630,*)'#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(630,113) inc,bbb,rrcb(inc),zzcb(inc) end do end if write(630,*) ' ' end do c return 113 format(i6,12(1x,e12.5)) end c c subroutine profdata(mrho, psijet, te, dne, zeff) implicit real*8 (a-h,o-z) c integer mrho real*8 psijet(mrho), te(mrho), dne(mrho), zeff(mrho) c 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,nrho,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 do i=1,mrho psrad(i)=psijet(i) terad(i)=te(i)*1.d-3*ffact**aat*factt derad(i)=dne(i)*1.d-19*ffact**aan*factn zfc(i)=zeff(i) wf(i)=1.0d0 end do psrad(1)=0.0d0 npp=mrho 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 gsplder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier) dnpp=densi(npp) ddnpp=ddensi(npp) c nu=2 call gsplder(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 common/cq/cq 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 difcsg(psinr,rhotnr,nr,iopt,crhot,ier) return end function fq_eq(psi) implicit real*8(a-h,o-z) parameter(nnw=501) dimension psinr(nnw),cq(nnw,4) common/psinr/psinr common/eqnn/nr,nz,nrho,npp,nintp common/cq/cq irt=int((nrho-1)*psi+1) if(irt.eq.0) irt=1 if(irt.eq.nrho) irt=nrho-1 dps=psi-psinr(irt) fq_eq=fspli(cq,nrho,irt,dps) 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,nrho,npp,nintp common/crhot/crhot c irt=int((nrho-1)*psi+1) if(irt.eq.0) irt=1 if(irt.eq.nrho) irt=nrho-1 dps=psi-psinr(irt) frhotor_eq=fspli(crhot,nrho,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,nrho,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=fspli(crhotq,nintp,ip,dps) return end subroutine rhopol implicit real*8(a-h,o-z) parameter(nnr=101,nrest=nnr+4) parameter(lwrkp=nnr*4+nrest*16) dimension rhop(nnr),rhot(nnr),rhopi(nnr) dimension trp(nrest),crp(nrest),wp(nrest) dimension wrkp(lwrkp),iwrkp(nrest) common/coffrtp/trp common/coffrn/nsrp common/coffrp/crp dr=1.0d0/dble(nnr-1) do i=1,nnr rhop(i)=(i-1)*dr psin=rhop(i)*rhop(i) rhot(i)=frhotor(psin) wp(i)=1.0d0 end do wp(1)=1.0d3 wp(nnr)=1.0d3 c spline interpolation of rhopol versus rhotor iopt=0 xb=0.0d0 xe=1.0d0 ss=0.00001 kspl=3 call curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,nsrp, . trp,crp,rp,wrkp,lwrkp,iwrkp,ier) c write(*,*) ier call splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier) do i=1,nnr write(644,*) rhop(i),rhot(i),rhopi(i) end do return end function frhopol(rhot) implicit real*8(a-h,o-z) parameter(nnr=101,nrest=nnr+4) dimension trp(nrest),crp(nrest),rrs(1),ffspl(1) common/coffrtp/trp common/coffrn/nsrp common/coffrp/crp rrs(1)=rhot call splev(trp,nsrp,crp,3,rrs,ffspl,1,ier) frhopol=ffspl(1) 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,nrho,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) use itm_constants, only : pi=>itm_pi 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 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 sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) write(*,*) ' sprofil =',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(631,111) ii,h,rcn(ii),zcn(ii) end do write(631,*) write(631,*) end if return 111 format(i6,12(1x,e12.5)) 99 format(12(1x,e12.5)) end c c c subroutine flux_average use itm_constants, only : pi=>itm_pi, itm_mu0 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(ccj=1.0d0/itm_mu0) 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 dadrhotv(nnintp),cdadrhot(nnintp,4) dimension dvdrhotv(nnintp),cdvdrhot(nnintp,4) 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,nrho,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 common/phitedge/phitedge common/cdadrhot/cdadrhot common/cdvdrhot/cdvdrhot c common/coffvl/ch common/coffdvl/ch01 common/coffvlt/tjp,tlm common/coffvln/njpt,nlmt c computation of flux surface averaged quantities write(631,*)' #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 dadrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) . *dadpsi/pi dvdrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2) . *dvdpsi/pi c c write(647,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(646,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// .' Area Vol |I_pl| qq fc ratioJa ratioJb' 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(646,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) . ,vratja(jp),vratjb(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 difcsg(rpstab,vvol,nintp,iopt,cvol,ier) iopt=0 call difcsg(rpstab,rbav,nintp,iopt,crbav,ier) iopt=0 call difcsg(rpstab,rri,nintp,iopt,crri,ier) iopt=0 call difcsg(rpstab,bmxpsi,nintp,iopt,cbmx,ier) iopt=0 call difcsg(rpstab,bmnpsi,nintp,iopt,cbmn,ier) iopt=0 call difcsg(rpstab,vratja,nintp,iopt,cratja,ier) iopt=0 call difcsg(rpstab,vratjb,nintp,iopt,cratjb,ier) iopt=0 call difcsg(rpstab,vratjpl,nintp,iopt,cratjpl,ier) iopt=0 call difcsg(rpstab,varea,nintp,iopt,carea,ier) iopt=0 call difcsg(rpstab,ffc,nintp,iopt,cfc,ier) iopt=0 call difcsg(rpstab,rhotqv,nintp,iopt,crhotq,ier) iopt=0 call difcsg(rpstab,dadrhotv,nintp,iopt,cdadrhot,ier) iopt=0 call difcsg(rpstab,dvdrhotv,nintp,iopt,cdvdrhot,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 function fdadrhot(rpsi) implicit real*8(a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp),cdadrhot(nnintp,4) common/pstab/rpstab common/eqnn/nr,nz,nrho,npp,nintp common/cdadrhot/cdadrhot ip=int((nintp-1)*rpsi+1) if(ip.eq.0) ip=1 if(ip.eq.nintp) ip=nintp-1 dps=rpsi-rpstab(ip) fdadrhot=fspli(cdadrhot,nintp,ip,dps) return end function fdvdrhot(rpsi) implicit real*8(a-h,o-z) parameter(nnintp=101) dimension rpstab(nnintp),cdvdrhot(nnintp,4) common/pstab/rpstab common/eqnn/nr,nz,nrho,npp,nintp common/cdvdrhot/cdvdrhot ip=int((nintp-1)*rpsi+1) if(ip.eq.0) ip=1 if(ip.eq.nintp) ip=nintp-1 dps=rpsi-rpstab(ip) fdvdrhot=fspli(cdvdrhot,nintp,ip,dps) return end 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) use itm_constants, only : pi=>itm_pi implicit real*8 (a-h,o-z) 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(lw01),cc10(lw10),cc02(lw02),cc20(lw20),cc11(lw11) dimension tfp(nrest),cfp(nrest),wrkfd(nrest) c common/eqnn/nr,nz,nrho,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) . write(*,'(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 gsplder(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) use itm_constants, only : pi=>itm_pi, itm_mu0 implicit real*8 (a-h,o-z) parameter(ccj=1.0d0/itm_mu0) 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 sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) write(*,*) ' sprofil =',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 gsplder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier) ddens=ffs(1) if(ier.gt.0) write(*,*) ier if(abs(dens).lt.1.0d-10) dens=0.0d0 end if if(dens.lt.0.0d0) write(*,*) ' 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,nrho,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 vlocate(psrad,npp,arg,k) if(k.eq.0) k=1 if(k.eq.npp) k=npp-1 dps=arg-psrad(k) temp=fspli(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,nrho,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 vlocate(psrad,npp,ps,k) if(k.eq.0) k=1 if(k.eq.npp) k=npp-1 dps=ps-psrad(k) fzeff=fspli(cz,npmx,k,dps) endif return end c c beam tracing initial conditions igrad=1 c subroutine ic_gb use itm_constants, only : pi=>itm_pi 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(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(633,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(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, . zero,zero,zero,zero, . zero,zero,zero,zero,one,zero,zero, . zero,zero,one,zero,zero,zero,zero,one ddr110=dd end if if(j.eq.nrayr.and.k.eq.1) then write(617,99) zero,ddr110,dd,ddi end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=608 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 use itm_constants, only : pi=>itm_pi 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(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(633,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(617,99) zero,zero,zero,zero write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, . zero,zero,zero,zero, . zero,zero,zero,zero,one,zero,zero, . zero,zero,one,zero,zero,zero,zero,one end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=608 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 use itm_constants, only : pi=>itm_pi 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(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(633,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(617,99) zero,zero,zero,zero write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, . zero,zero,zero,zero, . zero,zero,zero,zero,one,zero,zero, . zero,zero,one,zero,zero,zero,zero,one end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=608 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,nrho,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=fspli(carea,nintp,ip,dps) voli=fspli(cvol,nintp,ip,dps) dervoli=fsplid(cvol,nintp,ip,dps) rrii=fspli(crri,nintp,ip,dps) c if(intp.eq.0) return c rbavi=fspli(crbav,nintp,ip,dps) bmxi=fspli(cbmx,nintp,ip,dps) bmni=fspli(cbmn,nintp,ip,dps) fci=fspli(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,nrho,npp,nintp common/cratj/cratja,cratjb,cratjpl ip=int((nintp-1)*rpsi+1) if(ip.eq.0) ip=1 if(ip.eq.nintp) ip=nintp-1 dps=rpsi-rpstab(ip) ratjai=fspli(cratja,nintp,ip,dps) ratjbi=fspli(cratjb,nintp,ip,dps) ratjpli=fspli(cratjpl,nintp,ip,dps) return end c c c subroutine pabs_curr(i,j,k) use itm_constants, only : pi=>itm_pi 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 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 write(*,*) ' 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 write(*,*) ' 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) write(*,*) ' 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 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 write(*,*) '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(612,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(rhopin,nrho,pabs,currt,dpdvout,ajcdout) use itm_constants, only : pi=>itm_pi implicit real*8(a-h,o-z) parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000) parameter(rtbc=1.0d0) dimension rhopin(nrho),dpdvout(nrho),ajcdout(nrho) 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 if (ipec.ge.0) then nd=nnd else ! ipec=-1 uses rho_pol input grid nd=nrho end if 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 if (ipec.ge.0) then 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 else c radial coordinate of i-(i+1) interval left point rt=rhopin(it) if(it.lt.nd) then c radial coordinate of i-(i+1) interval mid point rt1=(rt+rhopin(it+1))/2.0d0 else rt1=rt end if 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 !!!!!!!!!! it should be isev(is)=i-1 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 !!!!!!!!!! is it safe? iind=ind !!!!!!!!!! iind reused in the loop! 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) !!!!!!!!!! equi-spaced grid expected in rhotpav=0.0d0 !!!!!!!!!! simpson subroutine drhotpav=0.0d0 !!!!!!!!!! wrong integrals until fixed rhotjav=0.0d0 !!!!!!!!!! (d)rhotpav,rhotjav,(d)rhotjava affected 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 rhpp=frhopol(rhotpav) rhpj=frhopol(rhotjava) dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp)) ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx, . drhotp,drhopp) if(ieccd.ne.0) then call profwidth(nd,rtab,ajphiv,rhotjfi,rhopfi,ajmxfi, . drhotjfi,drhopfi) xps=rhopfi c call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx) else rhotjfi=0.0d0 rhopfi=0.0d0 ajmxfi=0.0d0 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 vlocate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_02) xrhot=0.5d0 call vlocate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_05) xrhot=0.85d0 call vlocate(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(*,*)' ' write(*,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '// .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp' write(*,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp, . rhotjfi,rhotjava,rhotp,rhotpav, . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp write(607,99) currtka,pabstot,ajphip,dpdvp, . rhotjfi,rhotjava,rhotp,rhotpav, . drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol, . real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp 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(648,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i), . currins(i),pins(i),pinsr,real(index_rt) end do c fill output arguments only if ipec=-1 if (ipec.eq.-1) then do i=1,nd dpdvout(i)=dpdv(i) ajcdout(i)=ajcdbv(i) end do end if 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 vlocate(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) use itm_constants, only : pi=>itm_pi 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 beta0,alpha0,gam real*8 sox,psipol,chipol complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt parameter(ui=(0.0d0,1.0d0)) 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) use itm_constants, only : pi=>itm_pi 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 psipol,chipol,psitr,chitr complex*16 ui,extr,eytr,eztr,ext,eyt complex*16 evin(3),evrfl(3) parameter(ui=(0.0d0,1.0d0)) 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=5000) 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)