diff --git a/Makefile b/Makefile index ca70079..3d5f572 100644 --- a/Makefile +++ b/Makefile @@ -2,11 +2,11 @@ EXE=gray # Objects list -MAINOBJ=gray.o -OTHOBJ=conical.o const_and_precisions.o coreprofiles.o dierckx.o dispersion.o \ - eccd.o eierf.o graydata_anequil.o graydata_flags.o graydata_par.o \ - equilibrium.o magsurf_data.o math.o minpack.o numint.o quadpack.o \ - reflections.o simplespline.o utils.o beamdata.o +MAINOBJ=main.o +OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \ + dierckx.o dispersion.o eccd.o eierf.o graycore.o gray-externals.o \ + gray_params.o equilibrium.o magsurf_data.o math.o minpack.o numint.o \ + quadpack.o reflections.o simplespline.o utils.o # Alternative search paths vpath %.f90 src @@ -14,7 +14,8 @@ vpath %.f src # Fortran compiler name and flags FC=gfortran -FFLAGS=-O3 #-Wall -g -fcheck=all +FFLAGS=-O3 +#FFLAGS=-Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check DIRECTIVES = -DREVISION="'$(shell svnversion src)'" @@ -25,20 +26,24 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) $(FC) $(FFLAGS) -o $@ $^ # Dependencies on modules -gray.o: const_and_precisions.o coreprofiles.o dierckx.o dispersion.o eccd.o \ - graydata_anequil.o graydata_flags.o graydata_par.o \ +main.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \ + graycore.o gray_params.o reflections.o +graycore.o: const_and_precisions.o beamdata.o beams.o equilibrium.o \ + gray-externals.o gray_params.o reflections.o +gray-externals.o: const_and_precisions.o beams.o coreprofiles.o dierckx.o \ + dispersion.o eccd.o gray_params.o \ equilibrium.o magsurf_data.o math.o numint.o quadpack.o \ - reflections.o simplespline.o utils.o beamdata.o + reflections.o simplespline.o utils.o beamdata.o +beams.o: const_and_precisions.o simplespline.o utils.o +beamdata.o: const_and_precisions.o gray_params.o conical.o: const_and_precisions.o -coreprofiles.o: const_and_precisions.o dierckx.o graydata_anequil.o \ - graydata_flags.o simplespline.o utils.o +coreprofiles.o: const_and_precisions.o dierckx.o gray_params.o simplespline.o \ + utils.o dierckx.o: const_and_precisions.o dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o eccd.o: const_and_precisions.o conical.o magsurf_data.o dierckx.o numint.o eierf.o: const_and_precisions.o -graydata_anequil.o: const_and_precisions.o -graydata_flags.o: const_and_precisions.o -graydata_par.o: const_and_precisions.o +gray_params.o: const_and_precisions.o utils.o equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o utils.o magsurf_data.o: const_and_precisions.o math.o: const_and_precisions.o @@ -48,7 +53,6 @@ quadpack.o: const_and_precisions.o reflections.o: const_and_precisions.o utils.o simplespline.o: const_and_precisions.o utils.o: const_and_precisions.o -beamdata.o: const_and_precisions.o # General object compilation command %.o: %.f90 @@ -57,7 +61,7 @@ beamdata.o: const_and_precisions.o %.o: %.f $(FC) $(FFLAGS) -c $< -gray.o:gray.f +gray-externals.o:gray-externals.f $(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< .PHONY: clean install diff --git a/src/beamdata.f90 b/src/beamdata.f90 index 03831b5..83221ca 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -3,13 +3,14 @@ module beamdata implicit none integer, parameter :: jmx=31,kmx=36,nmx=8000 - integer, save :: nrayr,nrayth,nstep + integer, save :: nray,nrayr,nrayth,nstep + real(wp_) :: dst,rwmax real(wp_), dimension(:,:,:), allocatable, save :: psjki,ppabs,ccci,tauv,alphav real(wp_), dimension(:,:,:), allocatable, save :: pdjki,currj,didst integer, dimension(:,:), allocatable, save :: iiv,iop,iow,ihcd,istore real(wp_), dimension(:,:), allocatable, save :: tau1v real(wp_), dimension(:), allocatable, save :: q - real(wp_), dimension(:,:,:), allocatable, save :: yyrfl !(:,:,6) + real(wp_), dimension(:,:,:), allocatable, save :: yyrfl !(6,:,:) real(wp_), dimension(:,:,:), allocatable, save :: ywrk,ypwrk !(6,:,:) real(wp_), dimension(:,:,:), allocatable, save :: xc,xco,du1,du1o !(3,:,:) real(wp_), dimension(:,:,:), allocatable, save :: gri,dgrad2v !(3,:,:) @@ -20,9 +21,80 @@ module beamdata contains - subroutine alloc_beam(ierr) + subroutine init_rtr(rtrparam) + use gray_params, only : rtrparam_type + implicit none + type(rtrparam_type), intent(in) :: rtrparam + + dst=rtrparam%dst + rwmax=rtrparam%rwmax + nrayr=rtrparam%nrayr + nrayth=rtrparam%nrayth + if(nrayr==1) nrayth=1 + nray=(nrayr-1)*nrayth+1 + nstep=rtrparam%nstep + call alloc_beam +! call alloc_beam1 + end subroutine init_rtr + + function rayi2jk(i) result(jk) + implicit none + integer, intent(in) :: i + integer, dimension(2) :: jk + integer :: ioff + + if (i>1) then + ioff = i - 2 + jk(1) = ioff/nrayth ! jr-2 + jk(2) = ioff - (jk(1))*nrayth + 1 ! kt +! jk(2) = mod(ioff,nrayth) + 1 ! kt + jk(1) = jk(1) + 2 ! jr + else + jk = 1 + end if + end function rayi2jk + + function rayi2j(i) result(jr) + implicit none + integer, intent(in) :: i + integer :: jr + +! jr = max(1, (i-2)/nrayth + 2) + if (i>1) then + jr = (i-2)/nrayth + 2 + else + jr = 1 + end if + end function rayi2j + + function rayi2k(i) result(kt) + implicit none + integer, intent(in) :: i + integer :: kt + +! kt = max(1, mod(i-2,nrayth) + 1) + if (i>1) then + kt = mod(i-2,nrayth) + 1 + else + kt = 1 + end if + end function rayi2k + + function rayjk2i(jr,kt) result(i) + implicit none + integer, intent(in) :: jr,kt + integer :: i + +! i = max(1, (jr-2)*nrayth + kt + 1) + if (jr>1) then + i = (jr-2)*nrayth + kt + 1 + else + i = 1 + end if + end function rayjk2i + + subroutine alloc_beam implicit none - integer, intent(out) :: ierr call dealloc_beam allocate(psjki(nrayr,nrayth,nstep), ppabs(nrayr,nrayth,nstep), & @@ -32,18 +104,35 @@ contains iiv(nrayr,nrayth), iop(nrayr,nrayth), & iow(nrayr,nrayth), tau1v(nrayr,nrayth), & ihcd(nrayr,nrayth), istore(nrayr,nrayth), & - q(nrayr), yyrfl(nrayr,nrayth,6), & + q(nrayr), yyrfl(6,nrayr,nrayth), & ywrk(6,nrayr,nrayth), ypwrk(6,nrayr,nrayth), & xc(3,nrayr,nrayth), xco(3,nrayr,nrayth), & du1(3,nrayr,nrayth), du1o(3,nrayr,nrayth), & gri(3,nrayr,nrayth), dgrad2v(3,nrayr,nrayth), & ggri(3,3,nrayr,nrayth), grad2(nrayr,nrayth), & dffiu(nrayr), ddffiu(nrayr), & - ext(nrayr,nrayth,0:4), eyt(nrayr,nrayth,0:4), & - stat=ierr) - if (ierr/=0) call dealloc_beam + ext(nrayr,nrayth,0:4), eyt(nrayr,nrayth,0:4)) end subroutine alloc_beam +! subroutine alloc_beam1 +! implicit none +! integer, intent(out) :: ierr +! +! call dealloc_beam +! allocate(psjki(nray,nstep), ppabs(nray,nstep), pdjki(nray,nstep), & +! ccci(nray,nstep), currj(nray,nstep), didst(nray,nstep), & +! tauv(nray,nstep), alphav(nray,nstep), & +! ywrk(6,nray), ypwrk(6,nray), & +! xc(3,nray), xco(3,nray), & +! du1(3,nray), du1o(3,nray), & +! grad2(nray), dgrad2v(3,nray), & +! gri(3,nray), ggri(3,3,nray), & +! dffiu(nrayr), ddffiu(nrayr), q(nrayr), & +! iiv(nray), iop(nray), iow(nray), & +! ihcd(nray), istore(nray), tau1v(nray), & +! yyrfl(6,nray), ext(nray,0:4), eyt(nray,0:4)) +! end subroutine alloc_beam1 +! subroutine dealloc_beam implicit none if (allocated(psjki)) deallocate(psjki) diff --git a/src/beams.f90 b/src/beams.f90 new file mode 100644 index 0000000..8aef627 --- /dev/null +++ b/src/beams.f90 @@ -0,0 +1,217 @@ +module beams + use const_and_precisions, only : wp_ + implicit none + +contains + + subroutine read_beam0(file_beam,alpha0,beta0,fghz,x00,y00,z00, & + wcsi,weta,rcicsi,rcieta,phiw,phir,unit) + use const_and_precisions, only : pi,vc=>ccgs_ + use utils, only : get_free_unit + implicit none +! arguments + character(len=*), intent(in) :: file_beam + real(wp_), intent(in) :: alpha0,beta0 + real(wp_), intent(out) :: fGHz,x00,y00,z00 + real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw + integer, intent(in), optional :: unit +! local variables + integer :: u + real(wp_) :: ak0,zrcsi,zreta,w0csi,w0eta,d0csi,d0eta + + if (present(unit)) then + u=unit + else + u = get_free_unit() + end if + open(unit=u,file=trim(file_beam),status='OLD',action='READ') + +! fghz wave frequency (GHz) + read(u,*) fGHz +! x00,y00,z00 coordinates of launching point in cm + read(u,*) x00, y00, z00 +! beams parameters in local reference system +! w0 -> waist, d0 -> waist distance from launching point +! phiw angle of spot ellipse + read(u,*) w0csi,w0eta,d0csi,d0eta,phiw + close(u) + + ak0=2.0e9_wp_*pi*fghz/vc + zrcsi=0.5_wp_*ak0*w0csi**2 + zreta=0.5_wp_*ak0*w0eta**2 + + 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 + + end subroutine read_beam0 + + + + subroutine read_beam1(file_beam,alpha0,beta0,fghz,x00,y00,z00, & + wcsi,weta,rcicsi,rcieta,phiw,phir,unit) + use const_and_precisions, only : pi,vc=>ccgs_ + use simplespline, only : spli, difcs + use utils, only : get_free_unit,locate + implicit none +! arguments + character(len=*), intent(in) :: file_beam + real(wp_), intent(in) :: alpha0 + real(wp_), intent(out) :: fghz,x00,y00,z00,beta0 + real(wp_), intent(out) :: wcsi,weta,rcicsi,rcieta,phir,phiw + integer, intent(in), optional :: unit +! local variables + integer :: u,ierr,iopt,ier,nisteer,i,k,ii + real(wp_) :: steer,alphast,betast,x00mm,y00mm,z00mm + real(wp_) :: w1,w2,rci1,rci2,phi1,phi2,dal + real(wp_), dimension(:), allocatable :: alphastv,betastv,x00v,y00v, & + z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v, & + cbeta,cx0,cy0,cz0,cwaist1,cwaist2, & + crci1,crci2,cphi1,cphi2 + + if (present(unit)) then + u=unit + else + u = get_free_unit() + end if + open(unit=u,file=file_beam,status='OLD',action='READ') + + read(u,*) fghz + + read(u,*) nisteer + + allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), & + waist2v(nisteer),rci1v(nisteer),rci2v(nisteer), & + phi1v(nisteer),phi2v(nisteer),x00v(nisteer), & + y00v(nisteer),z00v(nisteer),cbeta(4*nisteer), & + cx0(4*nisteer),cy0(4*nisteer),cz0(4*nisteer), & + cwaist1(4*nisteer),cwaist2(4*nisteer),crci1(4*nisteer), & + crci2(4*nisteer),cphi1(4*nisteer),cphi2(4*nisteer), & + stat=ierr) + if (ierr/=0) then + close(u) + deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v, & + phi1v,phi2v,x00v,y00v,z00v,cbeta, & + cx0,cy0,cz0,cwaist1,cwaist2,crci1,crci2,cphi1,cphi2) + write(*,*) 'cannot allocate arrays for beam data' + stop + end if + + do i=1,nisteer + read(u,*) steer,alphast,betast,x00mm,y00mm,z00mm, & + w1,w2,rci1,rci2,phi1,phi2 +! initial beam data measured in mm -> transformed to cm + 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 + close(u) + + 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) + + if((alpha0 > alphastv(1)).and.(alpha0 < alphastv(nisteer))) then + call locate(alphastv,nisteer,alpha0,k) + dal=alpha0-alphastv(k) + beta0=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 + write(*,*) ' alpha0 outside table range !!!' + if(alpha0 >= alphastv(nisteer)) ii=nisteer + if(alpha0 <= alphastv(1)) ii=1 + beta0=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 + + deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v, & + phi1v,phi2v,x00v,y00v,z00v,cbeta, & + cx0,cy0,cz0,cwaist1,cwaist2,crci1,crci2,cphi1,cphi2) + + end subroutine read_beam1 + + subroutine launchangles2n(alpha,beta,x,y,z,anx,any,anz) + use const_and_precisions, only : degree + implicit none +! arguments + real(wp_), intent(in) :: alpha,beta,x,y,z + real(wp_), intent(out) :: anx,any,anz +! local variables + real(wp_) :: r,anr,anphi + + r=sqrt(x**2+y**2) +! phi=atan2(y,x) + print'(a,2f8.3)','alpha0, beta0 = ',alpha,beta + print'(a,4f8.3)','x00, y00, R00, z00 = ',x,y,r,z + print*,' ' +! +! angles alpha, beta in a local reference system as proposed by Gribov et al +! + anr = -cos(degree*beta)*cos(degree*alpha) + anphi = sin(degree*beta) +! anx = -cos(degree*beta)*cos(degree*alpha) +! any = sin(degree*beta) + + anx = (anr*x - anphi*y)/r + any = (anr*y + anphi*x)/r +! anr = (anx*x + any*y)/r +! anphi = (any*x - anx*y)/r + + anz =-cos(degree*beta)*sin(degree*alpha) + end subroutine launchangles2n + + subroutine xgygcoeff(fghz,ak0,bres,xgcn) + use const_and_precisions, only : qe=>ecgs_,me=>mecgs_,vc=>ccgs_,pi,wce1_ + implicit none +! arguments + real(wp_), intent(in) :: fghz + real(wp_), intent(out) :: ak0,bres,xgcn +! local variables + real(wp_) :: omega + + omega=2.0e9_wp_*pi*fghz ! [rad/s] + ak0=omega/vc ! [rad/cm] +! +! yg=btot/bres +! + bres=omega/wce1_ ! [T] +! +! xg=xgcn*dens19 +! + xgcn=1.0e-5_wp_*qe**2/(pi*me*fghz**2) ! [10^-19 m^3] + end subroutine xgygcoeff +end module beams diff --git a/src/const_and_precisions.f90 b/src/const_and_precisions.f90 index d9d0fdb..0789d2e 100755 --- a/src/const_and_precisions.f90 +++ b/src/const_and_precisions.f90 @@ -102,7 +102,7 @@ REAL(wp_), PARAMETER :: mu0inv = 1._wp_/mu0_ ! ! REAL(wp_), PARAMETER :: mc_ = me_*c_ ! [kg*m/s] ! ! f_ce = fce1_*B (B in Tesla): ! -! REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s] + REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s] ! REAL(wp_), PARAMETER :: fce1_ = wce1_/(2*pi) ! [1/s] ! ! f_pl = fpe1_*sqrt(Ne) (Ne in 1/m**3): ! ! REAL(wp_), PARAMETER :: wpe1_ = 56.4049201 ! [rad/s] diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index 00c271d..95bcfe6 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -11,8 +11,7 @@ module coreprofiles contains subroutine density(psin,dens,ddens) - use graydata_flags, only : iprof -! use graydata_anequil, only : dens0,aln1,aln2 + use gray_params, only : iprof use dierckx, only : splev,splder implicit none ! arguments @@ -68,18 +67,17 @@ contains if(ier > 0) print*,ier if(abs(dens) < 1.0e-10_wp_) dens=zero end if -! if(dens < zero) print*,' DENSITY NEGATIVE',dens - if(dens < zero) then - dens=zero - ddens=zero - end if + if(dens < zero) print*,' DENSITY NEGATIVE',dens +! if(dens < zero) then +! dens=zero +! ddens=zero +! end if end if end subroutine density function temp(psin) use const_and_precisions, only : wp_,zero,one - use graydata_flags, only : iprof -! use graydata_anequil, only : te0,dte0,alt1,alt2 + use gray_params, only : iprof use utils, only : locate use simplespline, only :spli implicit none @@ -105,8 +103,7 @@ contains function fzeff(psin) use const_and_precisions, only : wp_,zero,one - use graydata_flags, only : iprof -! use graydata_anequil, only : zeffan + use gray_params, only : iprof use utils, only : locate use simplespline, only :spli implicit none @@ -144,7 +141,7 @@ contains else u=get_free_unit() end if - open(file=trim(filenm),status='old',unit=u) + open(file=trim(filenm),status='old',action='read',unit=u) read(u,*) n if(allocated(psin)) deallocate(psin) if(allocated(te)) deallocate(te) @@ -179,21 +176,22 @@ contains if(allocated(zeff)) deallocate(zeff) allocate(te(4),ne(3),zeff(1)) - open(file=trim(filenm),status='old',unit=u) + open(file=trim(filenm),status='old',action='read',unit=u) read(u,*) ne(1:3) ! dens0,aln1,aln2 read(u,*) te(1:4) ! te0,dte0,alt1,alt2 read(u,*) zeff(1) ! zeffan close(u) end subroutine read_profiles_an - subroutine tene_scal(te,ne,tfact,nfact,bfact,iscal) + subroutine tene_scal(te,ne,tfact,nfact,bfact,iscal,iprof) implicit none ! arguments real(wp_), dimension(:), intent(inout) :: te,ne real(wp_), intent(in) :: tfact,nfact,bfact - integer, intent(in) :: iscal + integer, intent(in) :: iscal,iprof ! local variables real(wp_) :: aat,aan,ffact + integer :: lastte,lastne if (iscal==0) then aat=2.0_wp_/3.0_wp_ @@ -207,21 +205,28 @@ contains else ffact=bfact end if - te(:)=te(:)*ffact**aat*tfact - ne(:)=ne(:)*ffact**aan*nfact + if (iprof==0) then + lastte=2 + lastne=1 + else + lastte=size(te) + lastne=size(ne) + end if + te(1:lastte)=te(1:lastte)*ffact**aat*tfact + ne(1:lastne)=ne(1:lastne)*ffact**aan*nfact end subroutine tene_scal - subroutine set_prfspl(psin,te,ne,zeff,ssplne) + subroutine set_prfspl(psin,te,ne,zeff,ssplne,psdbndmx) use simplespline, only : difcs use dierckx, only : curfit, splev, splder implicit none ! arguments real(wp_), dimension(:), intent(in) :: psin,te,ne,zeff - real(wp_), intent(in) :: ssplne + real(wp_), intent(in) :: ssplne,psdbndmx ! local variables integer, parameter :: iopt=0, kspl=3 integer :: n, npest, lwrkf, ier - real(wp_) :: xb, xe, fp, xnv, ynv + real(wp_) :: xb, xe, fp, xnv, xxp,xxm,delta2 real(wp_), dimension(:), allocatable :: wf, wrkf integer, dimension(:), allocatable :: iwrkf real(wp_), dimension(1) :: dedge,ddedge,d2dedge @@ -266,25 +271,29 @@ contains ! 2nd order derivative, extending the profile up to psi=psdbnd where ! ne=ne'=ne''=0 ! spline value and derivatives at the edge - call splev(tfn,nsfd,cfn,3,psin(n:n),dedge(1:1),1,ier) - call splder(tfn,nsfd,cfn,3,1,psin(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier) - call splder(tfn,nsfd,cfn,3,2,psin(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier) + call splev(tfn,nsfd,cfn,kspl,psin(n:n),dedge(1:1),1,ier) + call splder(tfn,nsfd,cfn,kspl,1,psin(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier) + call splder(tfn,nsfd,cfn,kspl,2,psin(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier) ! determination of the boundary + psdbnd=psdbndmx + psnpp=psin(n) denpp=dedge(1) ddenpp=ddedge(1) d2denpp=d2dedge(1) - psdbnd=psnpp - if(ddenpp < zero) then - xnv=psnpp-ddenpp/d2denpp - ynv=denpp-0.5_wp_*ddenpp**2/d2denpp - if((d2denpp > zero).and.(ynv >= zero)) then - psdbnd=xnv - else - psdbnd=xnv+sqrt((ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp) - end if - print*,psnpp,denpp,ddenpp,d2denpp,xnv,ynv - print*,psdbnd + + delta2=(ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp + xnv=psnpp-ddenpp/d2denpp + if(delta2 < zero) then + if(xnv > psnpp) psdbnd=min(psdbnd,xnv) + else + xxm=xnv-sqrt(delta2) + xxp=xnv+sqrt(delta2) + if(xxm > psnpp) then + psdbnd=min(psdbnd,xxm) + else if (xxp > psnpp) then + psdbnd=min(psdbnd,xxp) + end if end if deallocate(iwrkf,wrkf,wf) diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index fab68ef..c47a2db 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -2,12 +2,12 @@ module equilibrium use const_and_precisions, only : wp_ implicit none - REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis + REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis,sgnbphi REAL(wp_), SAVE :: btrcen ! used only for Jcd_ASTRA def. REAL(wp_), SAVE :: rcen ! computed as fpol(a)/btrcen REAL(wp_), SAVE :: rmnm,rmxm,zmnm,zmxm REAL(wp_), SAVE :: zbinf,zbsup - REAL(wp_), SAVE :: rup,zup,rlw,zlw +! REAL(wp_), SAVE :: rup,zup,rlw,zlw INTEGER, PARAMETER :: kspl=3,ksplp=kspl+1 ! === 2D spline psi(R,z), normalization and derivatives ========== @@ -61,7 +61,7 @@ contains end if ! open G EQDSK file (see http://fusion.gat.com/efit/g_eqdsk.html) - open(file=trim(filenm),status='old',unit=u) + open(file=trim(filenm),status='old',action='read',unit=u) ! get size of main arrays and allocate them if (id==1) then @@ -165,7 +165,7 @@ contains else u=get_free_unit() end if - open(file=trim(filenm),status='old',unit=u) + open(file=trim(filenm),status='old',action='read',unit=u) read(u,*) rr0m,zr0m,rpam read(u,*) b0 read(u,*) q0,qa,alq @@ -286,7 +286,7 @@ contains real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol real(wp_), dimension(:,:), intent(in) :: psin real(wp_), intent(in) :: psiwbrad - real(wp_), intent(inout) :: sspl,ssfp + real(wp_), intent(in) :: sspl,ssfp real(wp_), intent(in), optional :: r0,rax,zax real(wp_), dimension(:), intent(in), optional :: rbnd,zbnd integer, intent(in), optional :: ixp @@ -295,7 +295,8 @@ contains ! local variables integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup - real(wp_) :: fp,rax0,zax0,psinoptmp,psinxptmp,rbmin,rbmax,rbinf,rbsup,r1,z1 + real(wp_) :: sspln,fp,rax0,zax0,psinoptmp,psinxptmp + real(wp_) :: rbmin,rbmax,rbinf,rbsup,r1,z1 real(wp_), dimension(1) :: fpoli real(wp_), dimension(:), allocatable :: fvpsi,wf,wrk integer, dimension(:), allocatable :: iwrk @@ -332,14 +333,15 @@ contains ! reshape 2D psi array to 1D (transposed) array and compute spline coeffs allocate(fvpsi(nrz)) fvpsi=reshape(transpose(psin),(/nrz/)) + sspln=sspl call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & - kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & + kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) ! if ier=-1 data are re-fitted using sspl=0 if(ier==-1) then - sspl=0.0_wp_ + sspln=0.0_wp_ call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & - kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & + kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) end if deallocate(fvpsi) @@ -375,6 +377,7 @@ contains call splev(tfp,nsf,cfp,kspl,psinr(npsi:npsi),fpoli,1,ier) ! set vacuum value used outside 0<=psin<=1 range fpolas=fpoli(1) + sgnbphi=sign(one,fpolas) ! free temporary arrays deallocate(wrk,iwrk,wf) ! @@ -667,6 +670,7 @@ contains function fq(psin) use const_and_precisions, only : wp_ + use gray_params, only : iequil use simplespline, only :spli use utils, only : locate implicit none @@ -675,12 +679,17 @@ contains real(wp_) :: fq ! local variables integer :: i - real(wp_) :: dps + real(wp_) :: dps,rn - call locate(psinr,nq,psin,i) - i=min(max(1,i),nq-1) - dps=psin-psinr(i) - fq=spli(cq,nq,i,dps) + if (iequil<2) then + rn=frhotor(sqrt(psin)) + fq=q0+(qa-q0)*rn**alq + else + call locate(psinr,nq,psin,i) + i=min(max(1,i),nq-1) + dps=psin-psinr(i) + fq=spli(cq,nq,i,dps) + end if end function fq subroutine set_rhospl(rhop,rhot) @@ -733,6 +742,27 @@ contains frhopol=spli(crhop,nrho,i,dr) end function frhopol + function frhopolv(rhot) + use utils, only : locate + use simplespline, only : spli + implicit none +! arguments + real(wp_), dimension(:), intent(in) :: rhot + real(wp_), dimension(size(rhot)) :: frhopolv +! local variables + integer :: i,i0,j + real(wp_) :: dr + + i0=1 + do j=1,size(rhot) + call locate(rhotr(i0:),nrho-i0+1,rhot(j),i) + i=min(max(1,i),nrho-i0)+i0-1 + dr=rhot(j)-rhotr(i) + frhopolv(j)=spli(crhop,nrho,i,dr) + i0=i + end do + end function frhopolv + function frhotor(rhop) use utils, only : locate use simplespline, only : spli @@ -888,6 +918,7 @@ contains q0=qax qa=q1 alq=qexp + sgnbphi=sign(one,bax) if (present(n)) then nq=n @@ -927,7 +958,6 @@ contains real(wp_), intent(out), optional :: psinv,fpolv,dfpv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz ! local variables - integer :: sgn real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,snt,rhop,rhot ! real(wp_) :: frhopol @@ -948,21 +978,21 @@ contains rhot=rn if(rn <= 1.0_wp_) then rhop=frhopol(rhot) - psinv=rhop*rhop + psinv=rhop**2 else - psinv=1.0_wp_+btaxis*aminor**2/2.0_wp_/psia/qa*(rn*rn-1.0_wp_) + psinv=1.0_wp_+btaxis/(2.0_wp_*psia*qa)*(rpm**2-aminor**2) rhop=sqrt(psinv) end if end if if(rn <= 1.0_wp_) then qq=q0+(qa-q0)*rn**alq - dpsidrp=-btaxis*aminor*rn/qq + dpsidrp=btaxis*aminor*rn/qq dqq=alq*(qa-q0)*rn**(alq-1.0_wp_) - d2psidrp=-btaxis*(1.0_wp_-rn*dqq/qq)/qq + d2psidrp=btaxis/qq*(1.0_wp_-rn*dqq/qq) else - dpsidrp=-btaxis*aminor/qa*rn - d2psidrp=-btaxis/qa + dpsidrp=btaxis*rpm/qa + d2psidrp=btaxis/qa end if if(present(fpolv)) fpolv=btaxis*rmaxis diff --git a/src/gray.f b/src/gray-externals.f similarity index 77% rename from src/gray.f rename to src/gray-externals.f index 8ba7b8d..5f92df0 100644 --- a/src/gray.f +++ b/src/gray-externals.f @@ -1,132 +1,64 @@ - 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 program gray +C use gray_params, only : ipass,igrad +C implicit none +Cc local variables +C real(wp_) :: p0mw1 +Cc common/external functions/variables +C integer :: ierr,index_rt +C real(wp_) :: sox,p0mw,powrfl,taumn,taumx,pabstot,currtot, +Cc +C common/ierr/ierr +C common/mode/sox +C common/p0/p0mw +C common/powrfl/powrfl +C common/index_rt/index_rt +C common/taumnx/taumn,taumx,pabstot,currtot +Cc +C if (ipass.gt.1) then +Cc second pass into plasma +C p0mw1=p0mw +C igrad=0 +Cc +C index_rt=2 +C p0mw=p0mw1*powrfl +C call prfile +C call vectinit2 +C call paraminit +C call ic_rt2 +C call gray_integration +C call after_gray_integration +C pabstott=pabstott+pabstot +C currtott=currtott+currtot +Cc +C index_rt=3 +C sox=-sox +C p0mw=p0mw1*(1.0_wp_-powrfl) +C call prfile +C call vectinit2 +C call paraminit +C call ic_rt2 +C call gray_integration +C call after_gray_integration +C pabstott=pabstott+pabstot +C currtott=currtott+currtot +C end if +Cc +C end program gray 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 + subroutine after_gray_integration(pabs,currt,dpdv,ajphiv) 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 gray_params, only : iwarm,nnd,ipec use beamdata, only : nrayr implicit none +c arguments + real(wp_) :: pabs,currt + real(wp_), dimension(nnd) :: dpdv,ajphiv c local variables integer :: iproj,nfilp,i - real(wp_) :: currtka,pabs,currt - - real(wp_), dimension(nnd) :: dpdv,ajphiv,pins,currins + real(wp_) :: currtka + real(wp_), dimension(nnd) :: 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 @@ -163,15 +95,6 @@ c write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot currtka =currtot*1.0e3_wp_ write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka -c - if (index_rt.eq.1) then - if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq - if(iequil.eq.1) write(*,*) 'ANALTYCAL EQUILIBRIUM' - if(iprof.eq.1) write(*,*) 'PROFILE file : ',filenmprf - if(iprof.eq.0) write(*,*) 'ANALTYCAL PROFILES 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 @@ -201,19 +124,20 @@ c c c c - subroutine after_onestep(i,istop) + subroutine after_onestep(p0,sox,ak0,xgcn,bres,i,istop,ier) 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 gray_params, only : iwarm,istpr0,istpl0,ipass,igrad use equilibrium, only : zbinf,zbsup - use beamdata, only : nrayr,nrayth,psjki,ppabs,ccci,iiv,tauv,ihcd, - . istore,iop,iow,tau1v,yyrfl,ywrk,ext,eyt + use beamdata, only : nrayr,nrayth,dst,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 + real(wp_), intent(in) :: p0,sox,ak0,xgcn,bres + integer, intent(in) :: i + integer, intent(out) :: istop, ier c local variables integer :: iproj,nfilp,irfl,j,k,kkk,iopmin,iowmax,iowmin,ivac,j1, . k1,kkkk @@ -241,7 +165,6 @@ c 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 @@ -254,20 +177,19 @@ c iopmin=100 iowmin=100 iowmax=0 + istop=0 + ier=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) + call gwork(sox,xgcn,bres,j,k) c if(ierr.gt.0) then write(*,*) ' IERR = ', ierr + ier=ierr if(ierr.eq.97) then c igrad=0 c ierr=0 @@ -280,11 +202,10 @@ c exit 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) + call pabs_curr(p0,sox,ak0,i,j,k) iiv(j,k)=i else if(iiv(j,k).gt.1) tauv(j,k,i)=tauv(j,k,i-1) @@ -294,13 +215,13 @@ c exit pabstot=pabstot+ppabs(j,k,iiv(j,k)) currtot=currtot+ccci(j,k,iiv(j,k)) end if - call print_output(i,j,k) + call print_output(p0,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))) + call pol_limit(sox,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 @@ -312,7 +233,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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))) + call pol_limit(sox,ext(j,k,iop(j,k)),eyt(j,k,iop(j,k))) end if if (ipass.gt.1) then @@ -323,7 +244,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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))) + call pol_limit(sox,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) @@ -345,10 +266,10 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ywrk(1:3,j,k)=xv ywrk(4:6,j,k)=anv igrad=0 - call gwork(j,k) + call gwork(sox,xgcn,bres,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))) + call pol_limit(sox,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 @@ -364,13 +285,14 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 + print*,i,j,k + print*,anwcl,xwcl,anvjk(1:2,j,k) 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)) @@ -445,8 +367,8 @@ c and evaluate polarization if (ivac.eq.1) then y(1:3)=xvjk(:,j,k) y(4:6)=anvjk(:,j,k) - call fwork(y,dery) - call pol_limit(exin2,eyin2) + call fwork(sox,xgcn,bres,y,dery) + call pol_limit(sox,exin2,eyin2) call stokes(exin2,eyin2,qqin2,uuin2,vvin2) powloop: do j1=1,nrayr kkkk=nrayth @@ -475,7 +397,7 @@ c initial polarization (possibly reflected) strfl11=i*dst write(6,*) ' ' write(6,*) 'Reflected power fraction =',powrfl - write(66,*) psipol0,chipol0,powrfl + write(66,*) psipol,chipol,powrfl istop=1 end if @@ -484,14 +406,16 @@ c initial polarization (possibly reflected) c c c - subroutine print_output(i,j,k) + subroutine print_output(p0mw,i,j,k) use const_and_precisions, only : wp_,pi - use graydata_flags, only : istpl0 + use gray_params, 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 arguments + real(wp_), intent(in) :: p0mw c local constants integer, parameter :: ndim=6 real(wp_), parameter :: taucr=12.0_wp_ @@ -503,14 +427,13 @@ c local variables . 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, + real(wp_) :: 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 @@ -590,7 +513,7 @@ c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1 . alpha11,tau11,pt11,ajcd11,dids11, . dble(nhf),dble(iohkawa),dble(index_rt) -c call polarcold(exf,eyif,ezf,elf,etf) +c call polarcold(sox,exf,eyif,ezf,elf,etf) end if c central ray only end @@ -618,9 +541,6 @@ c 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' @@ -634,490 +554,27 @@ c . '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 + subroutine bres_anal(bres) 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 + use equilibrium, only : btaxis,aminor,rmaxis,zmaxis 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 +c arguments 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 + rres=btaxis*rmaxis/bres + zmx=zmaxis+aminor + zmn=zmaxis-aminor do i=1,21 zzres=zmn+(i-1)*(zmx-zmn)/2.0e1_wp_ write(70,111) i,bres,rres,zzres @@ -1128,177 +585,6 @@ c print resonant B iequil=1 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_ @@ -1358,17 +644,19 @@ c c c c - subroutine surfq(psinq,qpsi,nq,qval,psival) + subroutine surfq(psinq,qpsi,nq,qval) use const_and_precisions, only : wp_ + use equilibrium, only : rmaxis,zmaxis,zbinf,zbsup,frhotor 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 + real(wp_) :: qval c local variables integer :: ncnt,i1,ipr + real(wp_) :: rup,zup,rlw,zlw,rhot,psival real(wp_), dimension(npoints) :: rcn,zcn c ncnt=(npoints-1)/2 @@ -1376,23 +664,32 @@ 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) + if (i1>0.and.i1 qqv(jp)=qq rn=frhotor(sqrt(pstab(jp))) - qqan=q0+(qa-q0)*rn**alq + qqan=fq(pstab(jp)) - dadr=2.0_wp_*pi*rn*rpam*rpam - dvdr=4.0_wp_*pi*pi*rn*rmaxis*rpam*rpam + dadr=2.0_wp_*pi*rn*aminor**2 + dvdr=4.0_wp_*pi**2*rn*rmaxis*aminor**2 c dadrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dadpsi c dvdrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dvdpsi @@ -2477,91 +1763,11 @@ c spline interpolation of H(lambda,rhop) 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 + use equilibrium, only : frhotor,aminor,rmaxis,zmaxis implicit none c arguments integer :: ipr @@ -2575,8 +1781,8 @@ c 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)) + zcn(ic)=zmaxis+aminor*rn*sin(th*(ic-1)) + rcn(ic)=rmaxis+aminor*rn*cos(th*(ic-1)) c if (ipr.gt.0) write(71,111) ic,h,rcn(ic),zcn(ic) end do @@ -2591,24 +1797,18 @@ 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, + use gray_params, only : iwarm + use beamdata, only : psjki,tauv,alphav,pdjki,ppabs, . currj,didst,ccci,iiv,iop,iow,ihcd,istore,tau1v,nrayr,nrayth, - . nstep,alloc_beam + . nstep implicit none c local variables - integer :: i,j,k,ierr + integer :: i,j,k c common/external functions/variables - real(wp_) :: jclosest + integer :: 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_ @@ -2641,17 +1841,6 @@ c 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_ @@ -2689,18 +1878,20 @@ c use const_and_precisions, only : wp_ implicit none c common/external functions/variables - integer :: istep,istpr,istpl,ierr,istop + integer :: istep,istpr,istpl,ierr,istop,index_rt c common/istep/istep common/istgr/istpr,istpl common/ierr/ierr common/istop/istop + common/index_rt/index_rt c istpr=0 istpl=1 ierr=0 istep=0 istop=0 + index_rt=1 c return end @@ -2955,14 +2146,16 @@ c c c c - subroutine rkint4 + subroutine rkint4(sox,xgcn,bres) 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, + use gray_params, only : igrad + use beamdata, only : dst,nrayr,nrayth,ywrk,ypwrk,grad2,dgrad2v, . gri,ggri implicit none +c arguments + real(wp_), intent(in) :: sox,xgcn,bres c parameter integer, parameter :: ndim=6 c local variables @@ -2998,17 +2191,17 @@ c ddgr(iv,jv)=ggri(iv,jv,j,k) end do end do - call fwork(yy,fk2) + call fwork(sox,xgcn,bres,yy,fk2) c do ieq=1,ndim yy(ieq)=y(ieq)+fk2(ieq)*hh end do - call fwork(yy,fk3) + call fwork(sox,xgcn,bres,yy,fk3) c do ieq=1,ndim yy(ieq)=y(ieq)+fk3(ieq)*h end do - call fwork(yy,fk4) + call fwork(sox,xgcn,bres,yy,fk4) c do ieq=1,ndim ywrk(ieq,j,k)=y(ieq) @@ -3026,14 +2219,15 @@ c c c c - subroutine gwork(j,k) + subroutine gwork(sox,xgcn,bres,j,k) use const_and_precisions, only : wp_ - use graydata_flags, only : igrad + use gray_params, only : igrad use beamdata, only : ywrk,ypwrk,grad2,dgrad2v,gri,ggri implicit none c local constants integer, parameter :: ndim=6 c arguments + real(wp_), intent(in) :: sox,xgcn,bres integer :: j,k c local variables integer :: ieq,iv,jv @@ -3063,7 +2257,7 @@ c end do end if c - call fwork(yy,yyp) + call fwork(sox,xgcn,bres,yy,yyp) c do ieq=1,ndim ypwrk(ieq,j,k)=yyp(ieq) @@ -3076,14 +2270,15 @@ c c c c - subroutine fwork(y,dery) + subroutine fwork(sox,xgcn,bres,y,dery) use const_and_precisions, only : wp_ - use graydata_flags, only : idst,igrad + use gray_params, only : idst,igrad implicit none c local constants integer, parameter :: ndim=6 c arguments real(wp_), dimension(ndim) :: y,dery + real(wp_), intent(in) :: sox,xgcn,bres c local variables integer :: iv,jv real(wp_) :: xx,yy,zz,yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl, @@ -3092,7 +2287,7 @@ c local variables 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, + real(wp_) :: 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 @@ -3100,7 +2295,6 @@ 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 @@ -3131,7 +2325,7 @@ c if (yy.lt.0.0_wp_) then c phi=-phi c end if c - call plas_deriv(xx,yy,zz) + call plas_deriv(xgcn,bres,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) @@ -3198,6 +2392,7 @@ c 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) dbgr(iv)=dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv) end do end do @@ -3269,14 +2464,13 @@ c c c c - subroutine plas_deriv(xx,yy,zz) + subroutine plas_deriv(xgcn,bres,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 + use gray_params, only : iequil + use equilibrium, only : equinum_fpol, equinum_psi, sgnbphi, equian implicit none c arguments - real(wp_) :: xx,yy,zz + real(wp_) :: xgcn,bres,xx,yy,zz c local variables integer :: iv,jv real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy, @@ -3284,12 +2478,11 @@ c local variables 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, + real(wp_) :: 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 @@ -3347,7 +2540,7 @@ c call equinum_fpol(psinv,fpolv,dfpv) end if - call sub_xg_derxg(psinv) + call sub_xg_derxg(xgcn,psinv) yg=fpolv/(rrm*bres) bphi=fpolv/rrm btot=abs(bphi) @@ -3445,75 +2638,6 @@ c 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_ @@ -3584,20 +2708,19 @@ c c c c - subroutine sub_xg_derxg(psinv) + subroutine sub_xg_derxg(xgcn,psinv) use const_and_precisions, only : wp_ use equilibrium, only : psia use coreprofiles, only : density implicit none c arguments - real(wp_) :: psinv + real(wp_) :: xgcn,psinv c common/external functions/variables - real(wp_) :: xg,xgcn,dxgdpsi + real(wp_) :: xg,dxgdpsi real(wp_) :: dens,ddenspsin c common/xgxg/xg common/dxgdps/dxgdpsi - common/xgcn/xgcn common/dens/dens,ddenspsin c xg=0.0_wp_ @@ -3613,17 +2736,22 @@ c c c c - subroutine ic_gb + subroutine ic_gb(x00,y00,z00,anx0c,any0c,anz0c,ak0,xgcn,bres, + . wcsi,weta,rcicsi,rcieta,phiw,phir,sox,psipol0,chipol0) 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, + use gray_params, only : ipol + use beamdata, only : nrayr,nrayth,rwmax,ywrk0=>ywrk,ypwrk0=>ypwrk, . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt implicit none +c arguments + real(wp_), intent(in) :: x00,y00,z00,anx0c,any0c,anz0c + real(wp_), intent(in) :: ak0,xgcn,bres + real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir + real(wp_), intent(in) :: sox,psipol0,chipol0 c local constants integer, parameter :: ndim=6,ndimm=3 c local variables @@ -3640,20 +2768,18 @@ c local variables 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 + real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi, + . ddr11,psinv,dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi, + . psipol,chipol,psinv11 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 + common/polcof/psipol,chipol + common/psinv11/psinv11 c csth=anz0c snth=sqrt(1.0_wp_-csth**2) @@ -3671,8 +2797,8 @@ c c csphir=cos(phirrad) c snphir=sin(phirrad) c - wwcsi=2.0_wp_*akinv/wcsi**2 - wweta=2.0_wp_*akinv/weta**2 + wwcsi=2.0_wp_/(ak0*wcsi**2) + wweta=2.0_wp_/(ak0*weta**2) c if(phir.ne.phiw) then sk=(rcicsi+rcieta) @@ -3700,9 +2826,9 @@ c qi2=rci2-ui*ww2 end if -c w01=sqrt(2.0_wp_*akinv/ww1) +c w01=sqrt(2.0_wp_/(ak0*ww1)) c z01=-rci1/(rci1**2+ww1**2) -c w02=sqrt(2.0_wp_*akinv/ww2) +c w02=sqrt(2.0_wp_/(ak0*ww2)) c z02=-rci2/(rci2**2+ww2**2) qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2 @@ -3739,7 +2865,7 @@ c z02=-rci2/(rci2**2+ww2**2) if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) da=2.0_wp_*pi/dble(nrayth) c - ddfu=2.0_wp_*dr**2*akinv + ddfu=2.0_wp_*dr**2/ak0 do j=1,nrayr u=dble(j-1) c ffi=u**2*ddfu/2.0_wp_ @@ -3857,10 +2983,10 @@ c c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) - call fwork(ytmp,yptmp) + call fwork(sox,xgcn,bres,ytmp,yptmp) if(ipol.eq.0) then - call pol_limit(ext(j,k,0),eyt(j,k,0)) + call pol_limit(sox,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))) @@ -3878,6 +3004,8 @@ c eyt(j,k,0)= 0.0_wp_ end if endif + psipol=psipol0 + chipol=chipol0 c grad2(j,k)=gr2 dgrad2v(1,j,k)=dgr2x @@ -3907,6 +3035,7 @@ c . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one + psinv11=psinv ddr110=dd end if if(j.eq.nrayr.and.k.eq.1) then @@ -3930,16 +3059,21 @@ c c c c - subroutine ic_rt + subroutine ic_rt(x00,y00,z00,anx0c,any0c,anz0c,ak0,xgcn,bres, + . wcsi,weta,rcicsi,rcieta,phiw,phir,sox,psipol0,chipol0) 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, + use gray_params, only : ipol + use beamdata, only : nrayr,nrayth,rwmax,ywrk0=>ywrk,ypwrk0=>ypwrk, . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt implicit none +c arguments + real(wp_), intent(in) :: x00,y00,z00,anx0c,any0c,anz0c + real(wp_), intent(in) :: ak0,xgcn,bres + real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir + real(wp_), intent(in) :: sox,psipol0,chipol0 c local constants integer, parameter :: ndim=6,ndimm=3 c local variables @@ -3950,19 +3084,18 @@ c local variables . 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 + real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv,dens,ddens, + . tekev,anpl,anpr,brr,bphi,bzz,ajphi,psipol,chipol,psinv11 + 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 + common/polcof/psipol,chipol + common/psinv11/psinv11 c csth=anz0c snth=sqrt(1.0_wp_-csth**2) @@ -4048,10 +3181,10 @@ c c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) - call fwork(ytmp,yptmp) + call fwork(sox,xgcn,bres,ytmp,yptmp) if(ipol.eq.0) then - call pol_limit(ext(j,k,0),eyt(j,k,0)) + call pol_limit(sox,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))) @@ -4075,6 +3208,8 @@ c deltapol=phix-phiy, phix =0 end if end if endif + psipol=psipol0 + chipol=chipol0 c do iv=1,3 gri(iv,j,k)=0.0_wp_ @@ -4100,6 +3235,7 @@ c . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then + psinv11=psinv 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, @@ -4124,13 +3260,14 @@ c c c c - subroutine ic_rt2 + subroutine ic_rt2(sox,xgcn,bres) 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 arguments + real(wp_), intent(in) :: sox,xgcn,bres c local constants integer, parameter :: ndim=6,ndimm=3 c local variables @@ -4140,7 +3277,8 @@ c local variables 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 + . dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi,psipol,chipol, + . psinv11 c common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nplr/anpl,anpr @@ -4148,6 +3286,8 @@ c common/parpl/brr,bphi,bzz,ajphi common/dens/dens,ddens common/tete/tekev + common/polcof/psipol,chipol + common/psinv11/psinv11 c do j=1,nrayr do k=1,nrayth @@ -4180,13 +3320,13 @@ c c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) - call fwork(ytmp,yptmp) + call fwork(sox,xgcn,bres,ytmp,yptmp) - call pol_limit(ext(j,k,0),eyt(j,k,0)) + call pol_limit(sox,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) + call polellipse(qq,uu,vv,psipol,chipol) c do iv=1,3 gri(iv,j,k)=0.0_wp_ @@ -4212,6 +3352,7 @@ c . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then + psinv11=psinv 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, @@ -4240,8 +3381,7 @@ c 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 + use beamdata, only : nrayr,nrayth,rwmax,q implicit none c local variables integer :: j @@ -4319,19 +3459,19 @@ c c c c - subroutine pabs_curr(i,j,k) + subroutine pabs_curr(p0mw,sox,ak0,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 gray_params, only : iwarm,ilarm,ieccd,imx + use beamdata, only : dst,psjki,tauv,alphav,pdjki,ppabs,currj, + . didst,ccci,q,tau1v use coreprofiles, only : temp, fzeff use dispersion, only : harmnumber, warmdisp + use equilibrium, only : rmaxis,sgnbphi use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl implicit none c arguments integer, intent(in) :: i,j,k + real(wp_), intent(in) :: p0mw, sox, ak0 c local constants real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ c local variables @@ -4343,13 +3483,12 @@ c local variables real(wp_), dimension(:), allocatable :: eccdpar c common/external functions/variables integer :: nhmin,nhmax,iokhawa,ierr - real(wp_) :: p0mw,ak0,akinv,fhz,dersdst + real(wp_) :: dersdst real(wp_) :: psinv,xg,yg,tekev,dens,ddens,btot - real(wp_) :: anpl,anpr,vgm,derdnm,sox,anprre,anprim,alpha,effjcd, + real(wp_) :: anpl,anpr,vgm,derdnm,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 @@ -4361,8 +3500,6 @@ c 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 @@ -4370,7 +3507,7 @@ c c dvvl=1.0_wp_ rbavi=1.0_wp_ - rrii=rr0m + rrii=rmaxis rhop=sqrt(psinv) c c @@ -4586,7 +3723,7 @@ c 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 gray_params, only : nnd use equilibrium, only : frhotor,frhopol implicit none integer :: it,ipec @@ -4655,7 +3792,7 @@ c psi grid at mid points, dimension nnd+1, for use in pec_tab 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 gray_params, only : nnd use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv implicit none c local constants @@ -4734,7 +3871,7 @@ c now ajphiv is J_phi=dI/dA [MA/m^2] 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 gray_params, only : nnd use beamdata, only : nstep use utils, only : locatex,intlin c arguments @@ -4827,7 +3964,7 @@ c local variables . 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 gray_params, only : nnd,ieccd use equilibrium, only : frhopol implicit none real(wp_), intent(in) :: pabs,currt @@ -4965,20 +4102,19 @@ c c c c - subroutine polarcold(exf,eyif,ezf,elf,etf) + subroutine polarcold(sox,exf,eyif,ezf,elf,etf) use const_and_precisions, only : wp_ implicit none c arguments - real(wp_) :: exf,eyif,ezf,elf,etf + real(wp_) :: sox,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 + real(wp_) :: anpl,anpr,xg,yg c common/nplr/anpl,anpr common/xgxg/xg common/ygyg/yg - common/mode/sox c c dcold dispersion c dielectric tensor (transposed) @@ -5036,10 +4172,11 @@ c c c c - subroutine pol_limit(ext,eyt) + subroutine pol_limit(sox,ext,eyt) use const_and_precisions, only : wp_,ui=>im,pi implicit none c arguments + real(wp_), intent(in) :: sox complex(wp_) :: ext,eyt c local variables real(wp_) :: anx,any,anz,xe2om,ye2om,xe2xm,ye2xm,an2,an,anxy, @@ -5047,14 +4184,13 @@ c local variables . del0,gam complex(wp_) :: exom,eyom,exxm,eyxm c common/external functions/variables - real(wp_) :: anpl,anpr,yg,sox + real(wp_) :: anpl,anpr,yg 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) @@ -5142,9 +4278,9 @@ c c logical function inside_plasma(rrm,zzm) use const_and_precisions, only : wp_ - use graydata_flags, only : iequil + use gray_params, only : iequil use coreprofiles, only : psdbnd - use equilibrium, only : zbinf,zbsup,equinum_psi + use equilibrium, only : zbinf,zbsup,equinum_psi,equian implicit none c arguments real(wp_) :: rrm,zzm @@ -5174,7 +4310,7 @@ 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 : dst implicit none c arguments integer :: ivac diff --git a/src/gray_params.f90 b/src/gray_params.f90 new file mode 100644 index 0000000..05e1cf5 --- /dev/null +++ b/src/gray_params.f90 @@ -0,0 +1,208 @@ +module gray_params + use const_and_precisions, only : wp_ + implicit none + integer, parameter :: lenfnm=256 + + type antctrl_type + real(wp_) :: alpha, beta, power + real(wp_) :: psi, chi + integer :: iox + integer :: ibeam + character(len=lenfnm) :: filenm + end type antctrl_type + + type eqparam_type + real(wp_) :: ssplps, ssplf, factb + integer :: sgnb, sgni, ixp + integer :: iequil, icocos, ipsinorm, idesc, ifreefmt + character(len=lenfnm) :: filenm + end type eqparam_type + + type prfparam_type + real(wp_) :: psnbnd, sspld, factne, factte + integer :: iscal, irho !, icrho, icte, icne, iczf + integer :: iprof + character(len=lenfnm) :: filenm + end type prfparam_type + + type rtrparam_type + real(wp_) :: rwmax, dst + integer :: nrayr, nrayth, nstep + integer :: igrad, idst, ipass, ipol + end type rtrparam_type + + type hcdparam_type + integer :: iwarm, ilarm, imx, ieccd + end type hcdparam_type + + type outparam_type + integer :: ipec, nrho, istpr, istpl + end type outparam_type + + + integer, save :: iequil,iprof,ipol + integer, save :: iwarm,ilarm,imx,ieccd + integer, save :: igrad,idst,ipass + integer, save :: istpr0,istpl0 + integer, save :: ipec,nnd + +contains + + subroutine read_inputs(filenm,antctrl,eqparam,rwall,prfparam,outparam,unit) + use const_and_precisions, only : wp_ + use utils, only : get_free_unit + implicit none +! arguments + character(len=*), intent(in) :: filenm + type(antctrl_type), intent(out) :: antctrl + type(eqparam_type), intent(out) :: eqparam + real(wp_), intent(out) :: rwall + type(prfparam_type), intent(out) :: prfparam + type(outparam_type), intent(out) :: outparam + integer, intent(in), optional :: unit +! local variables + integer :: u + + if (present(unit)) then + u=unit + else + u = get_free_unit() + end if + open(u,file=filenm,status= 'old',action='read') + +! alpha0, beta0 (cartesian) launching angles + read(u,*) antctrl%alpha, antctrl%beta +! p0mw injected power (MW) + read(u,*) antctrl%power +! abs(iox)=1/2 OM/XM +! psipol0,chipol0 polarization angles at the antenna (if iox<0) + read(u,*) antctrl%iox, antctrl%psi, antctrl%chi + +! ibeam=0 :read data for beam as above +! ibeam=1 :read data from file simple astigmatic beam +! ibeam=2 :read data from file general astigmatic beam + read(u,*) antctrl%ibeam + read(u,*) antctrl%filenm + +! iequil=0 :vacuum +! iequil=1 :analytical equilibrium +! iequil=2 :read eqdsk + read(u,*) eqparam%iequil + read(u,*) eqparam%filenm +! icocos :index for equilibrium from COCOS - O. Sauter Feb 2012 +! ipsinorm :0 standard EQDSK format, 1 format Portone summer 2004 + read(u,*) eqparam%icocos, eqparam%ipsinorm, eqparam%idesc, eqparam%ifreefmt +! ixp=0,-1,+1 : no X point , bottom/up X point +! ssplps : spline parameter for psi interpolation + read(u,*) eqparam%ixp, eqparam%ssplps !, eqparam%ssplf + eqparam%ssplf=0.01_wp_ +! signum of toroidal B and I +! factb factor for magnetic field (only for numerical equil) +! scaling adopted: beta=const, qpsi=const, nustar=const + read(u,*) eqparam%sgnb, eqparam%sgni, eqparam%factb + + read(u,*) rwall + +! iprof=0 :analytical density and temp. profiles +! iprof>0 :numerical density and temp. profiles + read(u,*) prfparam%iprof, prfparam%irho ! irho=0,1,2 -> num profiles vs rhot,rhop,psin + read(u,*) prfparam%filenm +! psbnd value of psi ( > 1 ) of density boundary + read(u,*) prfparam%psnbnd !, prfparam%sspld + prfparam%sspld=0.001_wp_ +! iscal ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling +! factT factn factor for Te&ne scaling + read(u,*) prfparam%factte, prfparam%factne, prfparam%iscal + +! ipec=0/1 :pec profiles grid in psi/rhop +! nrho :number of grid steps for pec profiles +1 + read(u,*) outparam%ipec, outparam%nrho +! istpr0 projection step = dsdt*istprj +! istpl0 plot step = dsdt*istpl + read(u,*) outparam%istpr, outparam%istpl + close(u) + end subroutine read_inputs + + subroutine read_params(filenm,rtrparam,hcdparam,unit) + use utils, only : get_free_unit + implicit none +! arguments + character(len=*), intent(in) :: filenm + type(rtrparam_type), intent(out) :: rtrparam + type(hcdparam_type), intent(out) :: hcdparam + integer, intent(in), optional :: unit +! local variables + integer :: u + + if (present(unit)) then + u=unit + else + u = get_free_unit() + end if + open(u,file=filenm,status= 'old',action='read') + +! nrayr number of rays in radial direction +! nrayth number of rays in angular direction +! rwmax normalized maximum radius of beam power +! rwmax=1 -> last ray at radius = waist + read(u,*) rtrparam%nrayr, rtrparam%nrayth, rtrparam%rwmax +! igrad=0 optical ray-tracing, initial conditions as for beam +! igrad=1 quasi-optical ray-tracing +! igrad=-1 ray-tracing, init. condit. +! from center of mirror and with angular spread +! ipass=1/2 1 or 2 passes into plasma +! ipol=0 compute mode polarization at antenna, ipol=1 use polariz angles + read(u,*) rtrparam%igrad, rtrparam%ipass, rtrparam%ipol +! dst integration step +! nstep maximum number of integration steps +! idst=0/1/2 0 integration in s, 1 integr. in ct, 2 integr. in Sr + read(u,*) rtrparam%dst, rtrparam%nstep, rtrparam%idst + +! iwarm=0 :no absorption and cd +! iwarm=1 :weakly relativistic absorption +! iwarm=2 :relativistic absorption, n<1 asymptotic expansion +! iwarm=3 :relativistic absorption, numerical integration +! ilarm :order of larmor expansion +! imx :max n of iterations in dispersion, imx<0 uses 1st +! iteration in case of failure after |imx| iterations + read(u,*) hcdparam%iwarm,hcdparam%ilarm,hcdparam%imx + +! ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models + read(u,*) hcdparam%ieccd + close(u) + end subroutine read_params + + subroutine set_codepar(eqparam,prfparam,outparam,rtrparam,hcdparam) + implicit none + type(eqparam_type), intent(in) :: eqparam + type(prfparam_type), intent(in) :: prfparam + type(outparam_type), intent(in) :: outparam + type(rtrparam_type), intent(in) :: rtrparam + type(hcdparam_type), intent(in) :: hcdparam + + iequil=eqparam%iequil + iprof=prfparam%iprof + + ipec=outparam%ipec + nnd=outparam%nrho + istpr0=outparam%istpr + istpl0=outparam%istpl + + ipol=rtrparam%ipol + igrad=rtrparam%igrad + idst=rtrparam%idst + ipass=rtrparam%ipass + if (rtrparam%nrayr<5) then + igrad=0 + print*,' nrayr < 5 ! => OPTICAL CASE ONLY' + print*,' ' + end if + + iwarm=hcdparam%iwarm + ilarm=hcdparam%ilarm + imx=hcdparam%imx + ieccd=hcdparam%ieccd + + end subroutine set_codepar + +end module gray_params diff --git a/src/graycore.f90 b/src/graycore.f90 new file mode 100644 index 0000000..73d2220 --- /dev/null +++ b/src/graycore.f90 @@ -0,0 +1,131 @@ +module graycore + use const_and_precisions, only : wp_ + implicit none + +contains + +subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & + psrad,terad,derad,zfc,prfp, rlim,zlim, & + p0,fghz,alpha0,beta0,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir,iox0, & + psipol0,chipol0, dpdv,jcd,pec,icd, outp,rtrp,hcdp,ierr) + use coreprofiles, only : set_prfan, set_prfspl + use gray_params, only : eqparam_type, prfparam_type, outparam_type, & + rtrparam_type, hcdparam_type, set_codepar, igrad, iequil, iprof + use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff + use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl + use beamdata, only : init_rtr, dealloc_beam, nstep, dst + use reflections, only : set_lim + implicit none + type(eqparam_type), intent(in) :: eqp + type(prfparam_type), intent(in) :: prfp + type(outparam_type), intent(in) :: outp + type(rtrparam_type), intent(in) :: rtrp + type(hcdparam_type), intent(in) :: hcdp + real(wp_), dimension(:), allocatable, intent(in) :: psrad, terad, derad, zfc + real(wp_), dimension(:), allocatable, intent(in) :: rv, zv, psinr, fpol, qpsi + real(wp_), dimension(:), allocatable, intent(in) :: rbnd, zbnd, rlim, zlim + real(wp_), dimension(:,:), allocatable, intent(in) :: psin + real(wp_), intent(in) :: psia, rvac, rax, zax + integer, intent(in) :: iox0 + real(wp_), intent(in) :: p0, fghz, psipol0, chipol0 + real(wp_), intent(in) :: alpha0,beta0, x0,y0,z0, w1,w2, ri1,ri2, phiw,phir + real(wp_), intent(out) :: pec,icd + real(wp_), dimension(:), intent(out) :: dpdv,jcd + integer, intent(out) :: ierr + real(wp_), dimension(:), allocatable :: rhotn + + real(wp_) :: sox,ak0,bres,xgcn,anx0,any0,anz0 + integer :: i,iox,istop,index_rt + + integer :: istep + real(wp_) :: st + common/ss/st + common/istep/istep + +! ======= set environment BEGIN ====== + call set_codepar(eqp,prfp,outp,rtrp,hcdp) + + call set_lim(rlim,zlim) + + if(iequil<2) then + call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3)) + call flux_average_an + else + call set_eqspl(rv,zv,psin, psia, psinr,fpol, eqp%ssplps,eqp%ssplf, rvac, & + rax,zax, rbnd,zbnd, eqp%ixp) +! compute rho_pol/rho_tor mapping + allocate(rhotn(size(qpsi))) + call setqphi_num(psinr,abs(qpsi),abs(psia),rhotn) + call set_rhospl(sqrt(psinr),rhotn) + deallocate(rhotn) +! compute flux surface averaged quantities + call flux_average ! requires frhotor for dadrhot,dvdrhot +! print psi surface for q=1.5 and q=2 + call surfq(psinr,qpsi,size(qpsi),1.5_wp_) + call surfq(psinr,qpsi,size(qpsi),2.0_wp_) + end if + + if(iprof==0) then + call set_prfan(terad,derad,zfc) + else + call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd) + end if + + call xgygcoeff(fghz,ak0,bres,xgcn) + call launchangles2n(alpha0,beta0,x0,y0,z0,anx0,any0,anz0) + call init_rtr(rtrp) + +! ======= set environment END ====== + +! ======= pre-proc prints BEGIN ====== +! print Btot=Bres +! print ne, Te, q, Jphi versus psi, rhop, rhot + if (iequil==1) then + call bres_anal(bres) + call print_prof_an + else + call bfield_res(rv,zv,size(rv),size(zv),bres) + call print_prof + end if +! ======= pre-proc prints END ====== + +! ======= main loop BEGIN ====== + iox=iox0 + sox=-1.0_wp_ + if(iox.eq.2) sox=1.0_wp_ + index_rt=1 + call prfile + call paraminit + call vectinit + if(igrad==0) then + call ic_rt(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, & + w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0) + else + call ic_gb(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, & + w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0) + end if +! if(ierr==0) return +! beam/ray propagation +! st0=0.0_wp_ +! if(index_rt.gt.1) st0=strfl11 + do i=1,nstep + istep=i + st=i*dst !+st0 +! advance one step + call rkint4(sox,xgcn,bres) +! calculations after one step: + call after_onestep(p0,sox,ak0,xgcn,bres,i,istop,ierr) + if(istop.eq.1) exit + end do +! postprocessing + call after_gray_integration(pec,icd,dpdv,jcd) + +! ======= main loop END ====== + +! ======= free memory BEGIN ====== + call dealloc_beam +! call unset... +! ======= free memory END ====== +end subroutine gray + +end module graycore diff --git a/src/graydata_anequil.f90 b/src/graydata_anequil.f90 deleted file mode 100644 index 079b695..0000000 --- a/src/graydata_anequil.f90 +++ /dev/null @@ -1,12 +0,0 @@ -module graydata_anequil - use const_and_precisions, only : wp_ - implicit none - - real(wp_), save :: dens0,aln1,aln2 - real(wp_), save :: te0,dte0,alt1,alt2 - real(wp_), save :: rr0,zr0,rpa - real(wp_), save :: b0,rr0m,zr0m,rpam - real(wp_), save :: q0,qa,alq - real(wp_), save :: zeffan - -end module graydata_anequil diff --git a/src/graydata_flags.f90 b/src/graydata_flags.f90 deleted file mode 100644 index 356f8bd..0000000 --- a/src/graydata_flags.f90 +++ /dev/null @@ -1,14 +0,0 @@ -module graydata_flags - use const_and_precisions, only : wp_ - implicit none - - character*255, save :: filenmeqq,filenmprf,filenmbm - real(wp_), save :: sspl, dst - integer, save :: ibeam,iequil,ixp,iprof - integer, save :: iwarm,ilarm,imx,ieccd,ipec,idst - integer, save :: igrad,ipass - integer, save :: ipsinorm,iscal,icocos - integer, save :: nnd,istpr0,istpl0,ipol - integer, save :: neqdsk,nprof - -end module graydata_flags diff --git a/src/graydata_par.f90 b/src/graydata_par.f90 deleted file mode 100644 index 4e8a29c..0000000 --- a/src/graydata_par.f90 +++ /dev/null @@ -1,10 +0,0 @@ -module graydata_par - use const_and_precisions, only : wp_ - implicit none - - real(wp_), save :: rwmax,rwallm - real(wp_), save :: psipol0,chipol0 - real(wp_), save :: factb,factt,factn - real(wp_), save :: sgnbphi,sgniphi - -end module graydata_par diff --git a/src/main.f90 b/src/main.f90 new file mode 100644 index 0000000..ccbab13 --- /dev/null +++ b/src/main.f90 @@ -0,0 +1,126 @@ +program gray_main + use const_and_precisions, only : wp_,one + use graycore, only : gray + use gray_params, only : read_inputs,read_params, antctrl_type,eqparam_type, & + prfparam_type,outparam_type,rtrparam_type,hcdparam_type + use beams, only : read_beam0, read_beam1 + use equilibrium, only : read_equil_an,read_eqdsk,change_cocos,eq_scal, & + set_rhospl,setqphi_num,frhopolv + use coreprofiles, only : read_profiles_an,read_profiles,tene_scal + use reflections, only : range2rect + implicit none + type(antctrl_type) :: antp + type(eqparam_type) :: eqp + type(prfparam_type) :: prfp + type(outparam_type) :: outp + type(rtrparam_type) :: rtrp + type(hcdparam_type) :: hcdp + + real(wp_), dimension(:), allocatable :: psrad, terad, derad, zfc + real(wp_), dimension(:), allocatable :: rv, zv, psinr, fpol, qpsi + real(wp_), dimension(:), allocatable :: rbnd, zbnd, rlim, zlim + real(wp_), dimension(:,:), allocatable :: psin + real(wp_) :: psia, rvac, rax, zax + integer :: iox0 + real(wp_) :: p0mw, fghz, psipol0, chipol0 + real(wp_) :: alpha0, beta0, x0, y0, z0, w1, w2, ri1, ri2, phiw, phir + + real(wp_) :: pec,icd + + integer :: ierr + real(wp_), dimension(:), allocatable :: xrad, rhot, dpdv, jcd + real(wp_) :: rwallm, rmxm, r0m, z0m, dzmx + +! ======= read parameters BEGIN ======= + call read_inputs('graynew.data',antp,eqp,rwallm,prfp,outp) + call read_params('gray_params.data',rtrp,hcdp) +! ======= read parameters END ======= + +! ======= read input data BEGIN ======= +!------------ equilibrium ------------ + if(eqp%iequil<2) then + call read_equil_an(eqp%filenm, rv, zv, fpol, qpsi) +! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0) + psia = sign(one,qpsi(2)*fpol(1)) + else + call read_eqdsk(eqp%filenm, rv,zv,psin, psia, psinr,fpol,qpsi, rvac, & + rax,zax, rbnd,zbnd, rlim,zlim, eqp%ipsinorm,eqp%idesc,eqp%ifreefmt) + call change_cocos(psia, fpol, qpsi, eqp%icocos, 3) + end if +! re-scale B/I and/or force signs. If sgn=0 on input, set to fpol/-psia signs on output + call eq_scal(psia, fpol, eqp%sgni, eqp%sgnb, eqp%factb) + qpsi(1) = sign(qpsi(1),qpsi(1)*qpsi(2)*psia*fpol(1)) + qpsi(2) = sign(qpsi(2),psia*fpol(1)) +!------------- profiles ------------- + if(prfp%iprof==0) then + call read_profiles_an(prfp%filenm, terad, derad, zfc) + else + call read_profiles(prfp%filenm, xrad, terad, derad, zfc) + allocate(psrad(size(xrad))) + if(prfp%irho==0) then + call setqphi_num(psinr,qpsi,psia,rhot) + call set_rhospl(sqrt(psinr),rhot) + psrad=frhopolv(xrad) + else if(prfp%irho == 1) then + psrad=xrad**2 + else + psrad=xrad + end if + deallocate(xrad) + end if +! re-scale input data + call tene_scal(terad,derad,prfp%factte,prfp%factne,eqp%factb,prfp%iscal,prfp%iprof) +!------------- antenna -------------- +! interpolate beam table if antctrl%ibeam>0 + if(antp%ibeam>0) then + call read_beam1(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir) + else + call read_beam0(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir) + end if + alpha0=antp%alpha + beta0=antp%beta + p0mw=antp%power + psipol0=antp%psi + chipol0=antp%chi + iox0=antp%iox +!--------------- wall --------------- +! set simple limiter if not read from EQDSK +! need to clean up... + r0m=sqrt(x0**2+y0**2)*0.01_wp_ + dzmx=rtrp%dst*rtrp%nstep*0.01_wp_ + z0m=z0*0.01_wp_ + if (.not.allocated(rlim).or.rtrp%ipass<0) then + rtrp%ipass=abs(rtrp%ipass) + if(eqp%iequil<2) then + rmxm=(rv(1)+rv(2))*0.01_wp_ + else + rmxm=rv(size(rv)) + end if + call range2rect(rwallm,max(r0m,rmxm),z0m-dzmx,z0m+dzmx,rlim,zlim) + end if +! ======= read input data END ======= + +! ========================= MAIN SUBROUTINE CALL ========================= + allocate(dpdv(outp%nrho),jcd(outp%nrho)) + call gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & + psrad,terad,derad,zfc,prfp, rlim,zlim, & + p0mw,fghz,alpha0,beta0,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir,iox0, & + psipol0,chipol0, dpdv,jcd,pec,icd, outp,rtrp,hcdp,ierr) +! ======================================================================== + +! ======= control prints BEGIN ======= + if(ierr/=0) print*,' IERR = ', ierr + print*,' ' + print*,'Pabs (MW), Icd (kA) = ', pec,icd*1.0e3_wp_ +! ======= control prints END ======= + +! ======= free memory BEGIN ======= + if(allocated(psrad)) deallocate(psrad) + if(allocated(terad)) deallocate(terad, derad, zfc) + if(allocated(rv)) deallocate(rv, zv, fpol, qpsi) + if(allocated(psin)) deallocate(psin, psinr) + if(allocated(rbnd)) deallocate(rbnd,zbnd) + if(allocated(rlim)) deallocate(rlim,zlim) + if(allocated(dpdv)) deallocate(dpdv, jcd) +! ======= free memory END ====== +end program gray_main \ No newline at end of file diff --git a/src/reflections.f90 b/src/reflections.f90 index 42f5fa6..1a3391e 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -4,12 +4,13 @@ module reflections ! === 1D array limiter Rlim_i, Zlim_i integer, public, save :: nlim + real(wp_), public, save :: rwallm real(wp_), public, dimension(:), allocatable, save :: rlim,zlim private public :: reflect,inters_linewall,inside public :: linecone_coord,interssegm_coord,interssegm - public :: alloc_lim,wall_refl + public :: alloc_lim,wall_refl,range2rect,set_lim contains @@ -313,7 +314,30 @@ subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl) extr=dot_product(vv1,evrfl) eytr=dot_product(vv2,evrfl) eztr=dot_product(vv3,evrfl) -end +end subroutine wall_refl + + subroutine range2rect(xmin,xmax,ymin,ymax,xv,yv) + implicit none + real(wp_), intent(in) :: xmin,xmax,ymin,ymax + real(wp_), intent(out), dimension(:), allocatable :: xv,yv + if (allocated(xv)) deallocate(xv) + if (allocated(yv)) deallocate(yv) + allocate(xv(5),yv(5)) + xv=(/xmin,xmax,xmax,xmin,xmin/) + yv=(/ymin,ymin,ymax,ymax,ymin/) + end subroutine range2rect + + subroutine set_lim(rv,zv) + implicit none + real(wp_), intent(in), dimension(:) :: rv,zv + if (allocated(rlim)) deallocate(rlim) + if (allocated(zlim)) deallocate(zlim) + nlim=size(rv) + allocate(rlim(nlim),zlim(nlim)) + rlim=rv + zlim=zv + rwallm=minval(rlim) + end subroutine set_lim end module reflections