diff --git a/Makefile b/Makefile index 4933240..929579a 100644 --- a/Makefile +++ b/Makefile @@ -35,7 +35,7 @@ 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 -beams.o: const_and_precisions.o simplespline.o utils.o +beams.o: const_and_precisions.o dierckx.o reflections.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 gray_params.o simplespline.o \ diff --git a/src/beamdata.f90 b/src/beamdata.f90 index 7e89d3d..c4573fa 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -9,13 +9,13 @@ module beamdata contains subroutine init_rtr(rtrparam,ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) + psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) use gray_params, only : rtrparam_type use const_and_precisions, only : zero,half,two implicit none type(rtrparam_type), intent(in) :: rtrparam real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & - gri,psjki,tauv,alphav,ppabs,didst,ccci + gri,psjki,tauv,alphav,ppabs,dids,ccci real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri real(wp_), dimension(:), intent(out), allocatable :: p0jk complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt @@ -41,7 +41,7 @@ contains nstep=rtrparam%nstep call alloc_beam(ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) + psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) end subroutine init_rtr function rayi2jk(i) result(jk) @@ -101,31 +101,31 @@ contains end function rayjk2i subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) + psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) implicit none real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & - gri,psjki,tauv,alphav,ppabs,didst,ccci + gri,psjki,tauv,alphav,ppabs,dids,ccci real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri real(wp_), dimension(:), intent(out), allocatable :: p0jk complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt integer, dimension(:), intent(out), allocatable :: iiv call dealloc_beam(ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) + psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) allocate(ywork(6,nray), ypwork(6,nray), gri(3,nray), ggri(3,3,nray), & xc(3,nrayth,nrayr), du1(3,nrayth,nrayr), & psjki(nray,nstep), tauv(nray,nstep), alphav(nray,nstep), & - ppabs(nray,nstep), didst(nray,nstep), ccci(nray,nstep), & + ppabs(nray,nstep), dids(nray,nstep), ccci(nray,nstep), & p0jk(nray), ext(nray), eyt(nray), iiv(nray)) end subroutine alloc_beam subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) + psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) implicit none real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & - gri,psjki,tauv,alphav,ppabs,didst,ccci + gri,psjki,tauv,alphav,ppabs,dids,ccci real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri real(wp_), dimension(:), intent(out), allocatable :: p0jk complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt @@ -141,7 +141,7 @@ contains if (allocated(tauv)) deallocate(tauv) if (allocated(alphav)) deallocate(alphav) if (allocated(ppabs)) deallocate(ppabs) - if (allocated(didst)) deallocate(didst) + if (allocated(dids)) deallocate(dids) if (allocated(ccci)) deallocate(ccci) if (allocated(p0jk)) deallocate(p0jk) if (allocated(ext)) deallocate(ext) diff --git a/src/beams.f90 b/src/beams.f90 index 7b71504..0343e1f 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -4,14 +4,13 @@ module beams contains - subroutine read_beam0(file_beam,alpha0,beta0,fghz,x00,y00,z00, & + subroutine read_beam0(file_beam,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 @@ -64,8 +63,7 @@ contains 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_) :: steer,dal real(wp_), dimension(:), allocatable :: alphastv,betastv,x00v,y00v, & z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v, & cbeta,cx0,cy0,cz0,cwaist1,cwaist2, & @@ -93,29 +91,25 @@ contains if (ierr/=0) then close(u) deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v, & - phi1v,phi2v,x00v,y00v,z00v,cbeta, & + 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 + read(u,*) steer,alphastv(i),betastv(i),x00v(i),y00v(i),z00v(i), & + waist1v(i),waist2v(i),rci1v(i),rci2v(i),phi1v(i),phi2v(i) end do close(u) +! initial beam data measured in mm -> transformed to cm + x00v = 0.1_wp_*x00v + y00v = 0.1_wp_*y00v + z00v = 0.1_wp_*z00v + waist1v = 0.1_wp_*waist1v + waist2v = 0.1_wp_*waist2v + rci1v = 10._wp_*rci1v + rci2v = 10._wp_*rci2v iopt=0 call difcs(alphastv,betastv,nisteer,iopt,cbeta,ier) @@ -164,6 +158,556 @@ contains end subroutine read_beam1 + + subroutine read_beam2(file_beam,beamid,alpha0,beta0,fghz,iox,x00,y00,z00, & + wcsi,weta,rcicsi,rcieta,phiw,phir,unit) + use utils, only : get_free_unit, intlin, locate + use reflections, only : inside + use dierckx, only : curfit, splev, surfit, bispev + implicit none + character(len=*), intent(in) :: file_beam + integer, intent(in) :: beamid + real(wp_), intent(inout) :: alpha0,beta0 + real(wp_), intent(out) :: fghz,x00,y00,z00, wcsi,weta,rcicsi,rcieta,phir,phiw + integer, intent(out) :: iox + integer, intent(in), optional :: unit + + character(len=20) :: beamname + integer :: u + integer :: i, ier, nisteer, fdeg, jumprow, nbeam, nalpha, nbeta + integer :: iopt, incheck, nxcoord, nycoord, nxest, nyest, lwrk, kwrk + integer :: nxwaist1, nywaist1, nxwaist2, nywaist2, nxrci1, nyrci1, nxrci2 + integer :: nyrci2, nxphi1, nyphi1, nxphi2, nyphi2, nxx0, nyx0, nxy0, nyy0 + integer :: nxz0, nyz0, kx, ky, ii, npolyg, nmax, lwrk2, in + integer :: nxycoord + integer, DIMENSION(:), ALLOCATABLE :: iwrk + real(wp_) :: alphast,betast, waist1, waist2, rci1, rci2, phi1, phi2 + real(wp_) :: fp, minx, maxx, miny, maxy, eps, xcoord0, ycoord0 + real(wp_), DIMENSION(:), ALLOCATABLE :: x00v, y00v, z00v, alphastv, & + betastv, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, xcoord, & + ycoord, wrk, txwaist1, tywaist1, txwaist2, tywaist2, & + txrci1, tyrci1, txrci2, tyrci2, txphi1, typhi1, txphi2, typhi2, & + txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1, & + cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2, & + xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC, & + ypolygC, xpolygD, ypolygD + real(wp_), DIMENSION(4) :: xvert, yvert + real(wp_), dimension(1) :: fi + integer, parameter :: kspl=1 + real(wp_), parameter :: sspl=0.01_wp_ + + if (present(unit)) then + u=unit + else + u = get_free_unit() + end if + open(unit=u,file=file_beam,status='OLD',action='READ') +!======================================================================================= +! # of beams + read(u,*) nbeam +! +! unused beams' data + jumprow=0 +! c==================================================================================== + do i=1,beamid-1 + read(u,*) beamname, iox, fghz, nalpha, nbeta + jumprow = jumprow+nalpha*nbeta + end do +! c==================================================================================== +! +! beam of interest + read(u,*) beamname, iox, fghz, nalpha, nbeta +! +! c==================================================================================== +! unused beams' data grids + do i=1,(nbeam - beamid) + read(u,*) beamname + end do + do i=1,jumprow + read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2 + end do +! c==================================================================================== +! +! # of elements in beam data grid + nisteer = nalpha*nbeta +! + allocate(alphastv(nisteer),betastv(nisteer),waist1v(nisteer), & + waist2v(nisteer),rci1v(nisteer),rci2v(nisteer),phi1v(nisteer), & + phi2v(nisteer),x00v(nisteer),y00v(nisteer),z00v(nisteer), & + xcoord(nisteer),ycoord(nisteer)) +! +! c==================================================================================== +! beam data grid reading + do i=1,nisteer + read(u,*) alphast,betast,x00,y00,z00,waist1,waist2,rci1,rci2,phi1,phi2 +! +! initial beam data (x00, y00, z00) are measured in mm -> transformed to cm + x00v(i)=0.1d0*x00 + y00v(i)=0.1d0*y00 + z00v(i)=0.1d0*z00 + alphastv(i)=alphast + betastv(i)=betast + waist1v(i)=0.1d0*waist1 + rci1v(i)=1.0d1*rci1 + waist2v(i)=0.1d0*waist2 + rci2v(i)=1.0d1*rci2 + phi1v(i)=phi1 + phi2v(i)=phi2 + end do + close(u) +! c==================================================================================== +! +! fdeg = 0 alpha, beta free variables +! 1 alpha free variable +! 2 beta free variable +! 3 no free variables + fdeg = 2*(1/nalpha) + 1/nbeta + +!####################################################################################### +! +! no free variables + if(fdeg.eq.3) then + alpha0=alphastv(1) + beta0=betastv(1) + x00=x00v(1) + y00=y00v(1) + z00=z00v(1) + wcsi=waist1v(1) + weta=waist2v(1) + rcicsi=rci1v(1) + rcieta=rci2v(1) + phiw=phi1v(1) + phir=phi2v(1) + return + end if +!####################################################################################### +! +! +!####################################################################################### + if(fdeg.eq.2) then +! beta = independent variable +! alpha = dependent variable + xcoord = betastv + ycoord = alphastv + xcoord0 = beta0 + ycoord0 = alpha0 + kx=min(nbeta-1,kspl) +! c==================================================================================== + else +! c==================================================================================== +! alpha = independent variable +! beta = dependent/independent (fdeg = 1/0) + xcoord = alphastv + ycoord = betastv + xcoord0 = alpha0 + ycoord0 = beta0 + nxcoord = nalpha + nycoord = nbeta + kx=min(nalpha-1,kspl) + ky=min(nbeta-1,kspl) + end if +!####################################################################################### +! + iopt = 0 + incheck = 0 +! +!####################################################################################### + if(fdeg.ne.0) then + nxest = kx + nxcoord + 1 + lwrk = (nxcoord*(kx+1)+nxest*(7+3*kx)) + kwrk = nxest + allocate(cycoord(nxest), txycoord(nxest), cwaist1(nxest), & + txwaist1(nxest), cwaist2(nxest), txwaist2(nxest), & + crci1(nxest), txrci1(nxest), crci2(nxest), txrci2(nxest), & + cphi1(nxest), txphi1(nxest), cphi2(nxest), txphi2(nxest), & + cx0(nxest), txx0(nxest), cy0(nxest), txy0(nxest), & + cz0(nxest), txz0(nxest), w(nxcoord), wrk(lwrk), iwrk(kwrk)) +! + w = 1.d0 +! +! 2D interpolation + call curfit(iopt,nxcoord,xcoord,ycoord,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxycoord, & + txycoord,cycoord,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,waist1v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist1, & + txwaist1,cwaist1,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,waist2v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxwaist2, & + txwaist2,cwaist2,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,rci1v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci1, & + txrci1,crci1,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,rci2v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxrci2, & + txrci2,crci2,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,phi1v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi1, & + txphi1,cphi1,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,phi2v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxphi2, & + txphi2,cphi2,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,x00v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxx0, & + txx0,cx0,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,y00v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxy0, & + txy0,cy0,fp,wrk,lwrk,iwrk,ier) +! + call curfit(iopt,nxcoord,xcoord,z00v,w, & + xcoord(1),xcoord(nxcoord),kx,sspl,nxest,nxz0, & + txz0,cz0,fp,wrk,lwrk,iwrk,ier) +! +! check if xcoord0 is out of table range +! incheck = 1 inside / 0 outside + if(xcoord0.gt.xcoord(1).and.xcoord0.lt.xcoord(nisteer)) incheck=1 +! c==================================================================================== + else +! c==================================================================================== + npolyg = 2*(nxcoord+nycoord-2) + minx = minval(xcoord) + maxx = maxval(xcoord) + miny = minval(ycoord) + maxy = maxval(ycoord) + nxest = kx + 1 + int(sqrt(nisteer/2.)) + nyest = ky + 1 + int(sqrt(nisteer/2.)) + nmax = max(nxest,nyest) + eps = 10.**(-8) + lwrk = (nmax-2)**2*(7*nmax-2)+18*nmax+8*nisteer-19 + lwrk2 = (nmax-2)**2*(4*nmax-1)+4*nmax-2 + kwrk = nisteer+(nmax-3)*(nmax-3) + allocate(cwaist1(nxest*nyest), txwaist1(nmax), tywaist1(nmax), & + cwaist2(nxest*nyest), txwaist2(nmax), tywaist2(nmax), & + crci1(nxest*nyest), txrci1(nmax), tyrci1(nmax), & + crci2(nxest*nyest), txrci2(nmax), tyrci2(nmax), & + cphi1(nxest*nyest), txphi1(nmax), typhi1(nmax), & + cphi2(nxest*nyest), txphi2(nmax), typhi2(nmax), & + cx0(nxest*nyest), txx0(nmax), tyx0(nmax), & + cy0(nxest*nyest), txy0(nmax), tyy0(nmax), & + cz0(nxest*nyest), txz0(nmax), tyz0(nmax), & + wrk(lwrk), wrk2(lwrk2), iwrk(kwrk), & + xpolyg(npolyg), ypolyg(npolyg), w(nisteer)) +! + w = 1.d0 +! +! 3D interpolation + call surfit(iopt,nisteer,xcoord,ycoord,waist1v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxwaist1,txwaist1,nywaist1,tywaist1,cwaist1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! + call surfit(iopt,nisteer,xcoord,ycoord,waist2v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxwaist2,txwaist2,nywaist2,tywaist2,cwaist2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! + call surfit(iopt,nisteer,xcoord,ycoord,rci1v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxrci1,txrci1,nyrci1,tyrci1,crci1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! + call surfit(iopt,nisteer,xcoord,ycoord,rci2v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxrci2,txrci2,nyrci2,tyrci2,crci2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! + call surfit(iopt,nisteer,xcoord,ycoord,phi1v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxphi1,txphi1,nyphi1,typhi1,cphi1,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! + call surfit(iopt,nisteer,xcoord,ycoord,phi2v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxphi2,txphi2,nyphi2,typhi2,cphi2,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! + call surfit(iopt,nisteer,xcoord,ycoord,x00v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxx0,txx0,nyx0,tyx0,cx0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! + call surfit(iopt,nisteer,xcoord,ycoord,y00v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxy0,txy0,nyy0,tyy0,cy0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! + call surfit(iopt,nisteer,xcoord,ycoord,z00v,w, & + minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & + nxz0,txz0,nyz0,tyz0,cz0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) +! data range polygon + xpolyg(1:nxcoord) = xcoord(1:nxcoord) + ypolyg(1:nxcoord) = ycoord(1:nxcoord) +! +! c==================================================================================== + do i=1,nycoord-2 + xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord) + xpolyg(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1) + ypolyg(nxcoord+i) = ycoord((i+1)*nxcoord) + ypolyg(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1) + end do +! c==================================================================================== + do i=1,nxcoord + xpolyg(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1) + ypolyg(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1) + end do +! c==================================================================================== +! +! check if (xcoord0, ycoord0) is out of table range +! incheck = 1 inside / 0 outside + if(inside(xpolyg,ypolyg,npolyg,xcoord0,ycoord0)) incheck = 1 + end if + deallocate(wrk,iwrk) +!####################################################################################### +! +! +!####################################################################################### + if(fdeg.ne.0) then +! c==================================================================================== + if(incheck.eq.1) then + call splev(txycoord,nxycoord,cycoord,kx,(/xcoord0/),fi,1,ier) + ycoord0=fi(1) + call splev(txwaist1,nxwaist1,cwaist1,kx,(/xcoord0/),fi,1,ier) + wcsi=fi(1) + call splev(txwaist2,nxwaist2,cwaist2,kx,(/xcoord0/),fi,1,ier) + weta=fi(1) + call splev(txrci1,nxrci1,crci1,kx,(/xcoord0/),fi,1,ier) + rcicsi=fi(1) + call splev(txrci2,nxrci2,crci2,kx,(/xcoord0/),fi,1,ier) + rcieta=fi(1) + call splev(txphi1,nxphi1,cphi1,kx,(/xcoord0/),fi,1,ier) + phiw=fi(1) + call splev(txphi2,nxphi2,cphi2,kx,(/xcoord0/),fi,1,ier) + phir=fi(1) + call splev(txx0,nxx0,cx0,kx,(/xcoord0/),fi,1,ier) + x00=fi(1) + call splev(txy0,nxy0,cy0,kx,(/xcoord0/),fi,1,ier) + y00=fi(1) + call splev(txz0,nxz0,cz0,kx,(/xcoord0/),fi,1,ier) + z00=fi(1) +! c---------------------------------------------------------------------------------- + else +! c---------------------------------------------------------------------------------- + if(xcoord0.ge.xcoord(nisteer)) ii=nisteer + if(xcoord0.le.xcoord(1)) ii=1 +! + xcoord0=xcoord(ii) + ycoord0=ycoord(ii) + x00=x00v(ii) + y00=y00v(ii) + z00=z00v(ii) + wcsi=waist1v(ii) + weta=waist2v(ii) + rcicsi=rci1v(ii) + rcieta=rci2v(ii) + phiw=phi1v(ii) + phir=phi2v(ii) + end if +! c==================================================================================== + else +! c==================================================================================== + if(incheck.eq.0) then + allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord), & + ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord), & + xpolygD(nycoord), ypolygD(nycoord)) +! coordinates of vertices v1,v2,v3,v4 + xvert(1) = xpolyg(1) + xvert(2) = xpolyg(nxcoord) + xvert(3) = xpolyg(nxcoord+nycoord-1) + xvert(4) = xpolyg(2*nxcoord+nycoord-2) + yvert(1) = ypolyg(1) + yvert(2) = ypolyg(nxcoord) + yvert(3) = ypolyg(nxcoord+nycoord-1) + yvert(4) = ypolyg(2*nxcoord+nycoord-2) +! coordinates of side A,B,C,D + xpolygA = xpolyg(1:nxcoord) + ypolygA = ypolyg(1:nxcoord) + xpolygB = xpolyg(nxcoord:nxcoord+nycoord-1) + ypolygB = ypolyg(nxcoord:nxcoord+nycoord-1) + xpolygC = xpolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2) + ypolygC = ypolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2) + xpolygD(1:nycoord-1) = xpolyg(2*nxcoord+nycoord-2:npolyg) + xpolygD(nycoord) = xpolyg(1) + ypolygD(1:nycoord-1) = ypolyg(2*nxcoord+nycoord-2:npolyg) + ypolygD(nycoord) = ypolyg(1) +! c---------------------------------------------------------------------------------- +! search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid +! +! | | +! (6) (5) (4) +! | | +! _ _ _ v4 _________________v3_ _ _ _ +! | C | (1)->(8) outside regions +! | | +! (7) D | | B (3) v1->v4 grid vertices +! | | +! _ _ _ _ |_________________| _ _ _ _ A-D grid sides +! v1 A v2 +! | | +! (8) (1) (2) +! | | +! + if(xcoord0.gt.xvert(1).and.xcoord0.lt.xvert(2).and.ycoord0.le.maxval(ypolygA)) then + in=1 + else if(ycoord0.gt.yvert(2).and.ycoord0.lt.yvert(3).and.xcoord0.ge.minval(xpolygB)) then + in=3 + else if(xcoord0.lt.xvert(3).and.xcoord0.gt.xvert(4).and.ycoord0.ge.minval(ypolygC)) then + in=5 + else if(ycoord0.lt.yvert(4).and.ycoord0.gt.yvert(1).and.xcoord0.le.maxval(xpolygD)) then + in=7 + else if(xcoord0.ge.xvert(2).and.ycoord0.le.yvert(2)) then + in=2 + else if(xcoord0.ge.xvert(3).and.ycoord0.ge.yvert(3)) then + in=4 + else if(xcoord0.le.xvert(4).and.ycoord0.ge.yvert(4)) then + in=6 + else if(xcoord0.le.xvert(1).and.ycoord0.le.yvert(1)) then + in=8 + endif +! c---------------------------------------------------------------------------------- +! +! c---------------------------------------------------------------------------------- +! (xcoord0,ycoord0) is set to its nearest point on (alpha, beta) grid border +! depending on the region +! 1: xcoord0 unchanged, ycoord0 moved on side A +! 3: xcoord0 moved on side B, ycoord0 unchanged +! 5: xcoord0 unchanged, ycoord0 moved on side C +! 7: xcoord0 moved on side D, ycoord0 unchanged +! 2,4,6,8: (xcoord0,ycoord0) set to nearest vertex coordinates +! in 1,3,5,7 incheck is set back to 1 to evaluate x00,y00,z00,waist,rci,phi in +! new (xcoord0,ycoord0) +! in 2,4,6,8 incheck remains 0 and x00,y00,z00,waist,rci,phi values at the +! (xcoord0,ycoord0) vertex are used + alpha0 = xcoord0 + beta0 = ycoord0 + SELECT CASE (in) + CASE (1) + write(*,*) ' beta0 outside table range !!!' +! locate position of xcoord0 with respect to x coordinates of side A + call locate(xpolygA,nxcoord,xcoord0,ii) +! find corresponding y value on side A for xcoord position + call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0) + incheck = 1 + CASE (2) + write(*,*) ' alpha0 and beta0 outside table range !!!' +! xcoord0, ycoord0 set + xcoord0 = xvert(2) + ycoord0 = yvert(2) + ii = nxcoord !indice per assegnare valori waist, rci, phi + CASE (3) + write(*,*) ' alpha0 outside table range !!!' + call locate(ypolygB,nycoord,ycoord0,ii) + call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0) + incheck = 1 + CASE (4) + write(*,*) ' alpha0 and beta0 outside table range !!!' + xcoord0 = xvert(3) + ycoord0 = yvert(3) + ii = nxcoord+nycoord-1 + CASE (5) + write(*,*) ' beta0 outside table range !!!' + call locate(xpolygC,nxcoord,xcoord0,ii) + call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0) + incheck = 1 + CASE (6) + write(*,*) ' alpha0 and beta0 outside table range !!!' + xcoord0 = xvert(4) + ycoord0 = yvert(4) + ii = 2*nxcoord+nycoord-2 + CASE (7) + write(*,*) ' alpha0 outside table range !!!' + call locate(ypolygD,nycoord,ycoord0,ii) + call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0) + incheck = 1 + CASE (8) + write(*,*) ' alpha0 and beta0 outside table range !!!!' + xcoord0 = xvert(1) + ycoord0 = yvert(1) + ii = 1 + END SELECT +! c---------------------------------------------------------------------------------- +! + deallocate(xpolygA, ypolygA, xpolygC, ypolygC, xpolygB, ypolygB, xpolygD, ypolygD) + end if +! c==================================================================================== +! +! c==================================================================================== + if(incheck.eq.1) then + lwrk = 2*(kx+ky+2) + kwrk = 4 + allocate(wrk(lwrk),iwrk(kwrk)) + call bispev(txwaist1,nxwaist1,tywaist1,nywaist1,cwaist1, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + wcsi=fi(1) + call bispev(txwaist2,nxwaist2,tywaist2,nywaist2,cwaist2, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + weta=fi(1) + call bispev(txrci1,nxrci1,tyrci1,nyrci1,crci1, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + rcicsi=fi(1) + call bispev(txrci2,nxrci2,tyrci2,nyrci2,crci2, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + rcieta=fi(1) + call bispev(txphi1,nxphi1,typhi1,nyphi1,cphi1, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + phiw=fi(1) + call bispev(txphi2,nxphi2,typhi2,nyphi2,cphi2, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + phir=fi(1) + call bispev(txx0,nxx0,tyx0,nyx0,cx0, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + x00=fi(1) + call bispev(txy0,nxy0,tyy0,nyy0,cy0, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + y00=fi(1) + call bispev(txz0,nxz0,tyz0,nyz0,cz0, & + kx,ky,(/xcoord0/),1,(/ycoord0/),1,fi,wrk,lwrk,iwrk,kwrk,ier) + z00=fi(1) + deallocate(wrk,iwrk) +! c---------------------------------------------------------------------------------- + else +! c---------------------------------------------------------------------------------- + x00=x00v(ii) + y00=y00v(ii) + z00=z00v(ii) + wcsi=waist1v(ii) + weta=waist2v(ii) + rcicsi=rci1v(ii) + rcieta=rci2v(ii) + phiw=phi1v(ii) + phir=phi2v(ii) + end if +! c==================================================================================== + end if +!####################################################################################### +! + if(fdeg.ne.0) then + deallocate(cycoord, txycoord, cwaist1, txwaist1, cwaist2, & + txwaist2, crci1, txrci1, crci2, txrci2, cphi1, txphi1, & + cphi2, txphi2, cx0, txx0, cy0, txy0, cz0, txz0, w) + else + deallocate(cwaist1, txwaist1, tywaist1, cwaist2, txwaist2, tywaist2, & + crci1, txrci1, tyrci1, crci2, txrci2, tyrci2, & + cphi1, txphi1, typhi1, cphi2, txphi2, typhi2, & + cx0, txx0, tyx0, cy0, txy0, tyy0, cz0, txz0, tyz0, & + wrk2, xpolyg, ypolyg, w) + end if +! +!####################################################################################### +! set correct values for alpha, beta + if(fdeg.eq.2) then + alpha0 = ycoord0 + beta0 = xcoord0 + else + alpha0 = xcoord0 + beta0 = ycoord0 + end if +!####################################################################################### + deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, & + phi2v,x00v,y00v,z00v,xcoord,ycoord) +! + end subroutine read_beam2 + + subroutine launchangles2n(alpha,beta,xv,anv) use const_and_precisions, only : degree implicit none diff --git a/src/dierckx.f90 b/src/dierckx.f90 index 507a1df..6d1f484 100644 --- a/src/dierckx.f90 +++ b/src/dierckx.f90 @@ -80,9 +80,9 @@ contains integer, intent(in) :: nx, ny, kx, ky, mx, my, lwrk, kwrk integer, intent(out) :: ier integer, intent(inout) :: iwrk(kwrk) - real(wp_), intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)) + real(wp_), intent(in) :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(mx), y(my) real(wp_), intent(out) :: z(mx*my) - real(wp_), intent(inout) :: x(mx), y(my), wrk(lwrk) + real(wp_), intent(inout) :: wrk(lwrk) ! local variables integer :: i, iw, lwest ! .. diff --git a/src/eccd.f90 b/src/eccd.f90 index a8d5e89..ed71c1c 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -275,7 +275,8 @@ contains end if resj=0.0_wp_ - do i=0,iokhawa +! do i=0,iokhawa + do i=0,1 resji=0.0_wp_ xx1=amu*(anpl*uleft(i)+ygn-1.0_wp_) xx2=amu*(anpl*uright(i)+ygn-1.0_wp_) @@ -343,9 +344,10 @@ contains real(wp_) :: upl,fpp real(wp_), dimension(npar) :: extrapar ! local variables - integer :: ithn,nhn,nm,np - real(wp_) :: yg,anpl,amu,anprre,ygn,upr,upr2,gam,ee,thn2,thn2u,bb,cth, & - ajbnm,ajbnp,ajbn + integer :: ithn,nhn !,nm,np + real(wp_) :: yg,anpl,amu,anprre,ygn,upr,upr2,gam,ee,thn2,thn2u,bb,cth !, & +! ajbnm,ajbnp,ajbn + real(wp_), dimension(3) :: ajb complex(wp_) :: ex,ey,ez,emxy,epxy yg=extrapar(1) @@ -379,12 +381,14 @@ contains thn2u=upr2*thn2 else ! Full polarization term - nm=nhn-1 - np=nhn+1 - ajbnm=dbesjn(nm, bb) - ajbnp=dbesjn(np, bb) - ajbn=dbesjn(nhn, bb) - thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0_wp_))**2 +! nm=nhn-1 +! np=nhn+1 +! ajbnm=dbesjn(nm, bb) +! ajbnp=dbesjn(np, bb) +! ajbn=dbesjn(nhn, bb) +! thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0_wp_))**2 + ajb=bessel_jn(nhn-1, nhn+1, bb) + thn2u=(abs(ez*ajb(2)*upl+upr*(ajb(3)*epxy+ajb(1)*emxy)/2.0_wp_))**2 end if end if end if diff --git a/src/gray-externals.f90 b/src/gray-externals.f90 index 4b55dd2..5cc4be0 100644 --- a/src/gray-externals.f90 +++ b/src/gray-externals.f90 @@ -556,6 +556,8 @@ do n=1,5 bbb=bres/dble(n) if (bbb.ge.btmn.and.bbb.le.btmx) then + nconts=size(ncpts) + nctot=size(rrcb) call cniteq(rv,zv,btotal,nr,nz,bbb,nconts,ncpts,nctot,rrcb,zzcb) do inc=1,nctot write(70,'(i6,12(1x,e12.5))') inc,bbb,rrcb(inc),zzcb(inc) @@ -593,40 +595,34 @@ ! (based on an older code) use const_and_precisions, only : wp_ implicit none -! local constants - integer, parameter :: icmx=2002 ! arguments - integer :: nr,nz - real(wp_), dimension(nr) :: rqgrid - real(wp_), dimension(nz) :: zqgrid - real(wp_), dimension(nr,nz) :: matr2dgrid - integer :: ncon,icount - integer, dimension(10) :: npts - real(wp_) :: h - real(wp_), dimension(icmx) :: rcon,zcon + integer, intent(in) :: nr,nz + real(wp_), dimension(nr), intent(in) :: rqgrid + real(wp_), dimension(nz), intent(in) :: zqgrid + real(wp_), dimension(nr,nz), intent(in) :: matr2dgrid + real(wp_), intent(in) :: h + integer, intent(inout) :: ncon, icount + integer, dimension(ncon), intent(out) :: npts + real(wp_), dimension(icount), intent(out) :: rcon,zcon ! local variables - integer :: i,j,k,l,ico,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb + integer :: i,j,k,l,nrqmax,iclast,mpl,ix,jx,mxr,n1,jm,jfor,lda,ldb integer :: jabs,jnb,kx,ikx,itm,inext,in integer, dimension(3,2) :: ja - integer, dimension(1000) :: lx + integer, dimension(icount/2-1) :: lx real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y real(wp_), dimension(nr*nz) :: a - logical :: flag1, flag2 + logical :: flag1 - px=0.5_wp_ + px = 0.5_wp_ - a=reshape(matr2dgrid,(/nr*nz/)) + a = reshape(matr2dgrid,(/nr*nz/)) - do ico=1,icmx - rcon(ico)=0.0_wp_ - zcon(ico)=0.0_wp_ - enddo + rcon = 0.0_wp_ + zcon = 0.0_wp_ - nrqmax=nr - nr=nr - nz=nz - drgrd=rqgrid(2)-rqgrid(1) - dzgrd=zqgrid(2)-zqgrid(1) + nrqmax = nr + drgrd = rqgrid(2) - rqgrid(1) + dzgrd = zqgrid(2) - zqgrid(1) ncon = 0 @@ -634,8 +630,8 @@ iclast = 0 icount = 0 - mpl=0 - ix=0 + mpl = 0 + ix = 0 mxr = nrqmax * (nz - 1) n1 = nr - 1 diff --git a/src/graycore.f90 b/src/graycore.f90 index 1076170..2765e51 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -12,8 +12,8 @@ contains use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff use dispersion, only : expinit use gray_params, only : eqparam_type, prfparam_type, outparam_type, & - rtrparam_type, hcdparam_type, set_codepar, igrad, iequil, iprof, ieccd, & - iwarm, ipec, istpr0 + rtrparam_type, hcdparam_type, set_codepar, iequil, iprof, ieccd, & + iwarm, ipec, istpr0, igrad use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff use beamdata, only : pweight, print_projxyzt, rayi2jk use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, & @@ -54,8 +54,8 @@ contains real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,zeff,dersdst,derdnm,st,st0 - real(wp_) :: tau0,alphaabs0,didst0,ccci0 - real(wp_) :: tau,pow,dpdst,ddr,ddi,taumn,taumx + real(wp_) :: tau0,alphaabs0,dids0,ccci0 + real(wp_) :: tau,pow,ddr,ddi,taumn,taumx real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava real(wp_), dimension(3) :: xv,anv0,anv real(wp_), dimension(:,:), allocatable :: yw,ypw,gri @@ -63,7 +63,7 @@ contains integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,index_rt=1 logical :: ins_pl, somein, allout - real(wp_), dimension(:,:), allocatable :: psjki,tauv,alphav,ppabs,didst,ccci + real(wp_), dimension(:,:), allocatable :: psjki,tauv,alphav,ppabs,dids,ccci real(wp_), dimension(:), allocatable :: p0jk complex(wp_), dimension(:), allocatable :: ext, eyt integer, dimension(:), allocatable :: iiv @@ -104,7 +104,7 @@ contains call xgygcoeff(fghz,ak0,bres,xgcn) call launchangles2n(alpha0,beta0,xv0,anv0) call init_rtr(rtrp,yw,ypw,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) + psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) if(iwarm > 1) call expinit @@ -127,7 +127,7 @@ contains iox=iox0 sox=-1.0_wp_ if(iox==2) sox=1.0_wp_ - call vectinit(psjki,tauv,alphav,ppabs,didst,ccci,iiv) + call vectinit(psjki,tauv,alphav,ppabs,dids,ccci,iiv) call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri) psipol=psipol0 @@ -151,8 +151,7 @@ contains end do ! update position and grad -! if(igrad == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) - call gradi_upd(yw,ak0,xc,du1,gri,ggri) + if(igrad == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) ! test if the beam is completely out of the plsama allout = .true. @@ -163,8 +162,8 @@ contains call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), & psinv,dens,btot,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm) - if( abs(anpl) > anplth1) then ! anplth1=0.99_wp_ - if(abs(anpl) <= anplth2) then ! anplth2=1.05_wp_ + if( abs(anpl) > anplth1) then + if(abs(anpl) <= anplth2) then ierr=97 ! igrad=0 else @@ -175,16 +174,15 @@ contains ierr=0 end if - tekev=zero if(i==1) then tau0=zero alphaabs0=zero - didst0=zero + dids0=zero ccci0=zero else tau0=tauv(jk,i-1) alphaabs0=alphav(jk,i-1) - didst0=didst(jk,i-1) + dids0=dids(jk,i-1) ccci0=ccci(jk,i-1) end if zzm = xv(3)*0.01_wp_ @@ -196,13 +194,12 @@ contains if(ierr==0 .and. iwarm>0 .and. ins_pl .and. tau0<=taucr) then ! print*,i,jk,rayi2jk(jk),psinv,zzm,anpl tekev=temp(psinv) - if(tekev>zero) then - if (ieccd> 0) zeff=fzeff(psinv) - call alpha_effj(psinv,xg,yg,dens,tekev,zeff,ak0,bres,derdnm, & - anpl,anpr,sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierr) - iiv(jk)=i - end if + if (ieccd> 0) zeff=fzeff(psinv) + call alpha_effj(psinv,xg,yg,dens,tekev,zeff,ak0,bres,derdnm, & + anpl,anpr,sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierr) else + tekev=zero + zeff=zero alpha=zero didp=zero anprim=zero @@ -211,9 +208,12 @@ contains nhf=0 iokhawa=0 end if + if(nharm>0) iiv(jk)=i +! full storage required only for psjki,ppabs,ccci +! (jk,i) indexing can be removed from tauv,alphav,dids +! adding (jk) indexing to alphaabs0,tau0,dids0,ccci0 psjki(jk,i) = psinv - ! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s) tau=tau0+0.5_wp_*(alpha+alphaabs0)*dersdst*dst tauv(jk,i)=tau @@ -221,12 +221,11 @@ contains pow=p0jk(jk)*exp(-tau) !*exp(-tau1v(jk)) ppabs(jk,i)=p0jk(jk)-pow - dpdst=pow*alpha*dersdst - didst(jk,i)=didp*dpdst - ccci(jk,i)=ccci0+0.5_wp_*(didst0+didst(jk,i))*dst + dids(jk,i)=didp*pow*alpha + ccci(jk,i)=ccci0+0.5_wp_*(dids0+dids(jk,i))*dersdst*dst call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,ak0,anpl,anpr, & - anprim,dens,tekev,alphav(jk,i),tauv(jk,i),didst(jk,i),nhf,iokhawa, & + anprim,dens,tekev,alpha,tau,dids(jk,i),nhf,iokhawa, & index_rt,ddr,ddi) ! print error code @@ -274,7 +273,7 @@ contains write(*,'(a,f9.4)') 'I_tot (kA) = ',icd*1.0e3_wp_ ! compute power and current density profiles for all rays - call pec_init(ipec,sqrt(psinr)) + call pec_init(ipec) !,sqrt(psinr)) nnd=size(rhop_tab) allocate(jphi(nnd),pins(nnd),currins(nnd)) call spec(psjki,ppabs,ccci,iiv,pabs,icd,dpdv,jphi,jcd,pins,currins) @@ -291,7 +290,7 @@ contains ! ======= free memory BEGIN ====== call dealloc_beam(yw,ypw,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,didst,ccci,p0jk,ext,eyt,iiv) + psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) ! call unset_eqspl ! call unset_q ! call unset_rhospl @@ -302,13 +301,11 @@ contains end subroutine gray - subroutine vectinit(psjki,tauv,alphav,ppabs,didst,ccci,iiv) - use const_and_precisions, only : wp_, zero, one - use dispersion, only: expinit - use gray_params, only : iwarm + subroutine vectinit(psjki,tauv,alphav,ppabs,dids,ccci,iiv) + use const_and_precisions, only : wp_, zero implicit none ! arguments - real(wp_), dimension(:,:), intent(out) :: psjki,tauv,alphav,ppabs,didst,ccci + real(wp_), dimension(:,:), intent(out) :: psjki,tauv,alphav,ppabs,dids,ccci integer, dimension(:), intent(out) :: iiv !! common/external functions/variables ! integer :: jclosest @@ -324,9 +321,9 @@ contains tauv = zero alphav = zero ppabs = zero - didst = zero + dids = zero ccci = zero - iiv = one + iiv = 1 end subroutine vectinit @@ -337,7 +334,7 @@ contains ! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! use const_and_precisions, only : wp_,izero,zero,one,pi,half,two,degree,ui=>im use math, only : catand - use gray_params, only : ipol,idst + use gray_params, only : idst use beamdata, only : nray,nrayr,nrayth,rwmax implicit none ! arguments @@ -412,9 +409,9 @@ contains end if ! w01=sqrt(2.0_wp_/(ak0*ww1)) -! z01=-rci1/(rci1**2+ww1**2) +! d01=-rci1/(rci1**2+ww1**2) ! w02=sqrt(2.0_wp_/(ak0*ww2)) -! z02=-rci2/(rci2**2+ww2**2) +! d02=-rci2/(rci2**2+ww2**2) qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2 qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2 @@ -625,7 +622,7 @@ contains subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr) ! Runge-Kutta integrator use const_and_precisions, only : wp_ - use gray_params, only : igrad +! use gray_params, only : igrad use beamdata, only : h,hh,h6 implicit none real(wp_), intent(in) :: sox,bres,xgcn @@ -658,7 +655,6 @@ contains ! Compute right-hand side terms of the ray equations (dery) ! used in R-K integrator use const_and_precisions, only : wp_ - use gray_params, only : idst,igrad implicit none ! arguments real(wp_), dimension(6), intent(in) :: y @@ -686,7 +682,7 @@ contains xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update - use gray_params, only : igrad +! use gray_params, only : igrad implicit none ! arguments real(wp_), dimension(3), intent(in) :: xv,anv @@ -924,7 +920,7 @@ contains real(wp_), dimension(3), intent(out) :: bv,derxg,deryg real(wp_), dimension(3,3), intent(out) :: derbv ! local variables - integer :: iv,jv + integer :: jv real(wp_) :: xx,yy,zz real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm real(wp_), dimension(3) :: dbtot,bvc @@ -1076,7 +1072,7 @@ contains real(wp_), dimension(3,3), intent(in) :: ddgr,derbv real(wp_), dimension(6), intent(out) :: dery ! local variables - integer :: iv,jv + integer :: iv real(wp_) :: yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl,an2,an2s real(wp_) :: dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel real(wp_) :: derdom,dfdiadnpl,dfdiadxg,dfdiadyg,fdia,bdotgr !,vgm @@ -1192,8 +1188,7 @@ contains sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,ierr) use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_ use gray_params, only : iwarm,ilarm,ieccd,imx - use equilibrium, only : rmaxis,sgnbphi -! use beamdata, only : psjki,tauv,alphav,ppabs,didst,ccci,tau1v,pdjki,currj + use equilibrium, only : sgnbphi use dispersion, only : harmnumber, warmdisp use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl use magsurf_data, only : fluxval @@ -1219,6 +1214,10 @@ contains anprim=zero anprre=zero didp=zero + nhmin=0 + nhmax=0 + iokhawa=0 + ierr=0 if(tekev>zero) then ! absorption computation @@ -1334,7 +1333,7 @@ contains end subroutine set_pol subroutine print_output(i,jk,st,qj,xv,psinv,btot,ak0,anpl,anpr,anprim, & - dens,tekev,alpha,tau,didst,nhf,iokhawa,index_rt,ddr,ddi) + dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi) use const_and_precisions, only : degree,zero,one use equilibrium, only : frhotor use gray_params, only : istpl0 @@ -1344,9 +1343,9 @@ contains integer, intent(in) :: i,jk,nhf,iokhawa,index_rt real(wp_), dimension(3), intent(in) :: xv real(wp_), intent(in) :: st,qj,psinv,btot,ak0,anpl,anpr,anprim - real(wp_), intent(in) :: dens,tekev,alpha,tau,didst,ddr,ddi + real(wp_), intent(in) :: dens,tekev,alpha,tau,dids,ddr,ddi ! local variables - real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,dids + real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,didsn integer :: k stm=st*1.0e-2_wp_ @@ -1366,10 +1365,10 @@ contains end if akim=anprim*ak0*1.0e2_wp_ pt=exp(-tau) - dids=didst*1.0e2_wp_/qj + didsn=dids*1.0e2_wp_/qj write(4,'(30(1x,e16.8e3))') stm,rrm,zzm,phideg,psinv,rhot,dens,tekev, & - btot,anpr,anpl,akim,alpha,tau,pt,dids,dble(nhf),dble(iokhawa), & + btot,anpr,anpl,akim,alpha,tau,pt,didsn,dble(nhf),dble(iokhawa), & dble(index_rt),ddr end if ! central ray only end diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index 926eed6..eaf80b5 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -212,10 +212,10 @@ contains lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, & kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam ! local variables - integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr,u56,unit + integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr integer, dimension(kwrk) :: iwrk real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, & - ratio_cdbtor,ratio_pltor,fc,height,height2,r2iav,currp, & + ratio_cdbtor,ratio_pltor,fc,height,r2iav,currp, & area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, & dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, & shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, & @@ -231,9 +231,6 @@ contains ! common/external functions/variables real(wp_) :: fpolv,ddpsidrr,ddpsidzz - u56=get_free_unit() - open(file='f56.txt',unit=u56) - npsi=nnintp ninpr=(npsi-1)/10 npoints = 2*ncnt+1 @@ -478,7 +475,7 @@ contains end do end do - write(u56,*)' #rhop rhot || |Bmx| |Bmn| Area Vol |I_pl| fc ratJa ratJb' + write(56,*)' #rhop rhot || |Bmx| |Bmn| Area Vol |I_pl| fc ratJa ratJb' do jp=1,npsi if(jp.eq.npsi) then @@ -486,13 +483,11 @@ contains pstab(jp)=1.0_wp_ end if rhotjp=frhotor(rpstab(jp)) - write(u56,99) rpstab(jp),rhotjp,bav(jp),bmxpsi(jp),bmnpsi(jp), & + write(56,99) rpstab(jp),rhotjp,bav(jp),bmxpsi(jp),bmnpsi(jp), & varea(jp),vvol(jp),vcurrp(jp),vajphiav(jp),ffc(jp), & vratja(jp),vratjb(jp) end do - close(u56) - ! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs ! used for computations of dP/dV and J_cd iopt=0 diff --git a/src/main.f90 b/src/main.f90 index 0c697e1..3eb97e9 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -3,7 +3,7 @@ program gray_main 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 beams, only : read_beam0, read_beam1, read_beam2 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 @@ -69,14 +69,21 @@ program gray_main deallocate(xrad) end if ! re-scale input data - call tene_scal(terad,derad,prfp%factte,prfp%factne,eqp%factb,prfp%iscal,prfp%iprof) + 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 + select case (antp%ibeam) + case (2) +! to be completed: now 1st beamd always selected, iox read from table + call read_beam2(antp%filenm,1,antp%alpha,antp%beta,fghz,antp%iox,x0,y0,z0, & + w1,w2,ri1,ri2,phiw,phir) + case (1) + call read_beam1(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0, & + w1,w2,ri1,ri2,phiw,phir) + case default + call read_beam0(antp%filenm,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir) + end select alpha0=antp%alpha beta0=antp%beta p0mw=antp%power diff --git a/src/pec.f90 b/src/pec.f90 index 428b133..61d3eb6 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -15,18 +15,18 @@ contains implicit none ! arguments integer, intent(in) :: ipec - real(wp_), dimension(nnd), intent(in) :: rt_in + real(wp_), dimension(nnd), intent(in), optional :: rt_in ! local variables integer :: it real(wp_) :: drt,rt,rt1,rhop1 real(wp_) :: ratjai,ratjbi,ratjpli real(wp_) :: voli0,voli1,areai0,areai1 -! ipec positive build equidistant grid dimension nnd -! ipec negative read input grid +! rt_in present: read input grid +! else: build equidistant grid dimension nnd -! ipec=+/-1 rho_pol grid -! ipec=+/-2 rho_tor grid +! ipec=1 rho_pol grid +! ipec=2 rho_tor grid call dealloc_pec allocate(rhop_tab(nnd),rhot_tab(nnd),rtabpsi1(0:nnd),dvol(nnd),darea(nnd), & ratjav(nnd),ratjbv(nnd),ratjplv(nnd)) @@ -36,14 +36,16 @@ contains rtabpsi1(0) = zero do it=1,nnd - if(ipec > 0) then + if(present(rt_in)) then +! read radial grid from input + rt = rt_in(it) + if(it