use const_and_precisions, only : wp_ use graydata_flags, only : ipass,igrad implicit none c local variables real(wp_) :: p0mw1 c common/external functions/variables integer :: ierr,index_rt real(wp_) :: sox,p0mw,powrfl,taumn,taumx,pabstot,currtot, . pabstott,currtott c common/ierr/ierr common/mode/sox common/p0/p0mw common/powrfl/powrfl common/index_rt/index_rt common/taumnx/taumn,taumx,pabstot,currtot c c read data plus initialization c index_rt=1 call prfile call paraminit call read_data call vectinit if(igrad.eq.0) call ic_rt if(igrad.gt.0) call ic_gb if(ierr.gt.0) then print*,' IERR = ', ierr stop end if c c beam/ray propagation call gray_integration c c postprocessing call after_gray_integration c pabstott=pabstot currtott=currtot c if (ipass.gt.1) then c second pass into plasma p0mw1=p0mw igrad=0 c index_rt=2 p0mw=p0mw1*powrfl call prfile call vectinit2 call paraminit call ic_rt2 call gray_integration call after_gray_integration pabstott=pabstott+pabstot currtott=currtott+currtot c index_rt=3 sox=-sox p0mw=p0mw1*(1.0_wp_-powrfl) call prfile call vectinit2 call paraminit call ic_rt2 call gray_integration call after_gray_integration pabstott=pabstott+pabstot currtott=currtott+currtot end if call vectfree print*,' ' write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0e3_wp_ c end c c c subroutine gray_integration use const_and_precisions, only : wp_ use graydata_flags, only : dst use beamdata, only : nstep implicit none c local variables integer :: i real(wp_) :: st0 c common/external functions/variables integer :: istep,istop,index_rt real(wp_) :: st,strfl11 c common/ss/st common/istep/istep common/istop/istop common/strfl11/strfl11 common/index_rt/index_rt c c ray integration: begin st0=0.0_wp_ if(index_rt.gt.1) st0=strfl11 do i=1,nstep istep=i st=i*dst+st0 c c advance one step call rkint4 c c calculations after one step: call after_onestep(istep,istop) if(istop.eq.1) exit c end do c c ray integration: end c return end c c c subroutine after_gray_integration use const_and_precisions, only : wp_,zero use graydata_flags, only : ibeam,iwarm,iequil,iprof, . nnd,ipec,filenmeqq,filenmprf,filenmbm use graydata_anequil, only : dens0,te0 use beamdata, only : nrayr implicit none c local variables integer :: iproj,nfilp,i real(wp_) :: currtka,pabs,currt real(wp_), dimension(nnd) :: dpdv,ajphiv,pins,currins real(wp_), dimension(nnd) :: dvol,darea,ratjav,ratjbv,ratjplv real(wp_), dimension(nnd) :: rt_in,rhop_tab,rhot_tab real(wp_), dimension(0:nnd) :: rtabpsi1 real(wp_) :: rhotpav,drhotpav real(wp_) :: rhotjava,drhotjava c common/external functions/variables integer :: index_rt real(wp_) :: st,taumn,taumx,pabstot,currtot ! arguments c common/ss/st common/index_rt/index_rt common/taumnx/taumn,taumx,pabstot,currtot c c print all ray positions in local reference system c if(nrayr.gt.1) then iproj=1 nfilp=9 call projxyzt(iproj,nfilp) end if c c print final results on screen c 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.0e3_wp_ write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka write(7,99) currtka,pabstot 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 ne0,Te0',dens0,te0 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_init(ipec,rt_in,rhop_tab,rhot_tab,rtabpsi1, . dvol,darea,ratjav,ratjbv,ratjplv) call spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv,pins,currins) call postproc_profiles(pabs,currt,rhot_tab,dvol,darea,dpdv,ajphiv, . rhotpav,drhotpav,rhotjava,drhotjava) write(*,*) 'rhotpav,drhotpav,rhotjava,drhotjava' write(*,99) rhotpav,drhotpav,rhotjava,drhotjava do i=1,nnd write(48,99) rhop_tab(i),rhot_tab(i),ajphiv(i), . ajphiv(i)*ratjbv(i),dpdv(i),currins(i),pins(i) end do return 99 format(16(1x,e16.8e3)) end c c c subroutine after_onestep(i,istop) use const_and_precisions, only : wp_,pi use reflections, only : inside,rlim,zlim,nlim,wall_refl use graydata_flags, only : iwarm,istpr0,istpl0,dst,ipass,igrad use graydata_par, only : psipol0,chipol0 use equilibrium, only : zbinf,zbsup use beamdata, only : nrayr,nrayth,psjki,ppabs,ccci,iiv,tauv,ihcd, . istore,iop,iow,tau1v,yyrfl,ywrk,ext,eyt implicit none c local constants real(wp_), parameter :: taucr=12.0_wp_ c arguments integer :: i, istop c local variables integer :: iproj,nfilp,irfl,j,k,kkk,iopmin,iowmax,iowmin,ivac,j1, . k1,kkkk real(wp_) :: qqout,uuout,vvout,qqin2,uuin2,vvin2,rrm,zzm,anvjkr, . aknmin,akdotn,anwclr real(wp_), dimension(6) :: y,dery real(wp_), dimension(3) :: xvrfl,anvrfl,xvvac,anw real(wp_), dimension(3,nrayr,nrayth) :: xvjk,anvjk complex(wp_) :: extr,eytr,exin2,eyin2 logical :: ins_pl c common/external functions/variables integer :: istpr,istpl,jclosest,ierr,index_rt real(wp_) :: psinv,psinv11,taumn,taumx,pabstot,currtot,psipol, . chipol,powrfl,strfl11 real(wp_), dimension(3) :: anwcl,xwcl,xv,anv logical :: inside_plasma c external inside_plasma c common/istgr/istpr,istpl common/refln/anwcl,xwcl,jclosest common/psival/psinv common/psinv11/psinv11 common/ierr/ierr common/taumnx/taumn,taumx,pabstot,currtot common/xv/xv common/anv/anv common/polcof/psipol,chipol common/powrfl/powrfl common/strfl11/strfl11 common/index_rt/index_rt c pabstot=0.0_wp_ currtot=0.0_wp_ taumn=1.0e+30_wp_ taumx=-1.0e+30_wp_ psinv11=1.0_wp_ iopmin=100 iowmin=100 iowmax=0 if(i.eq.1) then psipol=psipol0 chipol=chipol0 end if 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.0_wp_ zzm=xv(3)/100.0_wp_ if (iwarm.gt.0.and.i.gt.1.and.ihcd(j,k).gt.0) then if(psinv.ge.0.and.psinv.le.1.0_wp_.and. . zzm.ge.zbinf.and.zzm.le.zbsup) 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) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ins_pl=inside_plasma(rrm,zzm) if (mod(iop(j,k),2).eq.0 .and. ins_pl) then iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) if (ipass.gt.1 .and. index_rt.eq.1 .and. . iowmax.gt.1 .and. istore(j,k).eq.0) then istore(j,k)=istore(j,k)+1 yyrfl(j,k,1:3)=xv yyrfl(j,k,4:6)=anv ihcd(j,k)=0 end if else if (mod(iop(j,k),2).eq.1.and. . .not.ins_pl) then iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) end if if (ipass.gt.1) then if (iow(j,k).eq.0 .and. inside(rlim,zlim,nlim,rrm,zzm)) then iow(j,k)=1 else if (iow(j,k).eq.1 .and. . .not.inside(rlim,zlim,nlim,rrm,zzm)) then iow(j,k)=2 if (ins_pl) then iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) end if call wall_refl(xv-dst*anv,anv,ext(j,k,iop(j,k)), . eyt(j,k,iop(j,k)),xvrfl,anvrfl,extr,eytr,anw,irfl) if(index_rt.eq.1) then istore(j,k)=istore(j,k)+1 yyrfl(j,k,1:3)=xvrfl yyrfl(j,k,4:6)=anvrfl tau1v(j,k)=tauv(j,k,iiv(j,k)) ext(j,k,iop(j,k))=extr eyt(j,k,iop(j,k))=eytr end if if (j.lt.jclosest) then jclosest=j anwcl=anw xwcl=xvrfl end if xv=xvrfl anv=anvrfl rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2) zzm=1.0e-2_wp_*xv(3) ywrk(1:3,j,k)=xv ywrk(4:6,j,k)=anv igrad=0 call gwork(j,k) if (ins_pl) then iop(j,k)=iop(j,k)+1 call pol_limit(ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) if (index_rt.eq.1) ihcd(j,k)=0 end if end if end if if(index_rt.eq.1 .and. j.eq.1) psinv11=psinv if(iop(j,k).lt.iopmin) iopmin=iop(j,k) if(iow(j,k).lt.iowmin) iowmin=iow(j,k) if(iow(j,k).gt.iowmax) iowmax=iow(j,k) xvjk(:,j,k)=xv anvjk(:,j,k)=anv end do end do if(jclosest.le.nrayr) then aknmin=1.0_wp_ do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk anwclr=(anwcl(1)*xwcl(1)+anwcl(2)*xwcl(2)) . /sqrt(xwcl(1)**2+xwcl(2)**2) anvjkr=(anvjk(1,j,k)*xvjk(1,j,k)+anvjk(2,j,k)*xvjk(2,j,k)) . /sqrt(xvjk(1,j,k)**2+xvjk(2,j,k)**2) akdotn=anwclr*anvjkr+anwcl(3)*anvjk(3,j,k) if(akdotn.lt.aknmin) aknmin=akdotn end do end do else aknmin=-1.0_wp_ end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c print ray positions for j=nrayr in local reference system istpr=istpr+1 if (istpr.eq.istpr0) then c print*,'istep = ',i if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) end if istpr=0 end if c if (istpl.eq.istpl0) istpl=0 istpl=istpl+1 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c single pass is stopped when all the rays have crossed the plasma c or complete absorption has occurred c same for successive passes of multi-pass simulations (here exit c from vessel is detected too c first pass in multi-pass simulation is stopped when at least one c ray has reflected and all rays are directed away from c reflection point, or when no reflection has occurred and c central ray re-enters the plasma if((ipass.eq.1 .and. ((iopmin.gt.1) .or. . (taumn.lt.1.0e+30_wp_.and.taumn.gt.taucr))) . .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or. . (taumn.lt.1.0e+30_wp_.and.taumn.gt.taucr)))) then istop=1 else if(ipass.gt.1 .and. index_rt.eq.1 .and. . ((iowmax.gt.1 .and. aknmin.gt.0) .or. . (iowmax.le.1 .and. iop(1,1).gt.2))) then c flag second pass mode coupling as unset powrfl=-1.0_wp_ qqout=0.0_wp_ uuout=0.0_wp_ vvout=0.0_wp_ do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk c store missing initial conditions for the second pass if (istore(j,k).eq.0) then istore(j,k)=istore(j,k)+1 yyrfl(j,k,1:3)=xvjk(:,j,k) yyrfl(j,k,4:6)=anvjk(:,j,k) tau1v(j,k)=tauv(j,k,iiv(j,k)) end if c determine mode coupling at the plasma boundary if (powrfl.lt.0.0_wp_) then call vacuum_rt(xvjk(:,j,k),anvjk(:,j,k),xvvac,ivac) c look for first ray hitting the plasma, starting from the central c and evaluate polarization if (ivac.eq.1) then y(1:3)=xvvac y(4:6)=anvjk(:,j,k) call fwork(y,dery) call pol_limit(exin2,eyin2) call stokes(exin2,eyin2,qqin2,uuin2,vvin2) powloop: do j1=1,nrayr kkkk=nrayth if(j1.eq.1) kkkk=1 do k1=1,kkkk c look for first ray which completed the first pass in the plasma if (iop(j1,k1).gt.1) then c if found, use its polarization state to compute mode coupling call stokes(ext(j1,k1,2),eyt(j1,k1,2), . qqout,uuout,vvout) exit powloop end if end do end do powloop c if no ray completed a first pass in the plasma, use central ray c initial polarization (possibly reflected) if (qqout.le.0.0_wp_) then call stokes(ext(1,1,0),eyt(1,1,0),qqout,uuout,vvout) end if powrfl=0.5_wp_*(1.0_wp_+vvout*vvin2+ . uuout*uuin2+qqout*qqin2) end if end if end do end do if(powrfl.lt.0.0_wp_) powrfl=0.0_wp_ strfl11=i*dst write(6,*) ' ' write(6,*) 'Reflected power fraction =',powrfl write(66,*) psipol0,chipol0,powrfl istop=1 end if return end c c c subroutine print_output(i,j,k) use const_and_precisions, only : wp_,pi use graydata_flags, only : istpl0 use coreprofiles, only : psdbnd use beamdata, only : ywrk,psjki,tauv,alphav,pdjki, . currj,didst,q,tau1v,nrayr !,nrayth use equilibrium, only : frhotor implicit none c local constants integer, parameter :: ndim=6 real(wp_), parameter :: taucr=12.0_wp_ c arguments integer :: i,j,k c local variables real(wp_) :: x,y,z,rr,rrm,zzm,stm,xxm,yym,anx,any,anz,anr,anphi, . phi,phideg,psi,rhot,bbr,bbz,bpol,dens11,dids11,dpdv11,ajcd11, . alpha11,taujk,tau11,pt11 c common/external functions/variables integer :: nharm,nhf,iohkawa,index_rt,istpr,istpl real(wp_) :: p0mw,st,brr,bphi,bzz,ajphi,btot,xg,yg,dens,ddens, . tekev,alpha,effjcd,akim,tau0,anpl,anpr,ddr,an2s,an2,fdia,bdotgr, . ddi,ddr11,anprr,anpri complex(wp_) :: ex,ey,ez c common/nharm/nharm,nhf common/iokh/iohkawa common/p0/p0mw common/index_rt/index_rt common/ss/st common/istgr/istpr,istpl 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/epolar/ex,ey,ez common/nplr/anpl,anpr common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nprw/anprr,anpri 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.0e-2_wp_ zzm=z*1.0e-2_wp_ stm=st*1.0e-2_wp_ xxm=x*1.0e-2_wp_ yym=y*1.0e-2_wp_ 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.0_wp_) phi=-phi phideg=phi*180.0_wp_/pi psi=psjki(j,k,i) rhot=1.0_wp_ bbr=0.0_wp_ bbz=0.0_wp_ bpol=0.0_wp_ rhot=1.0_wp_ dens11=0.0_wp_ if(psi.ge.0.0_wp_) then if (psi.le.1.0_wp_) rhot=frhotor(psi) bbr=brr bbz=bzz bpol=sqrt(bbr**2+bbz**2) else tekev=0.0_wp_ akim=0.0_wp_ 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.0_wp_) dens11=dens dids11=didst(j,k,i)*1.0e2_wp_/(p0mw*q(j)) dpdv11=pdjki(j,k,i)/(p0mw*q(j)) ajcd11=currj(j,k,i)/(p0mw*q(j)) alpha11=alphav(j,k,i) tau11=taujk pt11=exp(-taujk) write(4,99) stm,rrm,zzm,phideg, . psi,rhot,dens11,tekev, . bbr,bphi,bbz,ajphi*1.0e-6_wp_,sqrt(an2),anpl,akim*100, . alpha11,tau11,pt11,ajcd11,dids11, . dble(nhf),dble(iohkawa),dble(index_rt) c call polarcold(exf,eyif,ezf,elf,etf) end if c central ray only end if(k.eq.1.and.j.eq.nrayr) write(17,99) st,ddr11,ddr,ddi c print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3 if(istpl.eq.istpl0) then if(j.eq.nrayr) then c if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then write(33,111) i,j,k,stm,xxm,yym,rrm,zzm, . psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt) end if c if(k.eq.nrayth) write(33,*) ' ' end if c return 99 format(30(1x,e16.8e3)) 111 format(3i5,16(1x,e16.8e3)) end c c c subroutine prfile implicit none c common/external functions/variables integer :: index_rt c common/index_rt/index_rt c if(index_rt.eq.1) then write(4,*)' #sst R z phi psi rhot ne Te BR Bphi Bz Jphi '// . 'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' write(8,*) ' #istep j k xt yt zt rt psin' write(9,*) ' #istep j k xt yt zt rt psin' write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' write(12,*) ' #i sst psi w1 w2' write(7,*)'#Icd Pa Jphip dPdVp '// . 'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// . 'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj '// . 'drhotp' write(48,*) '#rhop rhot Jphi Jcdb dPdV Icdins Pins' write(66,*) "# psipol0 chipol0 powrfl" c else c c if(index_rt.eq.3) then write(4,*) ' ' write(8,*) ' ' write(9,*) ' ' write(17,*) ' ' write(12,*) ' ' write(48,*) ' ' c end if end if c return end c c c subroutine read_data use const_and_precisions, only : wp_,qe=>ecgs_,me=>mecgs_, . vc=>ccgs_,cvdr=>degree,pi use graydata_flags use graydata_par use coreprofiles, only : psdbnd,read_profiles,tene_scal, . set_prfspl,set_prfan use equilibrium, only : read_eqdsk,change_cocos,eq_scal,set_eqspl, . rmxm,set_rhospl,setqphi_num,set_equian use reflections, only : rlim,zlim,nlim,alloc_lim use beamdata, only : nrayr,nrayth,nstep implicit none c local variables integer :: nfil,iox,ierr,isgniphi,isgnbphi,ifreefmt,idesc real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,w0csi,w0eta,d0csi, . d0eta,phi0,zrcsi,zreta,ssplf,rax,zax,rvac,psia,rr0m,zr0m,rpam,b0, . q0,qa,alq,te0,dte0,alt1,alt2,dens0,aln1,aln2,zeffan real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc, . rv,zv,fpol,qpsi,psinr,rhotn,rbbbs,zbbbs real(wp_), dimension(:,:), allocatable :: psin character(len=8) :: wdat character(len=10) :: wtim c common/external functions/variables real(wp_) :: xgcn,ak0,akinv,fhz,wcsi,weta,rcicsi,rcieta,phiw, . phir,anx0c,any0c,anz0c,x00,y00,z00,bres,p0mw,sox,alpha0,beta0 c common/xgcn/xgcn common/parwv/ak0,akinv,fhz common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/parbres/bres common/p0/p0mw common/mode/sox common/angles/alpha0,beta0 c open(2,file='gray.data',status= 'unknown') c c alpha0, beta0 (cartesian) launching angles c fghz wave frequency (GHz) c p0mw injected power (MW) c read(2,*) alpha0,beta0 read(2,*) fghz read(2,*) p0mw c c nrayr number of rays in radial direction c nrayth number of rays in angular direction c rwmax normalized maximum radius of beam power c rwmax=1 -> last ray at radius = waist c read(2,*) nrayr,nrayth,rwmax if(nrayr.eq.1) nrayth=1 c c x00,y00,z00 coordinates of launching point c read(2,*) x00,y00,z00 c c beams parameters in local reference system c w0 -> waist, d0 -> waist distance from launching point c phiw angle of beam ellipse c read(2,*) w0csi,w0eta,d0csi,d0eta,phiw c c ibeam=0 :read data for beam as above c ibeam=1 :read data from file simple astigmatic beam c ibeam=2 :read data from file general astigmatic beam read(2,*) ibeam read(2,*) filenmbm c c iequil=0 :vacuum c iequil=1 :analytical equilibrium c iequil=2 :read eqdsk c ixp=0,-1,+1 : no X point , bottom/up X point c read(2,*) iequil,ixp c c iprof=0 :analytical density and temp. profiles c iprof>0 :numerical density and temp. profiles c read(2,*) iprof c c iwarm=0 :no absorption and cd c iwarm=1 :weakly relativistic absorption c iwarm=2 :relativistic absorption, n<1 asymptotic expansion c iwarm=3 :relativistic absorption, numerical integration c ilarm :order of larmor expansion c imx :max n of iterations in dispersion, imx<0 uses 1st c iteration in case of failure after |imx| iterations read(2,*) iwarm,ilarm,imx c if(iwarm.gt.2) iwarm=2 c c ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models c read(2,*) ieccd !if (ieccd.ge.10) call Setup_SpitzFunc c c ipec=0/1 :pec profiles grid in psi/rhop c nnd :number of grid steps for pec profiles +1 c read(2,*) ipec,nnd c c dst integration step c nstep maximum number of integration steps c istpr0 projection step = dsdt*istprj c istpl0 plot step = dsdt*istpl c igrad=0 optical ray-tracing, initial conditions as for beam c igrad=1 quasi-optical ray-tracing c igrad=-1 ray-tracing, init. condit. c from center of mirror and with angular spread c ipass=1/2 1 or 2 passes into plasma c iox=1/2 OM/XM c idst=0/1/2 0 integration in s, 1 integr. in ct, 2 integr. in Sr c psipol0,chipol0 polarization angles c ipol=0 compute mode polarization at antenna, ipol=1 use polariz angles c read(2,*) dst,nstep,istpr0,istpl0,idst read(2,*) igrad,ipass,rwallm read(2,*) iox, psipol0,chipol0,ipol c c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004 c sspl spline parameter for psi interpolation c read(2,*) filenmeqq read(2,*) ipsinorm,sspl,ifreefmt,idesc c c factb factor for magnetic field (only for numerical equil) c scaling adopted: beta=const, qpsi=const, nustar=const c iscal ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling c factT factn factor for Te&ne scaling c read(2,*) factb, iscal,factt,factn if (factb.eq.0.0_wp_) factb=1.0_wp_ c c psbnd value of psi ( > 1 ) of density boundary c read(2,*) filenmprf read(2,*) psdbnd if(iprof.eq.0) psdbnd=1.0_wp_ c c signum of toroidal B and I c icocos index for equilibrium from COCOS - O. Sauter Feb 2012 read(2,*) isgnbphi,isgniphi,icocos c c read parameters for analytical density and temperature c profiles if iprof=0 c if(iprof.eq.0) then ! === to be replaced with call read_profiles_an === if (allocated(terad)) deallocate(terad) if (allocated(derad)) deallocate(derad) if (allocated(zfc)) deallocate(zfc) allocate(terad(4),derad(3),zfc(1)) read(2,*) derad(1:3) ! dens0,aln1,aln2 read(2,*) terad(1:4) ! te0,dte0,alt1,alt2 read(2,*) zfc(1) ! zeffan ! === end === call tene_scal(terad(1:2),derad(1:1),factt,factn,factb,iscal) call set_prfan(terad,derad,zfc) te0=terad(1) ! only for print in headers.txt dte0=terad(2) ! only for print in headers.txt alt1=terad(3) ! only for print in headers.txt alt2=terad(4) ! only for print in headers.txt dens0=derad(1) ! only for print in headers.txt aln1=derad(2) ! only for print in headers.txt aln2=derad(3) ! only for print in headers.txt zeffan=zfc(1) ! only for print in headers.txt deallocate(terad,derad,zfc) else read(2,*) dummy,dummy,dummy read(2,*) dummy,dummy,dummy,dummy read(2,*) dummy end if c c read parameters for analytical simple cylindical equilibrium c if iequil=0,1 if(iequil.lt.2) then read(2,*) rr0m,zr0m,rpam read(2,*) b0 read(2,*) q0,qa,alq b0=b0*factb call set_equian(rr0m,zr0m,rpam,b0,q0,qa,alq) call flux_average_an c call print_prof_an else read(2,*) dummy,dummy,dummy read(2,*) dummy read(2,*) dummy,dummy,dummy end if c close(unit=2) c if(nrayr.eq.1) igrad=0 if (nrayr.lt.5) then igrad=0 print*,' nrayr < 5 ! => OPTICAL CASE ONLY' print*,' ' end if c fhz=fghz*1.0e9_wp_ ak0=2.0e9_wp_*pi*fghz/vc akinv=1.0_wp_/ak0 c bresg=2.0_wp_*pi*fhz*me*vc/qe bres=bresg*1.0e-4_wp_ c c xg=xgcn*dens19 c xgcn=1.0e-5_wp_*qe**2/(pi*me*fghz**2) c sox=-1.0_wp_ if(iox.eq.2) sox=1.0_wp_ c c read data for beam from file if ibeam>0 c if(ibeam.gt.0) then call read_beams else zrcsi=0.5_wp_*ak0*w0csi**2 zreta=0.5_wp_*ak0*w0eta**2 if(igrad.gt.0) then wcsi=w0csi*sqrt(1.0_wp_+(d0csi/zrcsi)**2) weta=w0eta*sqrt(1.0_wp_+(d0eta/zreta)**2) rcicsi=-d0csi/(d0csi**2+zrcsi**2) rcieta=-d0eta/(d0eta**2+zreta**2) phir=phiw else d0eta=d0csi wcsi=w0csi*abs(d0csi/zrcsi) weta=w0eta*abs(d0eta/zreta) rcicsi=w0csi/zrcsi rcieta=w0eta/zreta if(d0csi.gt.0) then rcicsi=-rcicsi rcieta=-rcieta end if phir=phiw end if end if c print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0 c r00=sqrt(x00**2+y00**2) phi0=acos(x00/r00) if(y00.lt.0) phi0=-phi0 print'(a,4f8.3)','x00, y00, R00, z00 = ',x00,y00,r00,z00 print*,' ' c c anx0c=(anr0c*x00-anphi0c*y00)/r00 c any0c=(anr0c*y00+anphi0c*x00)/r00 c c anr0c=(anx0c*x00+any0c*y00)/r00 c anphi0c=(any0c*x00-anx0c*y00)/r00 c c angles alpha0, beta0 in a local reference system as proposed by Gribov et al c c anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0) c anphi0c=sin(cvdr*beta0) c anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0) anr0c=-cos(cvdr*beta0)*cos(cvdr*alpha0) anphi0c=sin(cvdr*beta0) anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0) anx0c=(anr0c*x00-anphi0c*y00)/r00 any0c=(anr0c*y00+anphi0c*x00)/r00 c c read data for Te , ne , Zeff from file if iprof>0 c if (iprof.eq.1) then nprof=98 call read_profiles(trim(filenmprf)//'.prf',psrad,terad,derad, . zfc,nprof) call tene_scal(terad,derad,factt,factn,factb,iscal) call set_prfspl(psrad,terad,derad,zfc,0.001_wp_) deallocate(psrad,terad,derad,zfc) end if c c read equilibrium data from file if iequil=2 c if (iequil.eq.2) then c neqdsk=99 call read_eqdsk(trim(filenmeqq)//'.eqdsk',rv,zv,psin,psia, . psinr,fpol,qpsi,rvac,rax,zax,rbbbs,zbbbs,rlim,zlim, . ipsinorm,idesc,ifreefmt,neqdsk) !,idesc,ifreefmt,neqdsk) call change_cocos(psia,fpol,qpsi,icocos,3) call eq_scal(psia,fpol,isgniphi,isgnbphi,factb) c ssplf=0.01_wp_ call set_eqspl(rv,zv,psin,psia,psinr,fpol,sspl,ssplf,rvac, . rax,zax,rbbbs,zbbbs,ixp) allocate(rhotn(size(qpsi))) call setqphi_num(psinr,abs(qpsi),abs(psia),rhotn) call set_rhospl(sqrt(psinr),rhotn) deallocate(rhotn) c nlim=size(rlim) call equidata(psinr,qpsi,size(qpsi)) c c locate btot=bres c call bfield_res(rv,zv,size(rv),size(zv)) deallocate(rv,zv,psin,fpol) c print density, temperature, safecty factor, toroidal current dens c versus psi, rhop, rhot call print_prof end if sgnbphi=isgnbphi sgnbphi=isgniphi if (iequil.eq.1) call bres_anal if (iequil.ne.2.or.ipass.lt.0) then c set simple limiter as two cylindrical walls at rwallm and r00 nlim=5 call alloc_lim(ierr) if (ierr.ne.0) stop rlim(1)=rwallm rlim(2)=max(r00*1.0e-2_wp_,rmxm) rlim(3)=rlim(2) rlim(4)=rlim(1) rlim(5)=rlim(1) zlim(1)=(z00-dst*nstep)*1.0e-2_wp_ zlim(2)=zlim(1) zlim(3)=(z00+dst*nstep)*1.0e-2_wp_ zlim(4)=zlim(3) zlim(5)=zlim(1) ipass=abs(ipass) end if nfil=78 open(nfil,file='headers.txt',status= 'unknown') call date_and_time(wdat,wtim) write(nfil,916) wdat(1:4),wdat(5:6),wdat(7:8), . wtim(1:2),wtim(3:4),wtim(5:6) write(nfil,904) REVISION if (iequil.eq.2) then write(nfil,905) trim(filenmeqq) else if (iequil.eq.1) write(nfil,912) rr0m,zr0m,rpam,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 ipol: ',7i5) 902 format('# X0 Y0 Z0 launching point (cm) : ',3(1x,es12.5)) 903 format('# w0xi w0eta d0xi d0eta (cm) phiw (deg) : ',5(1x,es12.5)) 904 format('# GRAY revision : ',a) 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)) end c c c subroutine surf_anal use const_and_precisions, only : wp_,pi use graydata_anequil, only : b0,rr0m,zr0m,rpam use magsurf_data, only : npsi,npoints,psicon,rcon,zcon, . alloc_surf_anal implicit none c local variables integer :: i,k,ierr real(wp_) :: rni,rres,zzres,zmx,zmn,dal,drn,drrm,dzzm,rrm,zzm c common/external functions/variables real(wp_) :: bres c common/parbres/bres c npsi=10 npoints=101 call alloc_surf_anal(ierr) if (ierr.ne.0) stop c c print circular magnetic surfaces iequil=1 c write(71,*) '#i psin R z' dal=2.0e-2_wp_*pi drn=0.1_wp_ do i=1,npsi rni=i*drn psicon(i)=rni do k=1,npoints drrm=rpam*rni*cos((k-1)*dal) dzzm=rpam*rni*sin((k-1)*dal) rrm=rr0m+drrm zzm=zr0m+dzzm rcon(i,k)=rrm zcon(i,k)=zzm write(71,111) i,rni,rrm,zzm end do write(71,*) ' ' write(71,*) ' ' end do c c print resonant B iequil=1 c write(70,*)'#i Btot R z' rres=b0*rr0m/bres zmx=zr0m+rpam zmn=zr0m-rpam do i=1,21 zzres=zmn+(i-1)*(zmx-zmn)/2.0e1_wp_ write(70,111) i,bres,rres,zzres end do return 111 format(i4,20(1x,e16.8e3)) end c c c subroutine bres_anal use const_and_precisions, only : wp_,pi use graydata_anequil, only : b0,rr0m,zr0m,rpam implicit none c local variables integer :: i real(wp_) :: rres,zmx,zmn,zzres c common/external functions/variables real(wp_) :: bres c common/parbres/bres c c print resonant B iequil=1 write(70,*)'#i Btot R z' rres=b0*rr0m/bres zmx=zr0m+rpam zmn=zr0m-rpam do i=1,21 zzres=zmn+(i-1)*(zmx-zmn)/2.0e1_wp_ write(70,111) i,bres,rres,zzres end do return 111 format(i4,20(1x,e16.8e3)) end c c c subroutine read_beams use const_and_precisions, only : wp_ use graydata_flags, only : filenmbm, ibeam use utils, only : locate use simplespline, only : difcs,spli implicit none c local constants integer, parameter :: nstrmx=50 c local variables integer :: i,k,ii,nfbeam,lbm,nisteer,steer,iopt,ier real(wp_) :: alphast,betast,x00mm,y00mm,z00mm,waist1,waist2, . dvir1,dvir2,delta,phi1,phi2,zr1,zr2,w1,w2,rci1,rci2,dal,betst real(wp_), dimension(nstrmx) :: alphastv,betastv,x00v,y00v, . z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v real(wp_), dimension(nstrmx,4) :: cbeta,cx0,cy0,cz0,cwaist1, . cwaist2,crci1,crci2,cphi1,cphi2 c common/external functions/variables real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir,x00,y00,z00, . alpha0,beta0,ak0,akinv,fhz c common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/mirr/x00,y00,z00 common/angles/alpha0,beta0 common/parwv/ak0,akinv,fhz c c for given alpha0 -> beta0 + beam parameters c c ibeam=1 simple astigmatic beam c ibeam=2 general astigmatic beam c c initial beam data are measured in mm -> transformed to cm c nfbeam=97 lbm=len_trim(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.5e-1_wp_*ak0*waist1**2 zr2=0.5e-1_wp_*ak0*waist2**2 w1=waist1*sqrt(1.0_wp_+(dvir1/zr1)**2) w2=waist2*sqrt(1.0_wp_+(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.1_wp_*x00mm y00v(i)=0.1_wp_*y00mm z00v(i)=0.1_wp_*z00mm alphastv(i)=alphast betastv(i)=betast waist1v(i)=0.1_wp_*w1 rci1v(i)=1.0e1_wp_*rci1 waist2v(i)=0.1_wp_*w2 rci2v(i)=1.0e1_wp_*rci2 phi1v(i)=phi1 phi2v(i)=phi2 end do c iopt=0 call difcs(alphastv,betastv,nisteer,iopt,cbeta,ier) call difcs(alphastv,waist1v,nisteer,iopt,cwaist1,ier) call difcs(alphastv,rci1v,nisteer,iopt,crci1,ier) call difcs(alphastv,waist2v,nisteer,iopt,cwaist2,ier) call difcs(alphastv,rci2v,nisteer,iopt,crci2,ier) call difcs(alphastv,phi1v,nisteer,iopt,cphi1,ier) call difcs(alphastv,phi2v,nisteer,iopt,cphi2,ier) call difcs(alphastv,x00v,nisteer,iopt,cx0,ier) call difcs(alphastv,y00v,nisteer,iopt,cy0,ier) call difcs(alphastv,z00v,nisteer,iopt,cz0,ier) c if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then call locate(alphastv,nisteer,alpha0,k) dal=alpha0-alphastv(k) betst=spli(cbeta,nisteer,k,dal) x00=spli(cx0,nisteer,k,dal) y00=spli(cy0,nisteer,k,dal) z00=spli(cz0,nisteer,k,dal) wcsi=spli(cwaist1,nisteer,k,dal) weta=spli(cwaist2,nisteer,k,dal) rcicsi=spli(crci1,nisteer,k,dal) rcieta=spli(crci2,nisteer,k,dal) phiw=spli(cphi1,nisteer,k,dal) phir=spli(cphi2,nisteer,k,dal) else print*,' alpha0 outside table range !!!' if(alpha0.ge.alphastv(nisteer)) ii=nisteer if(alpha0.le.alphastv(1)) ii=1 betst=betastv(ii) x00=x00v(ii) y00=y00v(ii) z00=z00v(ii) wcsi=waist1v(ii) weta=waist2v(ii) rcicsi=rci1v(ii) rcieta=rci2v(ii) phiw=phi1v(ii) phir=phi2v(ii) end if beta0=betst c close(nfbeam) c return end c c c subroutine equidata(psinq,qpsi,nq) use const_and_precisions, only : wp_,pi use utils, only : vmaxmini use graydata_par, only : factb use equilibrium, only : rup,zup,rlw,zlw,rmaxis,zmaxis, . zbinf,zbsup,frhotor implicit none c arguments integer, intent(in) :: nq real(wp_), dimension(nq), intent(in) :: psinq,qpsi c local variables integer :: i,info,iqmn,iqmx real(wp_) :: rbmin,rbmax,q2,q15,qmin,qmax, . rhot15,psi15,rhot2,psi2,rhotsx c compute flux surface averaged quantities rup=rmaxis rlw=rmaxis zup=zmaxis+(zbsup-zmaxis)/10.0_wp_ zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_ call flux_average c ipr=1 c call contours_psi(1.0_wp_,rcon,zcon,ipr) c do ii=1,2*np+1 c write(52,*) rcon(ii), zcon(ii) c end do c c locate psi surface for q=1.5 and q=2 rup=rmaxis rlw=rmaxis zup=(zbsup+zmaxis)/2.0_wp_ zlw=(zmaxis+zbinf)/2.0_wp_ q2=2.0_wp_ q15=1.5_wp_ call vmaxmini(abs(qpsi),nq,qmin,qmax,iqmn,iqmx) if (q15.gt.qmin.and.q15.lt.qmax) then call surfq(psinq,qpsi,nq,q15,psi15) rhot15=frhotor(sqrt(psi15)) print'(3(a,f8.5))','psi_1.5 = ',psi15,' rhop_1.5 = ' . ,sqrt(psi15),' rhot_1.5 = ',rhot15 end if if (q2.gt.qmin.and.q2.lt.qmax) then call surfq(psinq,qpsi,nq,q2,psi2) rhot2=frhotor(sqrt(psi2)) print'(3(a,f8.5))','psi_2 = ',psi2,' rhop_2 = ' . ,sqrt(psi2),' rhot_2 = ',rhot2 end if c return end c c c subroutine print_prof use const_and_precisions, only : wp_ use equilibrium, only : psinr,nq,fq,frhotor use coreprofiles, only : density, temp implicit none c local constants real(wp_), parameter :: eps=1.e-4_wp_ c local variables integer :: i real(wp_) :: psin,rhop,rhot,ajphi,te,qq real(wp_) :: dens,ddens c write(55,*) ' #psi rhot ne Te q Jphi' do i=1,nq psin=psinr(i) rhop=sqrt(psin) c call density(psin,dens,ddens) te=temp(psin) qq=fq(psin) rhot=frhotor(rhop) call tor_curr_psi(max(eps,psin),ajphi) write(55,"(12(1x,e12.5))") psin,rhot,dens,te,qq,ajphi*1.e-6_wp_ end do c return end c c c subroutine print_prof_an use const_and_precisions, only : wp_ use coreprofiles, only : density, temp use equilibrium, only : frhotor implicit none c local constants integer, parameter :: nst=51 c local variables integer :: i real(wp_) :: psin,rhop,rhot,te real(wp_) :: dens,ddens c write(55,*) ' #psi rhot ne Te' do i=1,nst psin=dble(i-1)/dble(nst-1) rhop=sqrt(psin) rhot=frhotor(rhop) call density(psin,dens,ddens) te=temp(psin) write(55,111) psin,rhot,dens,te end do c return 111 format(12(1x,e12.5)) end c c c subroutine surfq(psinq,qpsi,nq,qval,psival) use const_and_precisions, only : wp_ use magsurf_data, only : npoints use utils, only : locate, intlin implicit none c arguments integer, intent(in) :: nq real(wp_), dimension(nq), intent(in) :: psinq,qpsi real(wp_) :: qval,psival c local variables integer :: ncnt,i1,ipr real(wp_), dimension(npoints) :: rcn,zcn c ncnt=(npoints-1)/2 c c locate psi surface for q=qval c call locate(abs(qpsi),nq,qval,i1) call intlin(abs(qpsi(i1)),psinq(i1),abs(qpsi(i1+1)),psinq(i1+1), . qval,psival) ipr=1 call contours_psi(psival,rcn,zcn,ipr) c return end c c c subroutine bfield_res(rv,zv,nr,nz) use const_and_precisions, only : wp_ use equilibrium, only : bfield implicit none c arguments integer, intent(in) :: nr, nz real(wp_), intent(in) :: rv(nr), zv(nz) c local constants integer, parameter :: icmx=2002 c local variables integer :: j,k,n,nconts,inc,nctot integer, dimension(10) :: ncpts real(wp_) :: btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb real(wp_), dimension(icmx) :: rrcb,zzcb real(wp_), dimension(nr,nz) :: btotal c common/external functions/variables real(wp_) :: bres c common/parbres/bres c c Btotal on psi grid c btmx=-1.0e30_wp_ btmn=1.0e30_wp_ do j=1,nr rrj=rv(j) do k=1,nz zzk=zv(k) call bfield(rrj,zzk,bbphi,bbr,bbz) btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2) if(btotal(j,k).ge.btmx) btmx=btotal(j,k) if(btotal(j,k).le.btmn) btmn=btotal(j,k) c write(90,113) j,rrj,zzk,btotal(j,k) enddo c write(90,*) ' ' enddo c c compute Btot=Bres and Btot=Bres/2 c write(70,*)'#i Btot R z' do n=1,5 bbb=bres/dble(n) if (bbb.ge.btmn.and.bbb.le.btmx) then call cniteq(rv,zv,btotal,nr,nz,bbb, . nconts,ncpts,nctot,rrcb,zzcb) do inc=1,nctot write(70,113) inc,bbb,rrcb(inc),zzcb(inc) end do end if write(70,*) ' ' end do c return 113 format(i6,12(1x,e12.5)) end c c subroutine cniteq(rqgrid,zqgrid,btotal,nreq,nzeq,h, . ncon,npts,icount,rcon,zcon) c v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. c (based on an older code) c use const_and_precisions, only : wp_ implicit none c local constants integer, parameter :: icmx=2002 c arguments integer :: nreq,nzeq real(wp_) :: rqgrid(nreq),zqgrid(nzeq),btotal(nreq,nzeq) integer :: ncon,icount integer, dimension(10) :: npts real(wp_) :: h real(wp_), dimension(icmx) :: rcon,zcon c local variables integer :: i,j,k,l,ico,nrqmax,iclast,mpl,ix,jx, . mxr,n1,jm,jfor,lda,ldb,jabs,jnb,kx,ikx,itm,inext,in integer, dimension(3,2) :: ja integer, dimension(1000) :: lx real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y real(wp_), dimension(nreq*nzeq) :: a c data px/0.5d0/ c a=reshape(btotal,(/nreq*nzeq/)) c do ico=1,icmx rcon(ico)=0.0_wp_ zcon(ico)=0.0_wp_ enddo c nrqmax=nreq drgrd=rqgrid(2)-rqgrid(1) dzgrd=zqgrid(2)-zqgrid(1) 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,rcn,zcn,ipr) use const_and_precisions, only : wp_,pi use graydata_par, only : rwallm use magsurf_data, only : npoints use equilibrium, only : psiant,psinop,nsr,nsz, . cc=>cceq,tr,tz,rup,zup,rlw,zlw,kspl,points_tgo use dierckx, only : profil,sproota implicit none c local constants integer, parameter :: mest=4 c arguments integer :: ipr real(wp_) :: h real(wp_), dimension(npoints) :: rcn,zcn c local variables integer :: np,info,ic,ier,ii,iopt,m real(wp_) :: ra,rb,za,zb,th,zc,val real(wp_), dimension(mest) :: zeroc real(wp_), dimension(nsr) :: czc c np=(npoints-1)/2 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) c th=pi/dble(np) rcn(1)=rlw zcn(1)=zlw rcn(npoints)=rlw zcn(npoints)=zlw rcn(np+1)=rup zcn(np+1)=zup do ic=2,np zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_ iopt=1 call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) print*,' profil =',ier val=h*psiant+psinop call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) if (zeroc(1).gt.rwallm) then rcn(ic)=zeroc(1) zcn(ic)=zc rcn(npoints+1-ic)=zeroc(2) zcn(npoints+1-ic)=zc else rcn(ic)=zeroc(2) zcn(ic)=zc rcn(npoints+1-ic)=zeroc(3) zcn(npoints+1-ic)=zc end if end do if (ipr.gt.0) then do ii=1,npoints write(71,111) ii,h,rcn(ii),zcn(ii) end do write(71,*) write(71,*) end if c return 111 format(i6,12(1x,e12.5)) end c c c subroutine flux_average use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use magsurf_data use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge, . equinum_fpol,bfield,fq,frhotor use simplespline, only : difcs use dierckx, only : regrid,coeff_parder implicit none c local constants integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, . njest=nnintp+ksp+1,nlest=nlam+ksp+1, . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam c local variables integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr integer, dimension(kwrk) :: iwrk real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height,height2,r2iav,currp, . area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, . dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, . shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, . bphi,brr,bzz,riav,fp real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl real(wp_), dimension(2*ncnt) :: dlpv real(wp_), dimension(2*ncnt+1) :: bv,bpv real(wp_), dimension(nlam) :: alam,weights real(wp_), dimension(nnintp,nlam) :: fhlam real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam real(wp_), dimension(lwrk) :: wrk real(wp_), dimension(:), allocatable :: rctemp,zctemp c common/external functions/variables real(wp_) :: fpolv c npsi=nnintp ninpr=(npsi-1)/10 npoints = 2*ncnt+1 c call alloc_surfvec(ierr) if(allocated(tjp)) deallocate(tjp) if(allocated(tlm)) deallocate(tlm) if(allocated(ch)) deallocate(ch) allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), . rctemp(npoints),zctemp(npoints),stat=ierr) if (ierr.ne.0) return c computation of flux surface averaged quantities write(71,*)' #i psin R z' dlam=1.0_wp_/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam fhlam(1,l)=sqrt(1.0_wp_-alam(l)) ffhlam(l)=fhlam(1,l) dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) weights(l)=1.0_wp_ end do weights(1)=0.5_wp_ weights(nlam)=0.5_wp_ alam(nlam)=1.0_wp_ fhlam(1,nlam)=0.0_wp_ ffhlam(nlam)=0.0_wp_ dffhlam(nlam)=-99999.0_wp_ jp=1 anorm=2.0_wp_*pi*rmaxis/abs(btaxis) b2av=btaxis**2 dvdpsi=2.0_wp_*pi*anorm dadpsi=2.0_wp_*pi/abs(btaxis) ratio_cdator=abs(btaxis/btrcen) ratio_cdbtor=1.0_wp_ ratio_pltor=1.0_wp_ qq=1.0_wp_ fc=1.0_wp_ psicon(1)=0.0_wp_ rcon(1,:)=rmaxis zcon(1,:)=zmaxis pstab(1)=0.0_wp_ rhot_eq(1)=0.0_wp_ rpstab(1)=0.0_wp_ rhotqv(1)=0.0_wp_ vcurrp(1)=0.0_wp_ vajphiav(1)=0.0_wp_ bmxpsi(1)=abs(btaxis) bmnpsi(1)=abs(btaxis) bav(1)=abs(btaxis) rbav(1)=1.0_wp_ rri(1)=rmaxis varea(1)=0.0_wp_ vvol(1)=0.0_wp_ vratjpl(1)=ratio_pltor vratja(1)=ratio_cdator vratjb(1)=ratio_cdbtor ffc(1)=fc rhot2q=0.0_wp_ qqv(1)=1.0_wp_ dadrhotv(1)=0.0_wp_ dvdrhotv(1)=0.0_wp_ c rup=rmaxis c rlw=rmaxis c zup=(zbmax+zmaxis)/2.0_wp_ c zlw=(zmaxis+zbmin)/2.0_wp_ do jp=2,npsi height=dble(jp-1)/dble(npsi-1) psicon(jp)=height if(jp.eq.npsi) height=0.9999_wp_ ipr=0 jpr=mod(jp,ninpr) if(jpr.eq.1) ipr=1 height2=height*height call contours_psi(height2,rctemp,zctemp,ipr) rcon(jp,:) = rctemp zcon(jp,:) = zctemp c r2iav=0.0_wp_ anorm=0.0_wp_ dadpsi=0.0_wp_ currp=0.0_wp_ b2av=0.0_wp_ area=0.0_wp_ volume=0.0_wp_ ajphiav=0.0_wp_ bbav=0.0_wp_ bmmx=-1.0e+30_wp_ bmmn=1.0e+30_wp_ call tor_curr(rctemp(1),zctemp(1),ajphi0) call bfield(rctemp(1),zctemp(1),br=brr,bz=bzz) call equinum_fpol(height2,fpolv) bphi=fpolv/rctemp(1) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) bv(1)=btot0 bpv(1)=bpoloid0 rpsim0=rctemp(1) do inc=1,npoints-1 inc1=inc+1 dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ . (zctemp(inc1)-zctemp(inc))**2) drc=(rctemp(inc1)-rctemp(inc)) c compute length, area and volume defined by psi=height^2 ph=0.5_wp_*(dla+dlb+dlp) area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) area=area+sqrt(area2) rzp=rctemp(inc1)*zctemp(inc1) rz=rctemp(inc)*zctemp(inc) volume=pi*(rzp+rz)*drc+volume c compute line integral on the contour psi=height^2 rpsim=rctemp(inc1) zpsim=zctemp(inc1) call bfield(rpsim,zpsim,br=brr,bz=bzz) call tor_curr(rpsim,zpsim,ajphi) ! call equinum_fpol(height2,fpolv) bphi=fpolv/rpsim 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.5_wp_*dlp anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) dadpsi=dadpsi+dlph* . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(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.0_wp_/(bpoloid*rpsim**2) . +1.0_wp_/(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.0_wp_*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(height)/fq(height2) . *dadpsi/pi dvdrhotv(jp)=phitedge*frhotor(height)/fq(height2) . *dvdpsi/pi c computation of rhot from calculated q profile rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) . /dble(npsi-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.0_wp_ shlam=0.0_wp_ do l=nlam,1,-1 lam=alam(l) srl=0.0_wp_ rl2=1.0_wp_-lam*bv(1)/bmmx rl0=0.0_wp_ if(rl2.gt.0) rl0=sqrt(rl2) do inc=1,npoints-1 rl2=1.0_wp_-lam*bv(inc+1)/bmmx rl=0.0_wp_ if(rl2.gt.0) rl=sqrt(rl2) srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) rl0=rl end do srl=srl/anorm dhlam=0.5_wp_/srl fc=fc+lam/srl*weights(l) if(l.eq.nlam) then fhlam(jp,l)=0.0_wp_ ffhlam(nlam*(jp-1)+l)=0.0_wp_ dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam else shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam fhlam(jp,l)=shlam dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam end if end do fc=0.75_wp_*b2av/bmmx**2*fc*dlam ffc(jp)=fc ccfh=bmmn/bmmx/fc do l=1,nlam ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) end do end do write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// .' Area Vol |I_pl| qq fc ratioJa ratioJb' qqv(1)=qqv(2) vajphiav(1)=vajphiav(2) do jp=1,npsi rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) if(jp.eq.npsi) then rhotqv(jp)=1.0_wp_ rpstab(jp)=1.0_wp_ pstab(jp)=1.0_wp_ end if rhot_eq(jp)=frhotor(sqrt(pstab(jp))) write(56,99) pstab(jp),rhot_eq(jp),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 c rarea=sqrt(varea(npsi)/pi) c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi c used for computations of dP/dV and J_cd c spline coefficients of rhot iopt=0 call difcs(rpstab,vvol,npsi,iopt,cvol,ier) iopt=0 call difcs(rpstab,rbav,npsi,iopt,crbav,ier) iopt=0 call difcs(rpstab,rri,npsi,iopt,crri,ier) iopt=0 call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) iopt=0 call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) iopt=0 call difcs(rpstab,vratja,npsi,iopt,cratja,ier) iopt=0 call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) iopt=0 call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) iopt=0 call difcs(rpstab,varea,npsi,iopt,carea,ier) iopt=0 call difcs(rpstab,ffc,npsi,iopt,cfc,ier) iopt=0 call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) iopt=0 call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) iopt=0 call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) c spline interpolation of H(lambda,rhop) and dH/dlambda iopt=0 s=0.0_wp_ call regrid(iopt,npsi,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 return 99 format(20(1x,e12.5)) end c c c function fdadrhot(rpsi) use const_and_precisions, only : wp_ use magsurf_data, only : rpstab,cdadrhot,npsi use simplespline, only :spli implicit none c arguments real(wp_) :: rpsi,fdadrhot c local variables integer :: ip real(wp_) :: dps c ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.npsi) ip=npsi-1 ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) fdadrhot=spli(cdadrhot,npsi,ip,dps) return end c c c function fdvdrhot(rpsi) use const_and_precisions, only : wp_ use magsurf_data, only : rpstab,cdvdrhot,npsi use simplespline, only :spli implicit none c arguments real(wp_) :: rpsi,fdvdrhot c local variables integer :: ip real(wp_) :: dps c ip=int((npsi-1)*rpsi+1) ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) fdvdrhot=spli(cdvdrhot,npsi,ip,dps) c return end c c c subroutine flux_average_an use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use graydata_anequil, only : q0,qa,alq,b0,rr0m,zr0m,rpam use magsurf_data use equilibrium, only : btrcen, frhotor use simplespline, only : difcs use dierckx, only : regrid,coeff_parder implicit none c local constants integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, . njest=nnintp+ksp+1,nlest=nlam+ksp+1, . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam c local variables integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr integer, dimension(kwrk) :: iwrk real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height2,r2iav,currp,area, . volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb,dlp, . drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0,shlam, . srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi,bphi,brr,bzz,riav,fp, . btaxis,dadr,dbvcdc13,dbvcdc31,dvdr,height,phitedge,qqan,rmaxis, . zmaxis,rn real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl real(wp_), dimension(2*ncnt) :: dlpv real(wp_), dimension(2*ncnt+1) :: bv,bpv real(wp_), dimension(nlam) :: alam,weights real(wp_), dimension(nnintp,nlam) :: fhlam real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam real(wp_), dimension(lwrk) :: wrk real(wp_), dimension(:), allocatable :: rctemp,zctemp real(wp_) :: psinv !unused. can be removed when equian is modified !to accept optional arguments c common/external functions/variables real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv c common/fpol/fpolv,dfpv common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz c npsi=nnintp ninpr=(npsi-1)/10 npoints = 2*ncnt+1 c call alloc_surfvec(ierr) if(allocated(tjp)) deallocate(tjp) if(allocated(tlm)) deallocate(tlm) if(allocated(ch)) deallocate(ch) allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)), . rctemp(npoints),zctemp(npoints),stat=ierr) if (ierr.ne.0) return c c computation of flux surface averaged quantities rmaxis=rr0m zmaxis=zr0m btaxis=b0 call rhopol_an phitedge=pi*rpam*rpam*btaxis write(71,*)' #i psin R z' dlam=1.0_wp_/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam fhlam(1,l)=sqrt(1.0_wp_-alam(l)) ffhlam(l)=fhlam(1,l) dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) weights(l)=1.0_wp_ end do weights(1)=0.5_wp_ weights(nlam)=0.5_wp_ alam(nlam)=1.0_wp_ fhlam(1,nlam)=0.0_wp_ ffhlam(nlam)=0.0_wp_ dffhlam(nlam)=-99999.0_wp_ jp=1 anorm=2.0_wp_*pi*rmaxis/abs(btaxis) b2av=btaxis**2 dvdpsi=2.0_wp_*pi*anorm dadpsi=2.0_wp_*pi/abs(btaxis) ratio_cdator=abs(btaxis/btrcen) ratio_cdbtor=1.0_wp_ ratio_pltor=1.0_wp_ qq=1.0_wp_ fc=1.0_wp_ psicon(1)=0.0_wp_ rcon(1,:)=rmaxis zcon(1,:)=zmaxis pstab(1)=0.0_wp_ rhot_eq(1)=0.0_wp_ rpstab(1)=0.0_wp_ rhotqv(1)=0.0_wp_ vcurrp(1)=0.0_wp_ vajphiav(1)=0.0_wp_ bmxpsi(1)=abs(btaxis) bmnpsi(1)=abs(btaxis) bav(1)=abs(btaxis) rbav(1)=1.0_wp_ rri(1)=rmaxis varea(1)=0.0_wp_ vvol(1)=0.0_wp_ vratjpl(1)=ratio_pltor vratja(1)=ratio_cdator vratjb(1)=ratio_cdbtor ffc(1)=fc rhot2q=0.0_wp_ dadrhotv(1)=0.0_wp_ dvdrhotv(1)=0.0_wp_ qqv(1)=1.0_wp_ c rup=rmaxis c rlw=rmaxis c zup=(zbmax+zmaxis)/2.0_wp_ c zlw=(zmaxis+zbmin)/2.0_wp_ do jp=2,npsi height=dble(jp-1)/dble(npsi-1) psicon(jp)=height if(jp.eq.npsi) height=0.9999_wp_ height2=height*height ipr=0 jpr=mod(jp,ninpr) if(jpr.eq.1) ipr=1 call contours_psi_an(height2,rctemp,zctemp,ipr) rcon(jp,:) = rctemp zcon(jp,:) = zctemp c r2iav=0.0_wp_ anorm=0.0_wp_ dadpsi=0.0_wp_ currp=0.0_wp_ b2av=0.0_wp_ area=0.0_wp_ volume=0.0_wp_ ajphiav=0.0_wp_ bbav=0.0_wp_ bmmx=-1.0e+30_wp_ bmmn=1.0e+30_wp_ call equian(rctemp(1),zctemp(1),psinv) dbvcdc13=-ddpsidzz/rctemp(1) dbvcdc31= ddpsidrr/rctemp(1)-bzz/rctemp(1) ajphi=ccj*(dbvcdc13-dbvcdc31) brr=-dpsidz/rctemp(1) bzz= dpsidr/rctemp(1) bphi=fpolv/rctemp(1) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) bv(1)=btot0 bpv(1)=bpoloid0 rpsim0=rctemp(1) do inc=1,npoints-1 inc1=inc+1 dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2) dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2) dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+ . (zctemp(inc1)-zctemp(inc))**2) drc=(rctemp(inc1)-rctemp(inc)) c compute length, area and volume defined by psi=height^2 ph=0.5_wp_*(dla+dlb+dlp) area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) area=area+sqrt(area2) rzp=rctemp(inc1)*zctemp(inc1) rz=rctemp(inc)*zctemp(inc) volume=pi*(rzp+rz)*drc+volume c compute line integral on the contour psi=height^2 rpsim=rctemp(inc1) zpsim=zctemp(inc1) call equian(rpsim,zpsim,psinv) brr=-dpsidz/rpsim bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim ajphi=ccj*(dbvcdc13-dbvcdc31) bphi=fpolv/rpsim 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.5_wp_*dlp anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) dadpsi=dadpsi+dlph* . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(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.0_wp_/(bpoloid*rpsim**2) . +1.0_wp_/(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.0_wp_*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.0_wp_*pi*pi)) qqv(jp)=qq rn=frhotor(sqrt(pstab(jp))) qqan=q0+(qa-q0)*rn**alq dadr=2.0_wp_*pi*rn*rpam*rpam dvdr=4.0_wp_*pi*pi*rn*rmaxis*rpam*rpam c dadrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dadpsi c dvdrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dvdpsi dadrhotv(jp)=phitedge*rn*dadpsi/pi/qqan dvdrhotv(jp)=phitedge*rn*dvdpsi/pi/qqan c computation of rhot from calculated q profile rhot2q=rhot2q+(qqv(jp)*rpstab(jp)+qqv(jp-1)*rpstab(jp-1)) . /dble(npsi-1) c print*,jp,rhot2q,qqv(jp),rpstab(jp),qqv(jp-1),rpstab(jp-1) rhotqv(jp)=sqrt(rhot2q) c rhotqv(jp)=rn c write(57,99) rpstab(jp),rn,rhotqv(jp),qqv(jp),qqan,r2iav, . dadr,dadrhotv(jp),dadpsi,dvdpsi,fpolv 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.0_wp_ shlam=0.0_wp_ do l=nlam,1,-1 lam=alam(l) srl=0.0_wp_ rl2=1.0_wp_-lam*bv(1)/bmmx rl0=0.0_wp_ if(rl2.gt.0) rl0=sqrt(rl2) do inc=1,npoints-1 rl2=1.0_wp_-lam*bv(inc+1)/bmmx rl=0.0_wp_ if(rl2.gt.0) rl=sqrt(rl2) srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) rl0=rl end do srl=srl/anorm dhlam=0.5_wp_/srl fc=fc+lam/srl*weights(l) if(l.eq.nlam) then fhlam(jp,l)=0.0_wp_ ffhlam(nlam*(jp-1)+l)=0.0_wp_ dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam else shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam fhlam(jp,l)=shlam dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam end if end do fc=0.75_wp_*b2av/bmmx**2*fc*dlam ffc(jp)=fc ccfh=bmmn/bmmx/fc do l=1,nlam ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l) dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) end do end do write(56,*)' #psi rhot_eq rhot_av || |Bmx| |Bmn|'// .' Area Vol |I_pl| qq fc ratioJa ratioJb dadr dvdr' qqv(1)=qqv(2) vajphiav(1)=vajphiav(2) do jp=1,npsi rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) if(jp.eq.npsi) then rhotqv(jp)=1.0_wp_ rpstab(jp)=1.0_wp_ pstab(jp)=1.0_wp_ end if rhot_eq(jp)=frhotor(sqrt(pstab(jp))) write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), . vratja(jp),vratjb(jp),dadrhotv(jp),dvdrhotv(jp) end do c rarea=sqrt(varea(npsi)/pi) c spline coefficients of vvol,rbav,rri,bmxpsi,bmnpsi c used for computations of dP/dV and J_cd c spline coefficients of rhot iopt=0 call difcs(rpstab,vvol,npsi,iopt,cvol,ier) iopt=0 call difcs(rpstab,rbav,npsi,iopt,crbav,ier) iopt=0 call difcs(rpstab,rri,npsi,iopt,crri,ier) iopt=0 call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) iopt=0 call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) iopt=0 call difcs(rpstab,vratja,npsi,iopt,cratja,ier) iopt=0 call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) iopt=0 call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) iopt=0 call difcs(rpstab,varea,npsi,iopt,carea,ier) iopt=0 call difcs(rpstab,ffc,npsi,iopt,cfc,ier) iopt=0 call difcs(rpstab,rhotqv,npsi,iopt,crhotq,ier) iopt=0 call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) iopt=0 call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) c spline interpolation of H(lambda,rhop) iopt=0 s=0.0_wp_ call regrid(iopt,npsi,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 return 99 format(20(1x,e12.5)) end c c c subroutine rhopol_an use const_and_precisions, only : wp_ use dierckx, only : curfit,splev use graydata_par, only : sgniphi use graydata_anequil, only : q0,qa,alq,b0,rr0m,rpam use equilibrium, only : psia implicit none c local constants integer, parameter :: nnr=101,nrest=nnr+4,lwrk=nnr*4+nrest*16 c local variables integer :: i,ier,iopt,kspl integer, dimension(nrest) :: iwrk real(wp_) :: dr,xb,xe,ss,rp,fq0,fq1,qq,res,rn real(wp_), dimension(nnr) :: rhop,rhot,rhopi,rhoti,psin real(wp_), dimension(lwrk) :: wrk real(wp_), dimension(nrest) :: wp c common/external functions/variables integer :: nsrp,nsrot real(wp_), dimension(nrest) :: trp,crp,trot,crot c common/coffrtp/trp common/coffrn/nsrp common/coffrp/crp common/coffrptt/trot common/coffrnt/nsrot common/coffrpt/crot c dr=1.0_wp_/dble(nnr-1) rhot(1)=0.0_wp_ psin(1)=0.0_wp_ res=0.0_wp_ fq0=0.0_wp_ do i=2,nnr rhot(i)=(i-1)*dr rn=rhot(i) qq=q0+(qa-q0)*rn**alq fq1=rn/qq res=res+0.5_wp_*(fq1+fq0)*dr fq0=fq1 psin(i)=res end do wp=1.0_wp_ psin=psin/res rhop=sqrt(psin) psia=-sgniphi*abs(res*rpam*rpam*b0) print*,psia,log(8.0_wp_*rr0m/rpam)-2.0_wp_ wp(1)=1.0e3_wp_ wp(nnr)=1.0e3_wp_ c spline interpolation of rhopol versus rhotor iopt=0 xb=0.0_wp_ xe=1.0_wp_ ss=0.00001_wp_ kspl=3 call curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,nsrp, . trp,crp,rp,wrk,lwrk,iwrk,ier) call splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier) c spline interpolation of rhotor versus rhopol iopt=0 xb=0.0_wp_ xe=1.0_wp_ ss=0.00001_wp_ kspl=3 call curfit(iopt,nnr,rhop,rhot,wp,xb,xe,kspl,ss,nrest,nsrot, . trot,crot,rp,wrk,lwrk,iwrk,ier) call splev(trot,nsrot,crot,3,rhop,rhoti,nnr,ier) do i=1,nnr write(54,*) rhop(i),rhot(i),rhopi(i),rhoti(i) end do return end c c c subroutine contours_psi_an(h,rcn,zcn,ipr) use const_and_precisions, only : wp_,pi use graydata_anequil, only : rr0m,zr0m,rpam use magsurf_data, only : npoints use equilibrium, only : frhotor implicit none c arguments integer :: ipr real(wp_) :: h real(wp_), dimension(npoints) :: rcn,zcn c local variables integer :: np,ic real(wp_) :: rn,th c np=(npoints-1)/2 th=pi/dble(np) rn=frhotor(sqrt(h)) do ic=1,npoints zcn(ic)=zr0m+rpam*rn*sin(th*(ic-1)) rcn(ic)=rr0m+rpam*rn*cos(th*(ic-1)) c if (ipr.gt.0) write(71,111) ic,h,rcn(ic),zcn(ic) end do write(71,*) c return 111 format(i6,12(1x,e12.5)) end c c c subroutine vectinit use const_and_precisions, only : wp_ use dispersion, only: expinit use graydata_flags, only : iwarm use beamdata, only : jmx,kmx,nmx,psjki,tauv,alphav,pdjki,ppabs, . currj,didst,ccci,iiv,iop,iow,ihcd,istore,tau1v,nrayr,nrayth, . nstep,alloc_beam implicit none c local variables integer :: i,j,k,ierr c common/external functions/variables real(wp_) :: jclosest real(wp_), dimension(3) :: anwcl,xwcl c common/refln/anwcl,xwcl,jclosest c if(nstep.gt.nmx) nstep=nmx if(nrayr.gt.jmx) nrayr=jmx if(nrayth.gt.kmx) nrayth=kmx call alloc_beam(ierr) if (ierr.ne.0) return c jclosest=nrayr+1 anwcl(1:3)=0.0_wp_ xwcl(1:3)=0.0_wp_ c do i=1,nstep do k=1,nrayth do j=1,nrayr psjki(j,k,i)=0.0_wp_ tauv(j,k,i)=0.0_wp_ alphav(j,k,i)=0.0_wp_ pdjki(j,k,i)=0.0_wp_ ppabs(j,k,i)=0.0_wp_ didst(j,k,i)=0.0_wp_ ccci(j,k,i)=0.0_wp_ currj(j,k,i)=0.0_wp_ iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 ihcd(j,k)=1 istore(j,k)=0 tau1v(j,k)=0.0_wp_ end do end do end do c if(iwarm.gt.1) call expinit c return end c c c subroutine vectfree use beamdata, only : dealloc_beam implicit none c call dealloc_beam c return end c c c subroutine vectinit2 use const_and_precisions, only : wp_ use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj, . didst,ccci,iiv,iop,iow,ihcd,nrayr,nrayth,nstep implicit none c local variables integer :: i,j,k c do i=1,nstep do k=1,nrayth do j=1,nrayr psjki(j,k,i)=0.0_wp_ tauv(j,k,i)=0.0_wp_ alphav(j,k,i)=0.0_wp_ pdjki(j,k,i)=0.0_wp_ ppabs(j,k,i)=0.0_wp_ didst(j,k,i)=0.0_wp_ ccci(j,k,i)=0.0_wp_ currj(j,k,i)=0.0_wp_ iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 ihcd(j,k)=1 end do end do end do return end c c c subroutine paraminit use const_and_precisions, only : wp_ implicit none c common/external functions/variables integer :: istep,istpr,istpl,ierr,istop c common/istep/istep common/istgr/istpr,istpl common/ierr/ierr common/istop/istop c istpr=0 istpl=1 ierr=0 istep=0 istop=0 c return end c c c subroutine updatepos use const_and_precisions, only : wp_ use beamdata, only : xc,xco,du1,du1o,ywrk,nrayr,nrayth implicit none c local variables integer :: j,k 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 use const_and_precisions, only : wp_ use beamdata, only : dffiu,ddffiu,xc,xco,du1,du1o,gri,ggri, . dgrad2v,grad2,nrayr,nrayth implicit none c local variables integer :: iv,j,jm,jp,k,km,kp real(wp_) :: ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz, . dfu,dfuu,gx,gxx,gxy,gxz,gy,gyy,gyz,gz,gzz real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgu, . dgg1,dgg2,dgg3,df1,df2,df3 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.0_wp_ gri(2,j,k)=0.0_wp_ gri(3,j,k)=0.0_wp_ 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.0_wp_ uxz=(dgg1(3)+dgg3(1))/2.0_wp_ uyz=(dgg2(3)+dgg3(2))/2.0_wp_ 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.0_wp_*(gx*gxx+gy*gxy+gz*gxz) dgrad2v(2,j,k)=2.0_wp_*(gx*gxy+gy*gyy+gz*gyz) dgrad2v(3,j,k)=2.0_wp_*(gx*gxz+gy*gyz+gz*gzz) end do end do c return end c c c subroutine solg0(dxv1,dxv2,dxv3,dgg) c solution of the linear system of 3 eqs : dgg . dxv = dff c input vectors : dxv1, dxv2, dxv3, dff c output vector : dgg c dff=(1,0,0) c use const_and_precisions, only : wp_ implicit none c arguments real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgg c local variables real(wp_) :: denom,aa1,aa2,aa3 c 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 c return end c c c subroutine solg3(dxv1,dxv2,dxv3,df1,df2,df3,dg1,dg2,dg3) c three rhs vectors df1, df2, df3 c use const_and_precisions, only : wp_ implicit none c arguments real(wp_), dimension(3) :: dxv1,dxv2,dxv3,df1,df2,df3, . dg1,dg2,dg3 c local variables real(wp_) denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 c 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 c return end c c c subroutine rkint4 c Runge-Kutta integrator c use const_and_precisions, only : wp_ use graydata_flags, only : dst,igrad use beamdata, only : nrayr,nrayth,ywrk,ypwrk,grad2,dgrad2v, . gri,ggri implicit none c parameter integer, parameter :: ndim=6 c local variables integer :: ieq,iv,j,jv,k,kkk real(wp_) :: h,hh,h6 real(wp_), dimension(ndim) :: y,yy,fk1,fk2,fk3,fk4 c common/external functions/variables real(wp_) :: gr2 real(wp_), dimension(3) :: dgr2,dgr real(wp_), dimension(3,3) :: ddgr c common/gr/gr2 common/dgr/dgr2,dgr,ddgr c h=dst hh=h*0.5_wp_ h6=h/6.0_wp_ 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.0_wp_*fk2(ieq)+2.0_wp_*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) use const_and_precisions, only : wp_ use graydata_flags, only : igrad use beamdata, only : ywrk,ypwrk,grad2,dgrad2v,gri,ggri implicit none c local constants integer, parameter :: ndim=6 c arguments integer :: j,k c local variables integer :: ieq,iv,jv real(wp_), dimension(ndim) :: yy,yyp c common/external functions/variables real(wp_) :: gr2 real(wp_), dimension(3) :: dgr2,dgr real(wp_), dimension(3,3) :: ddgr c 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) use const_and_precisions, only : wp_ use graydata_flags, only : idst,igrad implicit none c local constants integer, parameter :: ndim=6 c arguments real(wp_), dimension(ndim) :: y,dery c local variables integer :: iv,jv real(wp_) :: xx,yy,zz,yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl, . dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel, . derdom,dfdiadnpl,dfdiadxg,dfdiadyg real(wp_), dimension(3) :: vgv,derdxv,danpldxv,derdnv,dbgr c common/external functions/variables integer :: ierr real(wp_) :: sox,gr2,dd,an2s,an2,fdia,bdotgr,ddi,ddr11,anpl, . anpr,xg,yg,vgm,derdnm,dersdst real(wp_), dimension(3) :: dgr2,dgr,xv,anv,bv,derxg,deryg real(wp_), dimension(3,3) :: ddgr,derbv 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/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 c ierr=0 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.0_wp_) 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.99_wp_) then if(abs(anpl).le.1.05_wp_) then ierr=97 else ierr=98 end if end if c anpl2=anpl*anpl dnl=1.0_wp_-anpl2 anpr2=an2-anpl2 anpr=0.0_wp_ if(anpr2.gt.0.0_wp_) anpr=sqrt(anpr2) yg2=yg**2 c an2s=1.0_wp_ dan2sdxg=0.0_wp_ dan2sdyg=0.0_wp_ dan2sdnpl=0.0_wp_ del=0.0_wp_ fdia=0.0_wp_ dfdiadnpl=0.0_wp_ dfdiadxg=0.0_wp_ dfdiadyg=0.0_wp_ c duh=1.0_wp_-xg-yg2 if(xg.gt.0.0_wp_) then del=sqrt(dnl*dnl+4.0_wp_*anpl*anpl*(1.0_wp_-xg)/yg2) an2s=1.0_wp_-xg -xg*yg2*(1.0_wp_+anpl2+sox*del)/duh/2.0_wp_ c dan2sdxg=-yg2*(1.0_wp_-yg2)*(1.0_wp_+anpl2+sox*del)/ . duh**2/2.0_wp_+sox*xg*anpl2/(del*duh)-1.0_wp_ dan2sdyg=-xg*yg*(1.0_wp_-xg)*(1.0_wp_+anpl2+sox*del)/duh**2 . +2.0_wp_*sox*xg*(1.0_wp_-xg)*anpl2/(yg*del*duh) dan2sdnpl=-xg*yg2*anpl/duh . -sox*xg*anpl*(2.0_wp_*(1.0_wp_-xg)-yg2*dnl)/(del*duh) c if(igrad.gt.0) then ddelnpl2=2.0_wp_*(2.0_wp_*(1.0_wp_-xg) . *(1.0_wp_+3.0_wp_*anpl2**2)-yg2*dnl**3)/yg2/del**3 fdia=-xg*yg2*(1.0_wp_+sox*ddelnpl2/2.0_wp_)/duh derdel=2.0_wp_*(1.0_wp_-xg)*anpl2*(1.0_wp_+3.0_wp_*anpl2**2) . -dnl**2*(1.0_wp_+3.0_wp_*anpl2)*yg2 derdel=4.0_wp_*derdel/(yg**5*del**5) ddelnpl2y=2.0_wp_*(1.0_wp_-xg)*derdel ddelnpl2x=yg*derdel dfdiadnpl=24.0_wp_*sox*xg*(1.0_wp_-xg)*anpl*(1.0_wp_-anpl**4) . /(yg2*del**5) dfdiadxg=-yg2*(1.0_wp_-yg2)/duh**2-sox*yg2*((1.0_wp_-yg2) . *ddelnpl2+xg*duh*ddelnpl2x)/(2.0_wp_*duh**2) dfdiadyg=-2.0_wp_*yg*xg*(1.0_wp_-xg)/duh**2 . -sox*xg*yg*(2.0_wp_*(1.0_wp_-xg)*ddelnpl2 . +yg*duh*ddelnpl2y)/(2.0_wp_*duh**2) c end if end if c bdotgr=0.0_wp_ do iv=1,3 bdotgr=bdotgr+bv(iv)*dgr(iv) dbgr(iv)=0.0_wp_ do jv=1,3 dbgr(iv)=dbgr(iv)+dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv) end do end do c derdnm=0.0_wp_ 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.5_wp_*bdotgr**2* . (derxg(iv)*dfdiadxg+deryg(iv)*dfdiadyg+danpldxv(iv)*dfdiadnpl) c derdnv(iv)=2.0_wp_*anv(iv)-bv(iv)*dan2sdnpl . +0.5_wp_*bdotgr**2*bv(iv)*dfdiadnpl c derdnm=derdnm+derdnv(iv)**2 end do c derdnm=sqrt(derdnm) c derdom=-2.0_wp_*an2+2.0_wp_*xg*dan2sdxg+yg*dan2sdyg+anpl*dan2sdnpl . +2.0_wp_*igrad*gr2 . -bdotgr**2*(fdia+xg*dfdiadxg+yg*dfdiadyg/2.0_wp_ . +anpl*dfdiadnpl/2.0_wp_) 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.5_wp_*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 const_and_precisions, only : wp_,pi,ccj=>mu0inv use graydata_flags, only : iequil use graydata_par, only : sgnbphi use equilibrium, only : equinum_fpol, equinum_psi implicit none c arguments real(wp_) :: xx,yy,zz c local variables integer :: iv,jv real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy, . rr,rr2,rrm,snphi,zzm real(wp_), dimension(3) :: dbtot,bvc real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv c common/external functions/variables real(wp_) :: bres,brr,bphi,bzz,ajphi,btot,xg,yg,dxgdpsi,dpsidr, . dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,psinv real(wp_), dimension(3) :: bv,derxg,deryg real(wp_), dimension(3,3) :: derbv 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/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/fpol/fpolv,dfpv common/psival/psinv xg=0.0_wp_ yg=9.9e1_wp_ c do iv=1,3 derxg(iv)=0.0_wp_ deryg(iv)=0.0_wp_ bv(iv)=0.0_wp_ dbtot(iv)=0.0_wp_ do jv=1,3 dbv(iv,jv)=0.0_wp_ derbv(iv,jv)=0.0_wp_ dbvcdc(iv,jv)=0.0_wp_ dbvcdc(iv,jv)=0.0_wp_ dbvdc(iv,jv)=0.0_wp_ 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.0e-2_wp_*zz rrm=1.0e-2_wp_*rr c if(iequil.eq.1) then call equian(rrm,zzm,psinv) end if c if(iequil.eq.2) then call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz, . ddpsidrz) call equinum_fpol(psinv,fpolv,dfpv) end if call sub_xg_derxg(psinv) yg=fpolv/(rrm*bres) bphi=fpolv/rrm btot=abs(bphi) if(psinv.lt.0.0_wp_) return 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 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)= dfpv*dpsidr/rrm-bphi/rrm dbvcdc(2,3)= dfpv*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.0e-2_wp_*dbtot(iv)/Bres bv(iv)=bv(iv)/btot do jv=1,3 derbv(iv,jv)=1.0e-2_wp_*(dbv(iv,jv)-bv(iv)*dbtot(jv))/btot end do end do c derxg(1)=1.0e-2_wp_*drrdx*dpsidr*dxgdpsi derxg(2)=1.0e-2_wp_*drrdy*dpsidr*dxgdpsi derxg(3)=1.0e-2_wp_*dpsidz*dxgdpsi c c current density computation in Ampere/m^2 c ccj==1/mu_0 c 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,psinv) use const_and_precisions, only : wp_ use graydata_par, only : sgnbphi,sgniphi use graydata_anequil, only : b0,rr0m,zr0m,rpam,q0,qa,alq use equilibrium, only : psia,frhopol implicit none c arguments real(wp_) :: rrm,zzm real(wp_) :: psinv c local variables real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,sgn,snt,rhop,rhot c common/external functions/variables real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv, . dfpv c common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/fpol/fpolv,dfpv c dpsidrp=0.0_wp_ d2psidrp=0.0_wp_ c c simple model for equilibrium: large aspect ratio c outside plasma: analytical continuation, not solution Maxwell eqs c rpm=sqrt((rrm-rr0m)**2+(zzm-zr0m)**2) rn=rpm/rpam c snt=0.0_wp_ cst=1.0_wp_ if (rpm.gt.0.0_wp_) then snt=(zzm-zr0m)/rpm cst=(rrm-rr0m)/rpm end if c rhot=rn if(rn.le.1) then rhop=frhopol(rhot) psinv=rhop*rhop else psinv=1.0_wp_+B0*rpam**2/2.0_wp_/abs(psia)/qa*(rn*rn-1.0_wp_) rhop=sqrt(psinv) end if c sgn=sgniphi*sgnbphi if(rn.le.1.0_wp_) then qq=q0+(qa-q0)*rn**alq dpsidrp=B0*rpam*rn/qq*sgn dqq=alq*(qa-q0)*rn**(alq-1.0_wp_) d2psidrp=B0*(1.0_wp_-rn*dqq/qq)/qq*sgn else dpsidrp=B0*rpam/qa*rn*sgn d2psidrp=B0/qa*sgn end if c fpolv=sgnbphi*b0*rr0m dfpv=0.0_wp_ 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 return end c c c subroutine tor_curr_psi(h,ajphi) use const_and_precisions, only : wp_ use equilibrium, only : zmaxis implicit none c arguments real(wp_) :: h,ajphi c local variables real(wp_) :: r1,r2 c call psi_raxis(h,r1,r2) call tor_curr(r2,zmaxis,ajphi) c return end c c c subroutine tor_curr(rpsim,zpsim,ajphi) use const_and_precisions, only : wp_,ccj=>mu0inv use equilibrium, only : equinum_psi implicit none c arguments real(wp_) :: rpsim,zpsim,ajphi c local variables real(wp_) :: bzz,dbvcdc13,dbvcdc31 real(wp_) :: dpsidr,ddpsidrr,ddpsidzz c call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,ddpsidrr=ddpsidrr, . ddpsidzz=ddpsidzz) bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim ajphi=ccj*(dbvcdc13-dbvcdc31) c return end c c c subroutine psi_raxis(h,r1,r2) use const_and_precisions, only : wp_ use equilibrium, only : psiant,psinop,zmaxis,nsr, . nsz,cc=>cceq,tr,tz,kspl use dierckx, only : profil,sproota implicit none c local constants integer, parameter :: mest=4 c arguments real(wp_) :: h,r1,r2 c local variables integer :: iopt,ier,m real(wp_) :: zc,val real(wp_), dimension(mest) :: zeroc real(wp_), dimension(nsr) :: czc c iopt=1 zc=zmaxis call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) print*,' profil =',ier val=h*psiant+psinop call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) r1=zeroc(1) r2=zeroc(2) c return end c c c subroutine sub_xg_derxg(psinv) use const_and_precisions, only : wp_ use equilibrium, only : psia use coreprofiles, only : density implicit none c arguments real(wp_) :: psinv c common/external functions/variables real(wp_) :: xg,xgcn,dxgdpsi real(wp_) :: dens,ddenspsin c common/xgxg/xg common/dxgdps/dxgdpsi common/xgcn/xgcn common/dens/dens,ddenspsin c xg=0.0_wp_ dxgdpsi=0.0_wp_ c if(psinv.le.psdbnd.and.psinv.ge.0) then call density(psinv,dens,ddenspsin) xg=xgcn*dens dxgdpsi=xgcn*ddenspsin/psia c end if c return end c c c subroutine ic_gb c beam tracing initial conditions igrad=1 c use const_and_precisions, only : wp_,izero,zero,one,pi, . cvdr=>degree,ui=>im use math, only : catand use graydata_flags, only : ipol use graydata_par, only : rwmax,psipol0,chipol0 use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk, . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt implicit none c local constants integer, parameter :: ndim=6,ndimm=3 c local variables integer :: j,k,iproj,nfilp real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw, . wwcsi,wweta,sk,sw,dk,dw,rci1,ww1,rci2,ww2,wwxx,wwyy,wwxy,rcixx, . rciyy,rcixy,dwwxx,dwwyy,dwwxy,d2wwxx,d2wwyy,d2wwxy,drcixx,drciyy, . drcixy,dr,da,ddfu,u,alfak,dcsiw,detaw,dx0t,dy0t,x0t,y0t,z0t,dx0, . dy0,dz0,x0,y0,z0,gxt,gyt,gzt,gr2,gxxt,gyyt,gzzt,gxyt,gxzt,gyzt, . dgr2xt,dgr2yt,dgr2zt,dgr2x,dgr2y,dgr2z,pppx,pppy,denpp,ppx,ppy, . anzt,anxt,anyt,anx,any,anz,an20,an0,anx0,any0,anz0,du1tx,du1ty, . du1tz,vgradi,r0,x0m,y0m,r0m,z0m,ddr110,deltapol,qq,uu,vv real(wp_), dimension(ndim) :: ytmp,yptmp complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1, . dqi2,dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy c common/external functions/variables real(wp_) :: ak0,akinv,fhz,wcsi,weta,rcicsi,rcieta,phiw,phir, . anx0c,any0c,anz0c,x00,y00,z00,dd,an2s,an2,fdia,bdotgr,ddi, . ddr11,psinv,dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi 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/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nplr/anpl,anpr common/psival/psinv common/parpl/brr,bphi,bzz,ajphi common/dens/dens,ddens common/tete/tekev c csth=anz0c snth=sqrt(1.0_wp_-csth**2) csps=1.0_wp_ snps=0.0_wp_ if(snth.gt.0.0_wp_) 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.0_wp_*akinv/wcsi**2 wweta=2.0_wp_*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.0_wp_*phirrad)-ui*dw*sin(2.0_wp_*phiwrad)) tc=(dk*cos(2.0_wp_*phirrad)-ui*dw*cos(2.0_wp_*phiwrad)) phic=0.5_wp_*catand(ts/tc) ddd=dk*cos(2*(phirrad+phic))-ui*dw*cos(2*(phiwrad+phic)) sss=sk-ui*sw qi1=0.5_wp_*(sss+ddd) qi2=0.5_wp_*(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.0_wp_*akinv/ww1) c z01=-rci1/(rci1**2+ww1**2) c w02=sqrt(2.0_wp_*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.0_wp_*phic) wwxx=-dimag(qqxx) wwyy=-dimag(qqyy) wwxy=-dimag(qqxy)/2.0_wp_ rcixx=dble(qqxx) rciyy=dble(qqyy) rcixy=dble(qqxy)/2.0_wp_ 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.0_wp_*phic) d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2 d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2 d2qqxy=-(d2qi1-d2qi2)*sin(2.0_wp_*phic) dwwxx=-dimag(dqqxx) dwwyy=-dimag(dqqyy) dwwxy=-dimag(dqqxy)/2.0_wp_ d2wwxx=-dimag(d2qqxx) d2wwyy=-dimag(d2qqyy) d2wwxy=-dimag(d2qqxy)/2.0_wp_ drcixx=dble(dqqxx) drciyy=dble(dqqyy) drcixy=dble(dqqxy)/2.0_wp_ dr=1.0_wp_ if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) da=2.0_wp_*pi/dble(nrayth) c ddfu=2.0_wp_*dr**2*akinv do j=1,nrayr u=dble(j-1) c ffi=u**2*ddfu/2.0_wp_ 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.5_wp_*(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.5_wp_*(x0t**2*dwwxx+y0t**2*dwwyy)+x0t*y0t*dwwxy gr2=gxt*gxt+gyt*gyt+gzt*gzt gxxt=wwxx gyyt=wwyy gzzt=0.5_wp_*(x0t**2*d2wwxx+y0t**2*d2wwyy)+x0t*y0t*d2wwxy gxyt=wwxy gxzt=x0t*dwwxx+y0t*dwwxy gyzt=x0t*dwwxy+y0t*dwwyy dgr2xt=2.0_wp_*(gxt*gxxt+gyt*gxyt+gzt*gxzt) dgr2yt=2.0_wp_*(gxt*gxyt+gyt*gyyt+gzt*gyzt) dgr2zt=2.0_wp_*(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.0_wp_*snth*csth*gyzt) . +2.0_wp_*snps*csps*(gxyt*csth+gxzt*snth) ggri(2,2,j,k)=gxxt*snps**2+csps**2 . *(gyyt*csth**2+gzzt*snth**2 . +2.0_wp_*snth*csth*gyzt) . -2.0_wp_*snps*csps*(gxyt*csth+gxzt*snth) ggri(3,3,j,k)=gzzt*csth**2+gyyt*snth**2 . -2.0_wp_*csth*snth*gyzt ggri(1,2,j,k)=csps*snps*(-gxxt+csth**2*gyyt+snth**2*gzzt . +2.0_wp_*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.5_wp_*u*(dx0t**2*dwwxx+dy0t**2*dwwyy . +2.0_wp_*dx0t*dy0t*dwwxy)/ddfu c pppx=x0t*rcixx+y0t*rcixy pppy=x0t*rcixy+y0t*rciyy denpp=pppx*gxt+pppy*gyt if (denpp.ne.0.0_wp_) then ppx=-pppx*gzt/denpp ppy=-pppy*gzt/denpp else ppx=0.0_wp_ ppy=0.0_wp_ end if c anzt=sqrt((1.0_wp_+gr2)/(1.0_wp_+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.0_wp_+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.0_wp_ ypwrk0(5,j,k) = dgr2y/an0/2.0_wp_ ypwrk0(6,j,k) = dgr2z/an0/2.0_wp_ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) if(ipol.eq.0) then call pol_limit(ext(j,k,0),eyt(j,k,0)) qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) call polellipse(qq,uu,vv,psipol0,chipol0) else qq=cos(2.0_wp_*chipol0*cvdr)*cos(2.0_wp_*psipol0*cvdr) uu=cos(2.0_wp_*chipol0*cvdr)*sin(2.0_wp_*psipol0*cvdr) vv=sin(2.0_wp_*chipol0*cvdr) if(qq**2.lt.1) then deltapol=asin(vv/sqrt(1.0_wp_-qq**2)) ext(j,k,0)= sqrt((1.0_wp_+qq)/2) eyt(j,k,0)= sqrt((1.0_wp_-qq)/2)*exp(-ui*deltapol) else ext(j,k,0)= 1.0_wp_ eyt(j,k,0)= 0.0_wp_ end if endif 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.0_wp_*vgradi c r0=sqrt(x0**2+y0**2) x0m=x0/1.0e2_wp_ y0m=y0/1.0e2_wp_ r0m=r0/1.0e2_wp_ z0m=z0/1.0e2_wp_ if(j.eq.nrayr) then write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, . psinv,one,dens,tekev,brr,bphi,bzz, . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one ddr110=dd end if if(j.eq.nrayr.and.k.eq.1) then write(17,99) zero,ddr110,dd,ddi end if end do end do call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end c c c subroutine ic_rt c ray tracing initial conditions igrad=0 c use const_and_precisions, only : wp_,izero,zero,one,pi, . cvdr=>degree,ui=>im use graydata_flags, only : ipol use graydata_par, only : rwmax,psipol0,chipol0 use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk, . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt implicit none c local constants integer, parameter :: ndim=6,ndimm=3 c local variables integer :: j,k,iv,jv,iproj,nfilp real(wp_) :: csth,snth,csps,snps,phiwrad,csphiw,snphiw,dr,da,u, . alfak,dcsiw,detaw,dx0t,dy0t,x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0, . anzt,anxt,anyt,anx,any,anz,an20,an0,anx0,any0,anz0,vgradi,r0, . x0m,y0m,r0m,z0m,ancsi,aneta,ppcsi,ppeta,deltapol,qq,uu,vv real(wp_), dimension(ndim) :: ytmp,yptmp c common/external functions/variables real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir,anx0c,any0c,anz0c, . x00,y00,z00,dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv,dens,ddens, . tekev,anpl,anpr,brr,bphi,bzz,ajphi c common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nplr/anpl,anpr common/psival/psinv common/parpl/brr,bphi,bzz,ajphi common/dens/dens,ddens common/tete/tekev c csth=anz0c snth=sqrt(1.0_wp_-csth**2) csps=1.0_wp_ snps=0.0_wp_ if(snth.gt.0.0_wp_) then csps=any0c/snth snps=anx0c/snth end if c phiwrad=phiw*cvdr csphiw=cos(phiwrad) snphiw=sin(phiwrad) c dr=1.0_wp_ if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) da=2.0_wp_*pi/dble(nrayth) z0t=0.0_wp_ c do j=1,nrayr u=dble(j-1) dffiu(j)=0.0_wp_ ddffiu(j)=0.0_wp_ 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.0_wp_/sqrt(1.0_wp_+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.0_wp_ 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.0_wp_ ypwrk0(5,j,k) = 0.0_wp_ ypwrk0(6,j,k) = 0.0_wp_ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) if(ipol.eq.0) then call pol_limit(ext(j,k,0),eyt(j,k,0)) qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) call polellipse(qq,uu,vv,psipol0,chipol0) else qq=cos(2.0_wp_*chipol0*cvdr)*cos(2.0_wp_*psipol0*cvdr) uu=cos(2.0_wp_*chipol0*cvdr)*sin(2.0_wp_*psipol0*cvdr) vv=sin(2.0_wp_*chipol0*cvdr) if(qq**2.lt.1.0_wp_) then c deltapol=phix-phiy, phix =0 deltapol=atan2(vv,uu) ext(j,k,0)= sqrt((1.0_wp_+qq)/2) eyt(j,k,0)= sqrt((1.0_wp_-qq)/2)*exp(-ui*deltapol) else if(qq.gt.0.0_wp_) then ext(j,k,0)= 1.0_wp_ eyt(j,k,0)= 0.0_wp_ else eyt(j,k,0)= 1.0_wp_ ext(j,k,0)= 0.0_wp_ end if end if endif c do iv=1,3 gri(iv,j,k)=0.0_wp_ dgrad2v(iv,j,k)=0.0_wp_ du10(iv,j,k)=0.0_wp_ do jv=1,3 ggri(iv,jv,j,k)=0.0_wp_ end do end do grad2(j,k)=0.0_wp_ c dd=anx0**2+any0**2+anz0**2-an20 vgradi=0.0_wp_ ddi=2.0_wp_*vgradi c r0=sqrt(x0**2+y0**2) x0m=x0/1.0e2_wp_ y0m=y0/1.0e2_wp_ r0m=r0/1.0e2_wp_ z0m=z0/1.0e2_wp_ if(j.eq.nrayr) then write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then write(17,99) zero,zero,zero,zero write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, . psinv,one,dens,tekev,brr,bphi,bzz, . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one end if end do end do call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end c c c subroutine ic_rt2 use const_and_precisions, only : wp_,izero,zero,one,pi, . cvdr=>degree use graydata_par, only : psipol0,chipol0 use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk, . xc0=>xc,du10=>du1,grad2,dgrad2v,gri,ggri,yyrfl,ext,eyt implicit none c local constants integer, parameter :: ndim=6,ndimm=3 c local variables integer :: j,k,iv,jv,iproj,nfilp real(wp_) :: x0,y0,z0,an20,an0,anx0,any0,anz0,vgradi, . r0,x0m,y0m,r0m,z0m,strfl11,qq,uu,vv real(wp_), dimension(ndim) :: ytmp,yptmp c common/external functions/variables real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv, . dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi c common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nplr/anpl,anpr common/psival/psinv common/parpl/brr,bphi,bzz,ajphi common/dens/dens,ddens common/tete/tekev c 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.0_wp_ ypwrk0(5,j,k) = 0.0_wp_ ypwrk0(6,j,k) = 0.0_wp_ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) call fwork(ytmp,yptmp) call pol_limit(ext(j,k,0),eyt(j,k,0)) qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) call polellipse(qq,uu,vv,psipol0,chipol0) c do iv=1,3 gri(iv,j,k)=0.0_wp_ dgrad2v(iv,j,k)=0.0_wp_ du10(iv,j,k)=0.0_wp_ do jv=1,3 ggri(iv,jv,j,k)=0.0_wp_ end do end do grad2(j,k)=0.0_wp_ c dd=anx0**2+any0**2+anz0**2-an20 vgradi=0.0_wp_ ddi=2.0_wp_*vgradi c r0=sqrt(x0**2+y0**2) x0m=x0/1.0e2_wp_ y0m=y0/1.0e2_wp_ r0m=r0/1.0e2_wp_ z0m=z0/1.0e2_wp_ if(j.eq.nrayr) then write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then write(17,99) zero,zero,zero,zero write(4,99) strfl11,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, . psinv,one,dens,tekev,brr,bphi,bzz, . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one end if end do end do c call pweigth c if(nrayr.gt.1) then iproj=0 nfilp=8 call projxyzt(iproj,nfilp) end if c return 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end c c c subroutine pweigth c ray power weigth coefficient q(j) c use const_and_precisions, only : wp_ use graydata_par, only : rwmax use beamdata, only : nrayr,nrayth,q implicit none c local variables integer :: j real(wp_) :: dr,r0,r,w,w0,summ c q(1)=1.0_wp_ if(nrayr.le.1) return dr=rwmax/dble(nrayr-1) summ=0.0_wp_ w0=1.0_wp_ do j=1,nrayr r=(dble(j)-0.5_wp_)*dr w=exp(-2.0_wp_*r**2) q(j)=w-w0 r0=r w0=w summ=summ+q(j) end do q=q/summ q(2:)=q(2:)/nrayth c return end c c c subroutine valpsisplpec(rhop,voli,areai) use const_and_precisions, only : wp_ use utils, only : locate use simplespline, only :spli use magsurf_data, only : rpstab,cvol,carea,npsi implicit none c arguments real(wp_), intent(in) :: rhop real(wp_), intent(out) :: voli,areai c local variables integer :: ip real(wp_) :: drh c call locate(rpstab,npsi,rhop,ip) ip=min(max(1,ip),npsi-1) drh=rhop-rpstab(ip) c voli=spli(cvol,npsi,ip,drh) areai=spli(carea,npsi,ip,drh) c return end c c c subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) use const_and_precisions, only : wp_ use magsurf_data, only : rpstab,cratja,cratjb,cratjpl,npsi use utils, only : locate use simplespline, only :spli implicit none c arguments real(wp_) :: rpsi,ratjai,ratjbi,ratjpli c local variables integer :: ip real(wp_) :: dps c call locate(rpstab,npsi,rpsi,ip) ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) ratjai=spli(cratja,npsi,ip,dps) ratjbi=spli(cratjb,npsi,ip,dps) ratjpli=spli(cratjpl,npsi,ip,dps) c return end c c c subroutine pabs_curr(i,j,k) use const_and_precisions, only : wp_,pi,mc2=>mc2_ use graydata_flags, only : dst,iwarm,ilarm,ieccd,imx use graydata_par, only : sgnbphi use graydata_anequil, only : rr0m use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,didst, . ccci,q,tau1v use coreprofiles, only : temp, fzeff use dispersion, only : harmnumber, warmdisp use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl implicit none c arguments integer, intent(in) :: i,j,k c local constants real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ c local variables real(wp_) :: dvvl,rbavi,rrii,rhop0,rhop,dervoli, . tau,effjcdav,dpdst,p0,pow,alpha0,didst0,ccci0 integer :: lrm,ithn real(wp_) :: amu,ratiovgr,rbn,rbx,zeff real(wp_) :: cst2,bmxi,bmni,fci real(wp_), dimension(:), allocatable :: eccdpar c common/external functions/variables integer :: nhmin,nhmax,iokhawa,ierr real(wp_) :: p0mw,ak0,akinv,fhz,dersdst real(wp_) :: psinv,xg,yg,tekev,dens,ddens,btot real(wp_) :: anpl,anpr,vgm,derdnm,sox,anprre,anprim,alpha,effjcd, . akim,tau0 complex(wp_) :: ex,ey,ez c common/p0/p0mw common/dersdst/dersdst common/absor/alpha,effjcd,akim,tau0 ! solo per print_output common/psival/psinv common/tete/tekev ! per print_output common/dens/dens,ddens common/btot/btot common/nharm/nhmin,nhmax common/xgxg/xg common/ygyg/yg common/nplr/anpl,anpr common/vgv/vgm,derdnm common/parwv/ak0,akinv,fhz common/mode/sox common/nprw/anprre,anprim common/epolar/ex,ey,ez common/iokh/iokhawa common/ierr/ierr c dvvl=1.0_wp_ rbavi=1.0_wp_ rrii=rr0m rhop=sqrt(psinv) c 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 rhop0=sqrt(psjki(j,k,i-1)) alpha0=alphav(j,k,i-1) tau0=tauv(j,k,i-1) didst0=didst(j,k,i-1) ccci0=ccci(j,k,i-1) c c c absorption computation c alpha=0.0_wp_ akim=0.0_wp_ effjcd=0.0_wp_ c tekev=temp(psinv) call valpsisplhcd(rhop,rrii,rbavi,dervoli,bmni,bmxi,fci) c if(tekev.gt.0.0_wp_.and.tau0.le.taucr) then c amu=mc2/tekev call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) if(nhmin.gt.0) then c lrm=max(ilarm,nhmax) call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim, * iwarm,imx,ex,ey,ez) c akim=ak0*anprim ratiovgr=2.0_wp_*anpr/derdnm!*vgm alpha=2.0_wp_*akim*ratiovgr if(alpha.lt.0.0_wp_) then ierr=94 print*,' IERR = ', ierr,' alpha negative', alpha end if c if(ieccd.gt.0) then ithn=1 if(lrm.gt.nhmin) ithn=2 zeff=fzeff(psinv) rbn=btot/bmni rbx=btot/bmxi select case(ieccd) case(1) c cohen model call setcdcoeff(zeff,rbn,rbx,cst2,eccdpar) call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, . ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierr) case(2) c no trapping call setcdcoeff(zeff,cst2,eccdpar) call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, . ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierr) case default c neoclassical model call setcdcoeff(zeff,rbx,fci,amu,rhop,cst2,eccdpar) call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, . ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierr) end select !deallocate(eccdpar) end if end if end if c alphav(j,k,i)=alpha tau=tau0+0.5_wp_*(alpha+alpha0)*dersdst*dst tauv(j,k,i)=tau p0=p0mw*q(j)*exp(-tau1v(j,k)) pow=p0*exp(-tau) ppabs(j,k,i)=p0-pow dpdst=pow*alpha*dersdst dvvl=dervoli*abs(rhop-rhop0) if(dvvl.le.0.0_wp_) dvvl=1.0_wp_ pdjki(j,k,i)=dpdst*dst/dvvl effjcdav=rbavi*effjcd currj(j,k,i)=sgnbphi*effjcdav*pdjki(j,k,i) didst(j,k,i)=sgnbphi*effjcdav*dpdst/(2.0_wp_*pi*rrii) ccci(j,k,i)=ccci0+0.5_wp_*(didst(j,k,i)+didst0)*dst return end c c c subroutine valpsisplhcd(rhop,rrii,rbavi,dervoli,bmni,bmxi,fci) use const_and_precisions, only : wp_ use utils, only : locate use simplespline, only :spli,splid use magsurf_data, only : rpstab,crri,crbav,cvol,cbmn,cbmx,cfc,npsi implicit none c arguments real(wp_), intent(in) :: rhop real(wp_), intent(out) :: rrii,rbavi,dervoli,bmni,bmxi,fci c local variables integer :: ip real(wp_) :: drh c call locate(rpstab,npsi,rhop,ip) ip=min(max(1,ip),npsi-1) drh=rhop-rpstab(ip) rrii=spli(crri,npsi,ip,drh) rbavi=spli(crbav,npsi,ip,drh) dervoli=splid(cvol,npsi,ip,drh) bmni=spli(cbmn,npsi,ip,drh) bmxi=spli(cbmx,npsi,ip,drh) fci=spli(cfc,npsi,ip,drh) c return end c c c subroutine projxyzt(iproj,nfile) use const_and_precisions, only : wp_ use beamdata, only : ywrk,nrayr,nrayth implicit none c arguments integer :: iproj,nfile c local variables integer :: jd,j,kkk,k real(wp_) :: dir,rtimn,rtimx,dx,dy,dz,dirx,diry,dirz, . csth1,snth1,csps1,snps1,xti,yti,zti,rti c common/external functions/variables integer :: istep real(wp_) :: psinv11,st c common/psinv11/psinv11 common/istep/istep common/ss/st c rtimn=1.0e+30_wp_ rtimx=-1.0e-30_wp_ 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.0_wp_-csth1**2) csps1=1.0_wp_ snps1=0.0_wp_ if(snth1.gt.0.0_wp_) 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) Cc c dirxt= (dirx*csps1-diry*snps1)/dir c diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir c dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir Cc c if(k.eq.1) then c xti1=xti c yti1=yti c zti1=zti c rti1=rti c end if if(.not.(iproj.eq.0.and.j.eq.1)) . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 c if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti c end do c c if(.not.(iproj.eq.0.and.j.eq.1)) c . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 if(iproj.eq.1) write(nfile,*) ' ' end do c write(nfile,*) ' ' c write(12,99) istep,st,psinv11,rtimn,rtimx return 99 format(i5,12(1x,e16.8e3)) 111 format(3i5,12(1x,e16.8e3)) end subroutine pec_init(ipec,rt_in,rhop_tab,rhot_tab,rtabpsi1, . dvol,darea,ratjav,ratjbv,ratjplv) use const_and_precisions, only : wp_,zero,one,izero use graydata_flags, only : nnd use equilibrium, only : frhotor,frhopol implicit none integer :: it,ipec real(wp_) :: drt,rt,rt1,rhop1 real(wp_) :: ratjai,ratjbi,ratjpli real(wp_) :: voli0,voli1,areai0,areai1 real(wp_), dimension(nnd), intent(in) :: rt_in real(wp_), dimension(nnd), intent(out) :: rhop_tab,rhot_tab real(wp_), dimension(0:nnd), intent(out) :: rtabpsi1 real(wp_), dimension(nnd), intent(out) :: dvol,darea real(wp_), dimension(nnd), intent(out) :: ratjav,ratjbv,ratjplv c ipec positive build equidistant grid dimension nnd c ipec negative read input grid c c ipec=+/-1 rho_pol grid c ipec=+/-2 rho_tor grid voli0=zero areai0=zero rtabpsi1(0)=zero do it=1,nnd if(ipec > 0) then c build equidistant radial grid drt=one/dble(nnd-1) rt=dble(it-1)*drt else c read radial grid from input rt=rt_in(it) drt=(rt_in(it+1)-rt)/2.0_wp_ end if c radial coordinate of i-(i+1) interval mid point if(it < nnd) then rt1=rt+drt/2.0_wp_ else rt1=one end if if (abs(ipec) == 1) then rhop_tab(it)=rt rhot_tab(it)=frhotor(rt) rhop1=rt1 else rhot_tab(it)=rt rhop_tab(it)=frhopol(rt) rhop1=frhopol(rt1) end if c psi grid at mid points, dimension nnd+1, for use in pec_tab rtabpsi1(it)=rhop1**2 call valpsisplpec(rhop1,voli1,areai1) dvol(it)=abs(voli1-voli0) darea(it)=abs(areai1-areai0) voli0=voli1 areai0=areai1 call ratioj(rhop_tab(it),ratjai,ratjbi,ratjpli) ratjav(it)=ratjai ratjbv(it)=ratjbi ratjplv(it)=ratjpli end do end subroutine pec_init subroutine spec(pabs,currt,rtabpsi1,dvol,darea,dpdv,ajphiv, . pins,currins) use const_and_precisions, only : wp_,one,zero,izero use graydata_flags, only : nnd use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv implicit none c local constants real(wp_), parameter :: rtbc=one c arguments real(wp_), intent(in) :: pabs,currt real(wp_), dimension(nstep):: xxi,ypt,yamp real(wp_), dimension(nnd), intent(out) :: dpdv,ajphiv real(wp_), dimension(nnd), intent(out) :: pins,currins real(wp_), dimension(nnd), intent(in) :: dvol,darea real(wp_), dimension(0:nnd), intent(in) :: rtabpsi1 c local variables integer :: i,ii,j,k,kkk real(wp_) :: spds,sccs,facpds,facjs real(wp_), dimension(nnd) :: wdpdv,wajphiv c calculation of dP and dI over radial grid dpdv=zero ajphiv=zero kkk=1 do j=1,nrayr if(j > 1) kkk=nrayth do k=1,kkk ii=iiv(j,k) if (ii < nstep ) then if(psjki(j,k,ii+1) /= zero) ii=ii+1 end if xxi=zero ypt=zero yamp=zero do i=1,ii xxi(i)=abs(psjki(j,k,i)) if(xxi(i) <= one) then ypt(i)=ppabs(j,k,i) yamp(i)=ccci(j,k,i) end if end do call pec_tab(xxi,ypt,yamp,ii,rtabpsi1, . wdpdv,wajphiv) dpdv=dpdv+wdpdv ajphiv=ajphiv+wajphiv end do end do c here dpdv is still dP (not divided yet by dV) c here ajphiv is still dI (not divided yet by dA) spds=zero sccs=zero do i=1,nnd spds=spds+dpdv(i) sccs=sccs+ajphiv(i) pins(i)=spds currins(i)=sccs dpdv(i)=dpdv(i)/dvol(i) ajphiv(i)=ajphiv(i)/darea(i) end do facpds=one facjs=one if(spds > zero) facpds=pabs/spds if(sccs /= zero) facjs=currt/sccs do i=1,nnd dpdv(i)=facpds*dpdv(i) ajphiv(i)=facjs*ajphiv(i) end do c now dpdv is dP/dV [MW/m^3] c now ajphiv is J_phi=dI/dA [MA/m^2] end subroutine spec subroutine pec_tab(xxi,ypt,yamp,ii,rtabpsi1,wdpdv,wajphiv) c Power and current projected on psi grid - mid points use const_and_precisions, only : wp_,one,zero use graydata_flags, only : nnd use beamdata, only : nstep use utils, only : locatex,intlin c arguments integer, intent(in) :: ii integer, parameter :: llmx = 21 integer, dimension(llmx) ::isev real(wp_), dimension(nstep), intent(in) :: xxi,ypt,yamp real(wp_), dimension(0:nnd), intent(in) :: rtabpsi1 real(wp_), dimension(nnd), intent(out) :: wdpdv,wajphiv c local variables real(wp_) :: ppa1,ppa2,cci1,cci2,dppa,didst,rt1 integer :: i,is,ise0,idecr,iise0,iise,iis,iis1 integer :: ind1,ind2,iind,ind,indi,itb1 isev=0 ise0=0 idecr=-1 is=1 wdpdv=zero wajphiv=zero do i=1,ii if(ise0 == 0) then if(xxi(i) < one) then ise0=i isev(is)=i-1 is=is+1 end if else if (idecr == -1) then if(xxi(i) > xxi(i-1)) then isev(is)=i-1 is=is+1 idecr=1 end if else if(xxi(i) > one) exit if(xxi(i) < xxi(i-1)) then isev(is)=i-1 is=is+1 idecr=-1 end if end if end if end do isev(is)=i-1 ppa1=zero cci1=zero do iis=1,is-1 iis1=iis+1 iise0=isev(iis) iise=isev(iis1) if (mod(iis,2) /= 0) then idecr=-1 ind1=nnd ind2=2 iind=-1 else idecr=1 ind1=1 ind2=nnd iind=1 end if do ind=ind1,ind2,iind indi=ind if (idecr == -1) indi=ind-1 rt1=rtabpsi1(indi) call locatex(xxi,iise,iise0,iise,rt1,itb1) if(itb1 >= iise0 .and. itb1 < 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 wdpdv(ind)=wdpdv(ind)+dppa didst=cci2-cci1 wajphiv(ind)=wajphiv(ind)+didst ppa1=ppa2 cci1=cci2 end if end do end do end subroutine pec_tab subroutine postproc_profiles(pabs,currt,rhot_tab,dvol,darea, . dpdv,ajphiv,rhotpav,drhotpav,rhotjava,drhotjava) c radial average values over power and current density profile use const_and_precisions, only : wp_,one,zero,izero,pi use graydata_flags, only : nnd,ieccd use equilibrium, only : frhopol implicit none real(wp_), intent(in) :: pabs,currt real(wp_), intent(out) :: rhotpav,rhotjava real(wp_), intent(out) :: drhotpav,drhotjava real(wp_) :: rhopjava,rhoppav real(wp_) :: dpdvp,dpdvmx,rhotp,drhotp real(wp_) :: ajphip,ajmxfi,rhotjfi,drhotjfi real(wp_) :: ratjamx,ratjbmx,ratjplmx real(wp_), dimension(nnd), intent(in) :: rhot_tab,dvol,darea real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv real(wp_), external :: fdadrhot,fdvdrhot real(wp_) :: sccsa real(wp_) :: rhotjav,rhot2pav,rhot2java rhotpav=zero rhot2pav=zero rhotjav=zero rhotjava=zero rhot2java=zero if (pabs > zero) then rhotpav=sum(rhot_tab(1:nnd)*dpdv(1:nnd)*dvol(1:nnd))/pabs rhot2pav=sum(rhot_tab(1:nnd)*rhot_tab(1:nnd)* . dpdv(1:nnd)*dvol(1:nnd))/pabs end if if (ieccd /= 0) then if (abs(currt) > zero) then rhotjav=sum(rhot_tab(1:nnd)*ajphiv(1:nnd)*darea(1:nnd))/currt end if sccsa=sum(abs(ajphiv(1:nnd))*darea(1:nnd)) if (sccsa > zero) then rhotjava=sum(rhot_tab(1:nnd)*abs(ajphiv(1:nnd)) . *darea(1:nnd))/sccsa rhot2java=sum(rhot_tab(1:nnd)*rhot_tab(1:nnd)* . abs(ajphiv(1:nnd))*darea(1:nnd))/sccsa end if end if c factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile drhotpav=sqrt(8._wp_*(rhot2pav-rhotpav**2)) drhotjava=sqrt(8._wp_*(rhot2java-rhotjava**2)) rhoppav=frhopol(rhotpav) rhopjava=frhopol(rhotjava) dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*fdvdrhot(rhoppav)) ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*fdadrhot(rhopjava)) call ratioj(rhopjava,ratjamx,ratjbmx,ratjplmx) call profwidth(nnd,rhot_tab,dpdv,rhotp,dpdvmx,drhotp) if(ieccd.ne.0) then call profwidth(nnd,rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi) else rhotjfi=0.0d0 ajmxfi=0.0d0 ajphip=0.0d0 drhotjfi=0.0d0 ratjamx=1.0d0 ratjbmx=1.0d0 ratjplmx=1.0d0 end if end subroutine postproc_profiles subroutine profwidth(nd,xx,yy,xpk,ypk,dxxe) use const_and_precisions, only : wp_,emn1 use utils, only : locatex, locate, intlin, vmaxmini implicit none c arguments integer :: nd real(wp_), dimension(nd) :: xx,yy real(wp_), intent(out) :: xpk,ypk,dxxe c local variables integer :: imn,imx,ipk,ie real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2 real(wp_) :: ypkp,ypkm c call vmaxmini(yy,nd,ymn,ymx,imn,imx) ypk=0.0_wp_ xmx=xx(imx) xmn=xx(imn) if (abs(ymx).gt.abs(ymn)) then ipk=imx ypkp=ymx xpkp=xmx if(abs(ymn/ymx).lt.1.0e-2_wp_) ymn=0.0_wp_ ypkm=ymn xpkm=xmn else ipk=imn ypkp=ymn xpkp=xmn if(abs(ymx/ymn).lt.1.0e-2_wp_) ymx=0.0_wp_ ypkm=ymx xpkm=xmx end if if(xpkp.gt.0.0_wp_) then xpk=xpkp ypk=ypkp yye=ypk*emn1 call locatex(yy,nd,1,ipk,yye,ie) if(ie.gt.0.and.ie.lt.nd) then call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1) else rte1=0.0_wp_ end if call locatex(yy,nd,ipk,nd,yye,ie) if(ie.gt.0.and.ie.lt.nd) then call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2) else rte2=0.0_wp_ end if else ipk=2 xpk=xx(2) ypk=yy(2) rte1=0.0_wp_ yye=ypk*emn1 call locate(yy,nd,yye,ie) if(ie.gt.0.and.ie.lt.nd) then call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2) else rte2=0.0_wp_ end if end if dxxe=rte2-rte1 if(ymx.ne.0.and.ymn.ne.0) dxxe=-dxxe return end c c c subroutine polarcold(exf,eyif,ezf,elf,etf) use const_and_precisions, only : wp_ implicit none c arguments real(wp_) :: exf,eyif,ezf,elf,etf c local variables real(wp_) :: anpl2,anpr2,an2,yg2,dy2,aa,e3,qq,p c common/external functions/variables real(wp_) :: anpl,anpr,xg,yg,sox c 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.0_wp_ c eyif=0.0_wp_ c ezf=0.0_wp_ c if(xg.le.0) return c anpl2=anpl*anpl anpr2=anpr*anpr an2=anpl2+anpr2 yg2=yg**2 dy2=1.0_wp_-yg2 aa=1.0_wp_-xg-yg2 e3=1.0_wp_-xg c if(xg.gt.0.0_wp_) then if (anpl.ne.0.0_wp_) then qq=xg*yg/(an2*dy2-aa) p=(anpr2-e3)/(anpl*anpr) ezf=1.0_wp_/sqrt(1.0_wp_+p*p*(1.0_wp_+qq*qq)) exf=p*ezf eyif=qq*exf else if(sox.lt.0.0_wp_) then ezf=1 exf=0 eyif=0 else ezf=0 qq=-aa/(xg*yg) exf=1.0_wp_/sqrt(1.0_wp_+qq*qq) eyif=qq*exf end if end if elf=(anpl*ezf+anpr*exf)/sqrt(an2) etf=sqrt(1.0_wp_-elf*elf) else if(sox.lt.0.0_wp_) then ezf=1 exf=0 eyif=0 else ezf=0 exf=0.0_wp_ eyif=1.0_wp_ end if elf=0 etf=1 end if c return end c c c subroutine pol_limit(ext,eyt) use const_and_precisions, only : wp_,ui=>im,pi implicit none c arguments complex(wp_) :: ext,eyt c local variables real(wp_) :: anx,any,anz,xe2om,ye2om,xe2xm,ye2xm,an2,an,anxy, . sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2,deno,denx,anpl2,dnl, . del0,gam complex(wp_) :: exom,eyom,exxm,eyxm c common/external functions/variables real(wp_) :: anpl,anpr,yg,sox real(wp_), dimension(3) :: bv,anv c common/anv/anv common/nplr/anpl,anpr common/ygyg/yg common/bb/bv common/mode/sox c anx=anv(1) any=anv(2) anz=anv(3) anpl2=anpl*anpl dnl=1.0_wp_-anpl2 del0=sqrt(dnl**2+4.0_wp_*anpl2/yg**2) ffo=0.5_wp_*yg*(dnl+del0) ffx=0.5_wp_*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) c c if anpl=0 the expressions for exxm and eyxm are 0/0 if (anpl.eq.0.0_wp_) then exxm=-ui*sngam eyxm=-ui*csgam else exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx) eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx) end if c if (sox.lt.0.0_wp_) then ext=exom eyt=eyom else ext=exxm eyt=eyxm endif gam=atan(sngam/csgam)*180.0_wp_/pi return end c c c subroutine stokes(ext,eyt,qq,uu,vv) use const_and_precisions, only : wp_ implicit none c arguments complex(wp_) :: ext,eyt real(wp_) :: qq,uu,vv c qq=abs(ext)**2-abs(eyt)**2 uu=2.0_wp_*dble(ext*dconjg(eyt)) vv=2.0_wp_*dimag(ext*dconjg(eyt)) c end subroutine stokes c c c subroutine polellipse(qq,uu,vv,psipol,chipol) use const_and_precisions, only : wp_,pi implicit none c arguments real(wp_) :: qq,uu,vv,psipol,chipol c c real*8 llm,aa,bb,ell c llm=sqrt(qq**2+uu**2) c aa=sqrt((1+llm)/2.0_wp_) c bb=sqrt((1-llm)/2.0_wp_) c ell=bb/aa psipol=0.5_wp_*atan2(uu,qq)*180.0_wp_/pi chipol=0.5_wp_*asin(vv)*180.0_wp_/pi c end subroutine polellipse c c c logical function inside_plasma(rrm,zzm) use const_and_precisions, only : wp_ use graydata_flags, only : iequil use coreprofiles, only : psdbnd use equilibrium, only : zbinf,zbsup,equinum_psi implicit none c arguments real(wp_) :: rrm,zzm c local variables real(wp_) :: psinv c if(iequil.eq.1) then call equian(rrm,zzm,psinv) else call equinum_psi(rrm,zzm,psinv) end if c if (psinv.ge.0.0_wp_.and.psinv.lt.psdbnd) then if (psinv.lt.1.0_wp_.and.zzm.lt.zbinf.or.zzm.gt.zbsup) then inside_plasma=.false. else inside_plasma=.true. end if else inside_plasma=.false. end if c end function inside_plasma c c c subroutine vacuum_rt(xvstart,anv,xvend,ivac) use const_and_precisions, only : wp_ use reflections, only : inters_linewall,inside,rlim,zlim,nlim use graydata_flags, only : dst use beamdata, only : nstep implicit none c arguments integer :: ivac real(wp_), dimension(3) :: xvstart,anv,xvend c local variables integer :: i real(wp_) :: st,rrm,zzm,smax real(wp_), dimension(3) :: walln,xv0,anv0 logical :: plfound c common/external functions/variables real(wp_) :: dstvac logical :: inside_plasma c external inside_plasma c common/dstvac/dstvac c c ivac=1 plasma hit before wall reflection c ivac=2 wall hit before plasma c ivac=-1 vessel (and thus plasma) never crossed c the real argument associated to xvstart is in a common block c used by fwork!!! xv0=xvstart anv0=anv call inters_linewall(xv0/1.0e2_wp_,anv,rlim(1:nlim),zlim(1:nlim), . nlim,smax,walln) smax=smax*1.0e2_wp_ rrm=1.0e-2_wp_*sqrt(xv0(1)**2+xv0(2)**2) zzm=1.0e-2_wp_*xv0(3) if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then ! first wall interface is outside-inside if (dot_product(walln,walln)