From 46e36a5792643965c664355e2502cefba7c2ee44 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Mon, 23 Nov 2015 17:55:27 +0000 Subject: [PATCH] re-added missing output files, file units specified in dedicated module, added option iequil=3 for partially filled psi(R,z) grid on input: required for integration in JETTO --- Makefile | 22 +- src/beamdata.f90 | 299 ++++++------- src/beams.f90 | 34 +- src/conical.f90 | 17 +- src/coreprofiles.f90 | 3 +- src/dispersion.f90 | 15 +- src/eccd.f90 | 2 +- src/equilibrium.f90 | 931 ++++++++++++++++++++++++----------------- src/gray-externals.f90 | 475 --------------------- src/graycore.f90 | 672 ++++++++++++++++++++++++++--- src/limiter.f90 | 33 ++ src/magsurf_data.f90 | 285 ++++++------- src/main.f90 | 15 +- src/pec.f90 | 19 +- src/reflections.f90 | 73 ++-- src/units.f90 | 13 + srcjetto/gray.f | 152 +++++++ 17 files changed, 1709 insertions(+), 1351 deletions(-) create mode 100644 src/limiter.f90 create mode 100644 src/units.f90 create mode 100644 srcjetto/gray.f diff --git a/Makefile b/Makefile index 9723231..3aac4ad 100644 --- a/Makefile +++ b/Makefile @@ -4,9 +4,9 @@ EXE=gray # Objects list MAINOBJ=main.o OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \ - dierckx.o dispersion.o eccd.o eierf.o errcodes.o graycore.o gray-externals.o \ - gray_params.o equilibrium.o magsurf_data.o math.o minpack.o numint.o \ - pec.o polarization.o quadpack.o reflections.o simplespline.o utils.o + dierckx.o dispersion.o eccd.o eierf.o errcodes.o graycore.o \ + gray_params.o equilibrium.o limiter.o magsurf_data.o math.o minpack.o numint.o \ + pec.o polarization.o quadpack.o reflections.o simplespline.o units.o utils.o # Alternative search paths vpath %.f90 src @@ -29,12 +29,8 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) 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 coreprofiles.o \ - dispersion.o equilibrium.o errcodes.o gray-externals.o gray_params.o \ - pec.o polarization.o reflections.o utils.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 + dispersion.o eccd.o equilibrium.o errcodes.o gray_params.o \ + pec.o polarization.o limiter.o units.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 @@ -47,10 +43,10 @@ eccd.o: const_and_precisions.o conical.o dierckx.o errcodes.o magsurf_data.o \ eierf.o: const_and_precisions.o errcodes.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 gray_params.o +equilibrium.o: const_and_precisions.o dierckx.o limiter.o minpack.o \ + reflections.o simplespline.o utils.o gray_params.o magsurf_data.o: const_and_precisions.o gray_params.o equilibrium.o dierckx.o \ - reflections.o simplespline.o utils.o + reflections.o simplespline.o units.o utils.o math.o: const_and_precisions.o minpack.o: const_and_precisions.o numint.o: const_and_precisions.o @@ -58,7 +54,7 @@ pec.o: const_and_precisions.o beamdata.o equilibrium.o gray_params.o \ magsurf_data.o utils.o polarization.o: const_and_precisions.o quadpack.o: const_and_precisions.o -reflections.o: const_and_precisions.o utils.o +reflections.o: const_and_precisions.o limiter.o utils.o simplespline.o: const_and_precisions.o utils.o: const_and_precisions.o diff --git a/src/beamdata.f90 b/src/beamdata.f90 index 4357276..22a8f7c 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -2,13 +2,12 @@ module beamdata use const_and_precisions, only : wp_ implicit none - integer, save :: nray,nrayr,nrayth,nstep,jray1 + integer, save :: nray,nrayr,nrayth,nstep,jkray1 real(wp_), save :: dst,h,hh,h6,rwmax,twodr2 - integer, parameter :: nfileproj0 = 8, nfilew = 12 contains - subroutine init_rtr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & + subroutine init_btr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) use gray_params, only : rtrparam_type use const_and_precisions, only : zero,half,two @@ -21,137 +20,41 @@ contains dids0,ccci0,p0jk complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt integer, dimension(:), intent(out), allocatable :: iiv + + integer :: jray1 - dst=rtrparam%dst - h=dst - hh=h*half - h6=h/6.0_wp_ + dst = rtrparam%dst + h = dst + hh = h*half + h6 = h/6.0_wp_ - nrayr=rtrparam%nrayr - nrayth=rtrparam%nrayth - if(nrayr==1) nrayth=1 - nray=(nrayr-1)*nrayth+1 - - rwmax=rtrparam%rwmax + nrayr = rtrparam%nrayr + nrayth = rtrparam%nrayth + rwmax = rtrparam%rwmax + + if (nrayr==1) then + nrayth = 1 + jray1 = 1 + else + jray1 = 1 + max(nint((nrayr-1)/rwmax),1) + rwmax = dble(nrayr-1)/dble(jray1-1) + end if + nray = (nrayr-1)*nrayth + 1 + jkray1 = (jray1-2)*nrayth + 2 + if(nrayr>1) then twodr2 = two*(rwmax/(nrayr-1))**2 else twodr2 = two end if - + nstep=rtrparam%nstep call alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - 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(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - implicit none - real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & - gri,psjki,ppabs,ccci - real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri - real(wp_), dimension(:), intent(out), allocatable :: tau0,alphaabs0, & - dids0,ccci0,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,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,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), ppabs(nray,nstep), ccci(nray,nstep), & - tau0(nray), alphaabs0(nray), dids0(nray), ccci0(nray), & - p0jk(nray), ext(nray), eyt(nray), iiv(nray)) - end subroutine alloc_beam + end subroutine init_btr - subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - implicit none - real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & - gri,psjki,ppabs,ccci - real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri - real(wp_), dimension(:), intent(out), allocatable :: tau0,alphaabs0, & - dids0,ccci0,p0jk - complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt - integer, dimension(:), intent(out), allocatable :: iiv - - if (allocated(ywork)) deallocate(ywork) - if (allocated(ypwork)) deallocate(ypwork) - if (allocated(xc)) deallocate(xc) - if (allocated(du1)) deallocate(du1) - if (allocated(gri)) deallocate(gri) - if (allocated(ggri)) deallocate(ggri) - if (allocated(psjki)) deallocate(psjki) - if (allocated(ppabs)) deallocate(ppabs) - if (allocated(ccci)) deallocate(ccci) - if (allocated(tau0)) deallocate(tau0) - if (allocated(alphaabs0)) deallocate(alphaabs0) - if (allocated(dids0)) deallocate(dids0) - if (allocated(ccci0)) deallocate(ccci0) - if (allocated(p0jk)) deallocate(p0jk) - if (allocated(ext)) deallocate(ext) - if (allocated(eyt)) deallocate(eyt) - if (allocated(iiv)) deallocate(iiv) - end subroutine dealloc_beam subroutine pweight(p0,p0jk) ! power associated to jk-th ray p0jk(j) for total beam power p0 @@ -192,61 +95,123 @@ contains end do end subroutine pweight - subroutine print_projxyzt(st,ywrk,iproj) - use const_and_precisions, only : wp_, comp_huge, zero, one + + + function rayi2jk(i) result(jk) implicit none -! arguments - real(wp_), intent(in) :: st - real(wp_), dimension(:,:), intent(in) :: ywrk - integer, intent(in) :: iproj -! local variables - integer :: jk,jkz,nfile - integer, dimension(2) ::jkv - real(wp_), dimension(3) :: xv1,dir,dxv - real(wp_) :: dirm,rtimn,rtimx,csth1,snth1,csps1,snps1,xti,yti,zti,rti -! common/external functions/variables + integer, intent(in) :: i + integer, dimension(2) :: jk + integer :: ioff - nfile = nfileproj0 + iproj - - xv1 = ywrk(1:3,1) - dir = ywrk(4:6,1) - dirm = sqrt(dir(1)**2 + dir(2)**2 + dir(3)**2) - dir = dir/dirm - csth1 = dir(3) - snth1 = sqrt(one - csth1**2) - if(snth1 > zero) then - csps1=dir(2)/snth1 - snps1=dir(1)/snth1 + 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 - csps1=one - snps1=zero + jk = 1 end if + end function rayi2jk - if(iproj==0) then - jkz = nray - nrayth + 1 + + + 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 - jkz = 1 + jr = 1 end if + end function rayi2j - rtimn = comp_huge - rtimx = zero - do jk = jkz, nray - dxv = ywrk(1:3,jk) - xv1 - xti = dxv(1)*csps1 - dxv(2)*snps1 - yti =(dxv(1)*snps1 + dxv(2)*csps1)*csth1 - dxv(3)*snth1 - zti =(dxv(1)*snps1 + dxv(2)*csps1)*snth1 + dxv(3)*csth1 - rti = sqrt(xti**2 + yti**2) - jkv=rayi2jk(jk) - if(.not.(iproj==0 .and. jk==1)) & - write(nfile,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') st,jkv,xti,yti,zti,rti - if(iproj==1 .and. jkv(2)==nrayth) write(nfile,*) ' ' - if(rti>=rtimx .and. jkv(1)==nrayr) rtimx = rti - if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti - end do - write(nfile,*) ' ' - write(nfilew,'(3(1x,e16.8e3))') st,rtimn,rtimx - end subroutine print_projxyzt + 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(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) + implicit none + real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & + gri,psjki,ppabs,ccci + real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri + real(wp_), dimension(:), intent(out), allocatable :: tau0,alphaabs0, & + dids0,ccci0,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,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,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), ppabs(nray,nstep), ccci(nray,nstep), & + tau0(nray), alphaabs0(nray), dids0(nray), ccci0(nray), & + p0jk(nray), ext(nray), eyt(nray), iiv(nray)) + end subroutine alloc_beam + + + + subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) + implicit none + real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & + gri,psjki,ppabs,ccci + real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri + real(wp_), dimension(:), intent(out), allocatable :: tau0,alphaabs0, & + dids0,ccci0,p0jk + complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt + integer, dimension(:), intent(out), allocatable :: iiv + + if (allocated(ywork)) deallocate(ywork) + if (allocated(ypwork)) deallocate(ypwork) + if (allocated(xc)) deallocate(xc) + if (allocated(du1)) deallocate(du1) + if (allocated(gri)) deallocate(gri) + if (allocated(ggri)) deallocate(ggri) + if (allocated(psjki)) deallocate(psjki) + if (allocated(ppabs)) deallocate(ppabs) + if (allocated(ccci)) deallocate(ccci) + if (allocated(tau0)) deallocate(tau0) + if (allocated(alphaabs0)) deallocate(alphaabs0) + if (allocated(dids0)) deallocate(dids0) + if (allocated(ccci0)) deallocate(ccci0) + if (allocated(p0jk)) deallocate(p0jk) + if (allocated(ext)) deallocate(ext) + if (allocated(eyt)) deallocate(eyt) + if (allocated(iiv)) deallocate(iiv) + end subroutine dealloc_beam end module beamdata diff --git a/src/beams.f90 b/src/beams.f90 index 0343e1f..6f360e4 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -62,7 +62,7 @@ contains 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 + integer :: u,iopt,ier,nisteer,i,k,ii real(wp_) :: steer,dal real(wp_), dimension(:), allocatable :: alphastv,betastv,x00v,y00v, & z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v, & @@ -86,16 +86,7 @@ contains 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 + crci2(4*nisteer),cphi1(4*nisteer),cphi2(4*nisteer)) do i=1,nisteer read(u,*) steer,alphastv(i),betastv(i),x00v(i),y00v(i),z00v(i), & @@ -137,7 +128,7 @@ contains phiw=spli(cphi1,nisteer,k,dal) phir=spli(cphi2,nisteer,k,dal) else - write(*,*) ' alpha0 outside table range !!!' + ! alpha0 outside table range if(alpha0 >= alphastv(nisteer)) ii=nisteer if(alpha0 <= alphastv(1)) ii=1 beta0=betastv(ii) @@ -581,45 +572,45 @@ contains beta0 = ycoord0 SELECT CASE (in) CASE (1) - write(*,*) ' beta0 outside table range !!!' + ! 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 !!!' + ! 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 !!!' + ! 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 !!!' + ! alpha0 and beta0 outside table range xcoord0 = xvert(3) ycoord0 = yvert(3) ii = nxcoord+nycoord-1 CASE (5) - write(*,*) ' beta0 outside table range !!!' + ! 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 !!!' + ! alpha0 and beta0 outside table range xcoord0 = xvert(4) ycoord0 = yvert(4) ii = 2*nxcoord+nycoord-2 CASE (7) - write(*,*) ' alpha0 outside table range !!!' + ! 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 !!!!' + ! alpha0 and beta0 outside table range xcoord0 = xvert(1) ycoord0 = yvert(1) ii = 1 @@ -719,9 +710,6 @@ contains r=sqrt(xv(1)**2+xv(2)**2) ! phi=atan2(y,x) - print'(a,2f8.3)','alpha0, beta0 = ',alpha,beta - print'(a,4f8.3)','x00, y00, R00, z00 = ',xv(1:2),r,xv(3) - print*,' ' a = degree*alpha b = degree*beta ! diff --git a/src/conical.f90 b/src/conical.f90 index fc9c606..81b2fee 100644 --- a/src/conical.f90 +++ b/src/conical.f90 @@ -21,7 +21,7 @@ contains integer :: jp,j,n real(wp_), parameter :: rpi=1.7724538509055_wp_,pi2=0.63661977236758_wp_ real(wp_), parameter :: eps=1.0e-14_wp_ - integer, parameter :: nout=2,nmax=200 + integer, parameter :: nmax=200 ! complex(wp_) a,b,c,ti,r,rr,q,u,u0,u1,u2,uu complex(wp_) v0,v1,v2,vv,w(19) @@ -32,7 +32,7 @@ contains lm0=m == 0 lm1=m == 1 if(.not.(lm0 .or. lm1)) then - write(nout,"(1x,'fconic ... illegal value for m = ',i4)") m + write(*,"(1x,'fconic ... illegal value for m = ',i4)") m return end if fm=m @@ -202,7 +202,7 @@ contains do n=n+1 if(n > nmax) then - write(nout,200) x,tau,m + write(*,200) x,tau,m return end if rr=r @@ -256,7 +256,7 @@ contains if(abs(r-rr) < eps) exit end do if (n > nmax) then - write(nout,200) x,tau,m + write(*,200) x,tau,m return end if end if @@ -299,7 +299,6 @@ contains complex(wp_) :: v,h,r integer :: i,n real(wp_) :: x,t,a,c,d,e,f - integer, parameter :: nout=2 real(wp_), parameter :: pi=3.1415926535898_wp_ real(wp_), dimension(10), parameter :: b= & (/+8.3333333333333e-2_wp_, -2.7777777777778e-3_wp_, & @@ -311,7 +310,7 @@ contains x=real(z) t=aimag(z) if(-abs(x) == aint(x) .and. t == 0.0_wp_) then - write(nout,'(1x,f20.2)') x + write(*,'(1x,f20.2)') x clogam=(0.0_wp_,0.0_wp_) return end if @@ -431,7 +430,6 @@ contains real(wp_) :: besy0,besy1 logical :: l real(wp_) :: v,f,a,b,p,q - integer, parameter :: nout=2 ! entry besj0l(x) ! @@ -607,7 +605,7 @@ contains go to 3 ! 9 besjy=0.0_wp_ - write(nout,"(1x,'besjy ... non-positive argument x = ',e15.4)") x + write(*,"(1x,'besjy ... non-positive argument x = ',e15.4)") x end function besjy function besik(x) @@ -616,7 +614,6 @@ contains real(wp_) :: besik,ebesi0,besi0,ebesi1,besi1,ebesk0,besk0,ebesk1,besk1 logical :: l,e real(wp_) :: v,f,a,b,z - integer, parameter :: nout=2 ! entry ebesi0(x) ! @@ -845,7 +842,7 @@ contains if(x < 180.0_wp_) besik=exp(-x)*z return 9 besik=0.0_wp_ - write(nout,"(1x,'besik ... non-positive argument x = ',e15.4)") x + write(*,"(1x,'besik ... non-positive argument x = ',e15.4)") x end function besik ! ! routines for conical function: end diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index 95bcfe6..c23aeb3 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -64,10 +64,9 @@ contains ier=0 call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier) ddens=ffs(1) - 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) print*,'psin = ',psin,': DENSITY NEGATIVE ne=',dens ! if(dens < zero) then ! dens=zero ! ddens=zero diff --git a/src/dispersion.f90 b/src/dispersion.f90 index e4d87d1..2dca274 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -194,8 +194,6 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez) end if end do -! -! if(i.gt.imx) print*,' i>imx ',yg,errnpr,i ! if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).ge.huge(one).or. & abs(npr2).le.tiny(one)) then @@ -274,18 +272,11 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast) end do end do ! - select case(fast) - - case(2:3) + if (fast<4) then call hermitian(rr,yg,mu,npl,cr,fast,lrm) - - case(4:) + else call hermitian_2(rr,yg,mu,npl,cr,fast,lrm) - - case default - write(*,*) "unexpected value for flag 'fast' in dispersion:", fast - - end select + end if ! call antihermitian(ri,yg,mu,npl,ci,lrm) ! diff --git a/src/eccd.f90 b/src/eccd.f90 index 309e944..48e3bf5 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -154,7 +154,7 @@ contains call SpitzFuncCoeff(amu,Zeff,fc) nlm=nlmt call profil(0,tjp,njpt,tlm,nlmt,ch,ksp,ksp,rhop,nlm,chlm,ierr) - if(ierr>0) write(*,*) ' Hlambda profil =',ierr + if(ierr>0) print*,' Hlambda profil =',ierr npar=3+2*nlm allocate(eccdpar(npar)) eccdpar(1)=zeff diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 3c66e0b..8e96d80 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -274,14 +274,16 @@ contains bsign=int(sign(one,fpol(size(fpol)))) end subroutine eq_scal - subroutine set_eqspl(rv,zv,psin,psiwbrad,psinr,fpol,sspl,ssfp, & + subroutine set_eqspl(rv,zv,psin,psiwbrad,psinr,fpol,qpsi,sspl,ssfp, & r0,rax,zax,rbnd,zbnd,ixp) use const_and_precisions, only : zero,one use dierckx, only : regrid,coeff_parder,curfit,splev + use gray_params, only : iequil + use reflections, only : inside use utils, only : vmaxmin,vmaxmini implicit none ! arguments - real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol + real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol,qpsi real(wp_), dimension(:,:), intent(in) :: psin real(wp_), intent(in) :: psiwbrad real(wp_), intent(in) :: sspl,ssfp @@ -295,10 +297,11 @@ contains integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup real(wp_) :: sspln,fp,rax0,zax0,psinoptmp,psinxptmp real(wp_) :: rbmin,rbmax,rbinf,rbsup,r1,z1 + real(wp_), dimension(size(psinr)) :: rhotn real(wp_), dimension(1) :: fpoli - real(wp_), dimension(:), allocatable :: fvpsi,wf,wrk + real(wp_), dimension(:), allocatable :: rv1d,zv1d,fvpsi,wf,wrk integer, dimension(:), allocatable :: iwrk - integer :: ier,ixploc,info + integer :: ier,ixploc,info,i,j,ij ! compute array sizes and prepare working space arrays nr=size(rv) @@ -317,32 +320,95 @@ contains ! spline fitting/interpolation of psin(i,j) and derivatives +! allocate knots and spline coefficients arrays + if (allocated(tr)) deallocate(tr) + if (allocated(tz)) deallocate(tz) + if (allocated(cceq)) deallocate(cceq) + allocate(tr(nrest),tz(nzest),cceq(nrz)) + ! length in m !!! rmnm=rv(1) rmxm=rv(nr) zmnm=zv(1) zmxm=zv(nz) -! allocate knots and spline coefficients arrays - if (allocated(tr)) deallocate(tr) - if (allocated(tz)) deallocate(tz) - allocate(tr(nrest),tz(nzest),cceq(nrz)) -! allocate work arrays + + if (iequil>2) then +! data valid only inside boundary (psin=0 outside), e.g. source==ESCO +! presence of boundary anticipated here to filter invalid data + if(present(rbnd).and.present(zbnd)) then + nbnd=min(size(rbnd),size(zbnd)) + else + nbnd=0 + end if +! determine number of valid grid points + nrz=0 + do j=1,nz + do i=1,nr + if (nbnd.gt.0) then + if(.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle + else + if(psin(i,j).le.0.0d0) cycle + end if + nrz=nrz+1 + end do + end do +! store valid data + allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz)) + ij=0 + do j=1,nz + do i=1,nr + if (nbnd.gt.0) then + if(.not.inside(rbnd,zbnd,nbnd,rv(i),zv(j))) cycle + else + if(psin(i,j).le.0.0d0) cycle + end if + ij=ij+1 + rv1d(ij)=rv(i) + zv1d(ij)=zv(j) + fvpsi(ij)=psin(i,j) + wf(ij)=1.0d0 + end do + end do + +! fit as a scattered set of points +! use reduced number of knots to limit memory comsumption ? + nsr=nr/4+4 + nsz=nz/4+4 + sspln=sspl + call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln, & + rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier) +! if ier=-1 data are fitted using sspl=0 + if(ier.eq.-1) then + sspln=0.0_wp_ + nsr=nr/4+4 + nsz=nz/4+4 + call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln, & + rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier) + end if + deallocate(rv1d,zv1d,wf,fvpsi) +! reset nrz to the total number of grid points for next allocations + nrz=nr*nz + else +! iequil==2: data are valid on the full R,z grid + ! 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,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 - sspln=0.0_wp_ - call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & - kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & + allocate(fvpsi(nrz)) + fvpsi=reshape(transpose(psin),(/nrz/)) + sspln=sspl + call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & + 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 + sspln=0.0_wp_ + call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & + kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & + wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) + end if + deallocate(fvpsi) end if - deallocate(fvpsi) + ! compute spline coefficients for psi partial derivatives lw10 = nr*(ksplp-1) + nz*ksplp + nrz lw01 = nr*ksplp + nz*(ksplp-1) + nrz @@ -486,24 +552,188 @@ contains end if print'(a,f8.4)','BT_centr= ',btrcen print'(a,f8.4)','BT_axis = ',btaxis + +! compute rho_pol/rho_tor mapping based on input q profile + call setqphi_num(psinr,abs(qpsi),abs(psia),rhotn) + call set_rhospl(sqrt(psinr),rhotn) + end subroutine set_eqspl - subroutine unset_eqspl + + + subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, & + tx,nknt_x,ty,nknt_y,coeff,ierr) + use const_and_precisions, only : wp_, comp_eps + use dierckx, only : surfit implicit none - if(allocated(tr)) deallocate(tr) - if(allocated(tz)) deallocate(tz) - if(allocated(tfp)) deallocate(tfp) - if(allocated(cfp)) deallocate(cfp) - if(allocated(cceq)) deallocate(cceq) - if(allocated(cceq01)) deallocate(cceq01) - if(allocated(cceq10)) deallocate(cceq10) - if(allocated(cceq02)) deallocate(cceq02) - if(allocated(cceq20)) deallocate(cceq20) - if(allocated(cceq11)) deallocate(cceq11) - nsr=0 - nsz=0 - nsf=0 - end subroutine unset_eqspl +! arguments + integer, intent(in) :: n + real(wp_), dimension(n), intent(in) :: x, y, z + real(wp_), dimension(n), intent(in) :: w + integer, intent(in) :: kspl + real(wp_), intent(in) :: sspl + real(wp_), intent(in) :: xmin, xmax, ymin, ymax + real(wp_), dimension(nknt_x), intent(inout) :: tx + real(wp_), dimension(nknt_y), intent(inout) :: ty + integer, intent(inout) :: nknt_x, nknt_y + real(wp_), dimension(nknt_x*nknt_y), intent(out) :: coeff + integer, intent(out) :: ierr +! local variables + integer :: iopt + real(wp_) :: resid + integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest + real(wp_), dimension(:), allocatable :: wrk1, wrk2 + integer, dimension(:), allocatable :: iwrk + + nxest=nknt_x + nyest=nknt_y + ne = max(nxest,nyest) + + km = kspl+1 + u = nxest-km + v = nyest-km + b1 = kspl*min(u,v)+kspl+1 + b2 = (kspl+1)*min(u,v)+1 + lwrk1 = u*v*(2+b1+b2)+2*(u+v+km*(n+ne)+ne-2*kspl)+b2+1 + lwrk2 = u*v*(b2+1)+b2 + kwrk = n+(nknt_x-2*kspl-1)*(nknt_y-2*kspl-1) + allocate(wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk)) + + iopt=0 + call surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl, & + sspl,nxest,nyest,ne,comp_eps,nknt_x,tx,nknt_y,ty, & + coeff,resid,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ierr) + + deallocate(wrk1,wrk2,iwrk) + + end subroutine scatterspl + + + + subroutine setqphi_num(psinq,q,psia,rhotn) + use const_and_precisions, only : pi + use simplespline, only : difcs + implicit none +! arguments + real(wp_), dimension(:), intent(in) :: psinq,q + real(wp_), intent(in) :: psia + real(wp_), dimension(:), intent(out), optional :: rhotn +! local variables + real(wp_), dimension(size(q)) :: phit + real(wp_) :: dx + integer, parameter :: iopt=0 + integer :: k,ier + + nq=size(q) + if(allocated(psinr)) deallocate(psinr) + if(allocated(cq)) deallocate(cq) + allocate(psinr(nq),cq(nq,4)) + + psinr=psinq + call difcs(psinr,q,nq,iopt,cq,ier) + +! Toroidal flux phi = 2*pi*Integral q dpsi + phit(1)=0.0_wp_ + do k=1,nq-1 + dx=psinr(k+1)-psinr(k) + phit(k+1)=phit(k) + dx*(cq(k,1) + dx*(cq(k,2)/2.0_wp_ + & + dx*(cq(k,3)/3.0_wp_ + dx* cq(k,4)/4.0_wp_) ) ) + end do + phitedge=phit(nq) + if(present(rhotn)) rhotn(1:nq)=sqrt(phit/phitedge) + phitedge=2*pi*psia*phitedge + end subroutine setqphi_num + + + + subroutine set_equian(rax,zax,a,bax,qax,q1,qexp,n) + use const_and_precisions, only : pi,zero,one + implicit none +! arguments + real(wp_), intent(in) :: rax,zax,a,bax,qax,q1,qexp + integer, intent(in), optional :: n +! local variables + integer, parameter :: nqdef=101 + integer :: i + real(wp_) :: dr,fq0,fq1,qq,res,rn + real(wp_), dimension(:), allocatable :: rhotn,rhopn + + btaxis=bax + rmaxis=rax + zmaxis=zax + btrcen=bax + rcen=rax + aminor=a + zbinf=zmaxis-a + zbsup=zmaxis+a + q0=qax + qa=q1 + alq=qexp + sgnbphi=sign(one,bax) + + rmxm=rmaxis+aminor + rmnm=rmaxis-aminor + zmxm=zbsup + zmnm=zbinf + + if (present(n)) then + nq=n + else + nq=nqdef + end if + if (allocated(psinr)) deallocate(psinr) + allocate(psinr(nq),rhotn(nq),rhopn(nq)) + + dr=one/(nq-1) + rhotn(1)=zero + psinr(1)=zero + res=zero + fq0=zero + do i=2,nq + rn=(i-1)*dr + qq=q0+(q1-q0)*rn**qexp + fq1=rn/qq + res=res+0.5_wp_*(fq1+fq0)*dr + fq0=fq1 + rhotn(i)=rn + psinr(i)=res + end do + + phitedge=btaxis*aminor**2 ! temporary + psia=res*phitedge + phitedge=pi*phitedge ! final + psinr=psinr/res + rhopn=sqrt(psinr) + + call set_rhospl(rhopn,rhotn) + end subroutine set_equian + + + + subroutine set_rhospl(rhop,rhot) + use simplespline, only : difcs + implicit none +! arguments + real(wp_), dimension(:), intent(in) :: rhop, rhot +! local variables + integer, parameter :: iopt=0 + integer :: ier + + nrho=size(rhop) + + if(allocated(rhopr)) deallocate(rhopr) + if(allocated(rhotr)) deallocate(rhotr) + if(allocated(crhop)) deallocate(crhop) + if(allocated(crhot)) deallocate(crhot) + allocate(rhopr(nrho),rhotr(nrho),crhop(nrho,4),crhot(nrho,4)) + + rhopr=rhop + rhotr=rhot + call difcs(rhotr,rhopr,nrho,iopt,crhop,ier) + call difcs(rhopr,rhotr,nrho,iopt,crhot,ier) + end subroutine set_rhospl + + subroutine equinum_psi(rpsim,zpsim,psinv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz) @@ -558,6 +788,8 @@ contains end if end subroutine equinum_psi + + subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc) use dierckx, only : fpbisp implicit none @@ -579,6 +811,8 @@ contains derpsi=ffspl(1)*psia end subroutine sub_derpsi + + subroutine equinum_fpol(psinv,fpolv,dfpv) use dierckx, only : splev,splder implicit none @@ -605,354 +839,7 @@ contains end if end subroutine equinum_fpol - subroutine bfield(rpsim,zpsim,bphi,br,bz) - use gray_params, only : iequil - implicit none -! arguments - real(wp_), intent(in) :: rpsim,zpsim - real(wp_), intent(out), optional :: bphi,br,bz -! local variables - real(wp_) :: psin,fpol - if (iequil < 2) then - call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br) - if (present(bphi)) bphi=bphi/rpsim - else - call equinum_psi(rpsim,zpsim,psinv=bphi,dpsidr=bz,dpsidz=br) - if (present(bphi)) then - psin=bphi - call equinum_fpol(psin,fpol) - bphi=fpol/rpsim - end if - end if - if (present(br)) br=-br/rpsim - if (present(bz)) bz= bz/rpsim - end subroutine bfield - - subroutine setqphi_num(psinq,q,psia,rhotn) - use const_and_precisions, only : pi - use simplespline, only : difcs - implicit none -! arguments - real(wp_), dimension(:), intent(in) :: psinq,q - real(wp_), intent(in) :: psia - real(wp_), dimension(:), intent(out), optional :: rhotn -! local variables - real(wp_), dimension(size(q)) :: phit - real(wp_) :: dx - integer, parameter :: iopt=0 - integer :: k,ier - - nq=size(q) - if(allocated(psinr)) deallocate(psinr) - if(allocated(cq)) deallocate(cq) - allocate(psinr(nq),cq(nq,4)) - - psinr=psinq - call difcs(psinr,q,nq,iopt,cq,ier) - -! Toroidal flux phi = 2*pi*Integral q dpsi - phit(1)=0.0_wp_ - do k=1,nq-1 - dx=psinr(k+1)-psinr(k) - phit(k+1)=phit(k) + dx*(cq(k,1) + dx*(cq(k,2)/2.0_wp_ + & - dx*(cq(k,3)/3.0_wp_ + dx* cq(k,4)/4.0_wp_) ) ) - end do - phitedge=phit(nq) - if(present(rhotn)) rhotn(1:nq)=sqrt(phit/phitedge) - phitedge=2*pi*psia*phitedge - end subroutine setqphi_num - - subroutine unset_q - implicit none - - if(allocated(psinr)) deallocate(psinr) - if(allocated(cq)) deallocate(cq) - nq=0 - end subroutine unset_q - - function fq(psin) - use const_and_precisions, only : wp_ - use gray_params, only : iequil - use simplespline, only :spli - use utils, only : locate - implicit none -! arguments - real(wp_), intent(in) :: psin - real(wp_) :: fq -! local variables - integer :: i - real(wp_) :: dps,rn - - 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) - use simplespline, only : difcs - implicit none -! arguments - real(wp_), dimension(:), intent(in) :: rhop, rhot -! local variables - integer, parameter :: iopt=0 - integer :: ier - - nrho=size(rhop) - - if(allocated(rhopr)) deallocate(rhopr) - if(allocated(rhotr)) deallocate(rhotr) - if(allocated(crhop)) deallocate(crhop) - if(allocated(crhot)) deallocate(crhot) - allocate(rhopr(nrho),rhotr(nrho),crhop(nrho,4),crhot(nrho,4)) - - rhopr=rhop - rhotr=rhot - call difcs(rhotr,rhopr,nrho,iopt,crhop,ier) - call difcs(rhopr,rhotr,nrho,iopt,crhot,ier) - end subroutine set_rhospl - - subroutine unset_rhospl - implicit none - - if(allocated(rhopr)) deallocate(rhopr) - if(allocated(rhotr)) deallocate(rhotr) - if(allocated(crhop)) deallocate(crhop) - if(allocated(crhot)) deallocate(crhot) - nrho=0 - end subroutine unset_rhospl - - function frhopol(rhot) - use utils, only : locate - use simplespline, only : spli - implicit none -! arguments - real(wp_), intent(in) :: rhot - real(wp_) :: frhopol -! local variables - integer :: i - real(wp_) :: dr - - call locate(rhotr,nrho,rhot,i) - i=min(max(1,i),nrho-1) - dr=rhot-rhotr(i) - 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 - implicit none -! arguments - real(wp_), intent(in) :: rhop - real(wp_) :: frhotor -! local variables - integer :: i - real(wp_) :: dr - - call locate(rhopr,nrho,rhop,i) - i=min(max(1,i),nrho-1) - dr=rhop-rhopr(i) - frhotor=spli(crhot,nrho,i,dr) - end function frhotor - - subroutine points_ox(rz,zz,rf,zf,psinvf,info) - use const_and_precisions, only : comp_eps - use minpack, only : hybrj1 - implicit none -! local constants - integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 -! arguments - real(wp_), intent(in) :: rz,zz - real(wp_), intent(out) :: rf,zf,psinvf - integer, intent(out) :: info -! local variables - real(wp_) :: tol - real(wp_), dimension(n) :: xvec,fvec - real(wp_), dimension(lwa) :: wa - real(wp_), dimension(ldfjac,n) :: fjac - - xvec(1)=rz - xvec(2)=zz - tol = sqrt(comp_eps) - call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) - if(info.gt.1) then - print'(a,i2,a,2f8.4)',' info subr points_ox =',info, & - ' O/X coord.',xvec - end if - rf=xvec(1) - zf=xvec(2) - call equinum_psi(rf,zf,psinvf) - end subroutine points_ox - - subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) - implicit none -! arguments - integer, intent(in) :: n,iflag,ldfjac - real(wp_), dimension(n), intent(in) :: x - real(wp_), dimension(n), intent(inout) :: fvec - real(wp_), dimension(ldfjac,n), intent(inout) :: fjac -! local variables - real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz - - select case(iflag) - case(1) - call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz) - fvec(1) = dpsidr/psia - fvec(2) = dpsidz/psia - case(2) - call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, & - ddpsidrz=ddpsidrz) - fjac(1,1) = ddpsidrr/psia - fjac(1,2) = ddpsidrz/psia - fjac(2,1) = ddpsidrz/psia - fjac(2,2) = ddpsidzz/psia - case default - print*,'iflag undefined' - end select - end subroutine fcnox - - subroutine points_tgo(rz,zz,rf,zf,psin0,info) - use const_and_precisions, only : comp_eps - use minpack, only : hybrj1mv - implicit none -! local constants - integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 -! arguments - real(wp_), intent(in) :: rz,zz,psin0 - real(wp_), intent(out) :: rf,zf - integer, intent(out) :: info -! local variables - real(wp_) :: tol - real(wp_), dimension(n) :: xvec,fvec,f0 - real(wp_), dimension(lwa) :: wa - real(wp_), dimension(ldfjac,n) :: fjac - - xvec(1)=rz - xvec(2)=zz - f0(1)=psin0 - f0(2)=0.0_wp_ - tol = sqrt(comp_eps) - call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa) - if(info.gt.1) then - print'(a,i2,a,5f8.4)',' info subr points_tgo =',info, & - ' R,z coord.',xvec,rz,zz,psin0 - end if - rf=xvec(1) - zf=xvec(2) - end - - subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag) - use const_and_precisions, only : wp_ - implicit none -! arguments - integer, intent(in) :: n,ldfjac,iflag - real(wp_), dimension(n), intent(in) :: x,f0 - real(wp_), dimension(n), intent(inout) :: fvec - real(wp_), dimension(ldfjac,n), intent(inout) :: fjac -! internal variables - real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz - - select case(iflag) - case(1) - call equinum_psi(x(1),x(2),psinv,dpsidr) - fvec(1) = psinv-f0(1) - fvec(2) = dpsidr/psia-f0(2) - case(2) - call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, & - ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) - fjac(1,1) = dpsidr/psia - fjac(1,2) = dpsidz/psia - fjac(2,1) = ddpsidrr/psia - fjac(2,2) = ddpsidrz/psia - case default - print*,'iflag undefined' - end select - end subroutine fcntgo - - subroutine set_equian(rax,zax,a,bax,qax,q1,qexp,n) - use const_and_precisions, only : pi,zero,one - implicit none -! arguments - real(wp_), intent(in) :: rax,zax,a,bax,qax,q1,qexp - integer, intent(in), optional :: n -! local variables - integer, parameter :: nqdef=101 - integer :: i - real(wp_) :: dr,fq0,fq1,qq,res,rn - real(wp_), dimension(:), allocatable :: rhotn,rhopn - - btaxis=bax - rmaxis=rax - zmaxis=zax - btrcen=bax - rcen=rax - aminor=a - zbinf=zmaxis-a - zbsup=zmaxis+a - q0=qax - qa=q1 - alq=qexp - sgnbphi=sign(one,bax) - - if (present(n)) then - nq=n - else - nq=nqdef - end if - if (allocated(psinr)) deallocate(psinr) - allocate(psinr(nq),rhotn(nq),rhopn(nq)) - - dr=one/(nq-1) - rhotn(1)=zero - psinr(1)=zero - res=zero - fq0=zero - do i=2,nq - rn=(i-1)*dr - qq=q0+(q1-q0)*rn**qexp - fq1=rn/qq - res=res+0.5_wp_*(fq1+fq0)*dr - fq0=fq1 - rhotn(i)=rn - psinr(i)=res - end do - - phitedge=btaxis*aminor**2 ! temporary - psia=res*phitedge - phitedge=pi*phitedge ! final - psinr=psinr/res - rhopn=sqrt(psinr) - - call set_rhospl(rhopn,rhotn) - end subroutine set_equian subroutine equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz) @@ -1013,6 +900,119 @@ contains + function frhopol(rhot) + use utils, only : locate + use simplespline, only : spli + implicit none +! arguments + real(wp_), intent(in) :: rhot + real(wp_) :: frhopol +! local variables + integer :: i + real(wp_) :: dr + + call locate(rhotr,nrho,rhot,i) + i=min(max(1,i),nrho-1) + dr=rhot-rhotr(i) + 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 + implicit none +! arguments + real(wp_), intent(in) :: rhop + real(wp_) :: frhotor +! local variables + integer :: i + real(wp_) :: dr + + call locate(rhopr,nrho,rhop,i) + i=min(max(1,i),nrho-1) + dr=rhop-rhopr(i) + frhotor=spli(crhot,nrho,i,dr) + end function frhotor + + + + function fq(psin) + use const_and_precisions, only : wp_ + use gray_params, only : iequil + use simplespline, only :spli + use utils, only : locate + implicit none +! arguments + real(wp_), intent(in) :: psin + real(wp_) :: fq +! local variables + integer :: i + real(wp_) :: dps,rn + + 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 bfield(rpsim,zpsim,bphi,br,bz) + use gray_params, only : iequil + implicit none +! arguments + real(wp_), intent(in) :: rpsim,zpsim + real(wp_), intent(out), optional :: bphi,br,bz +! local variables + real(wp_) :: psin,fpol + + if (iequil < 2) then + call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br) + if (present(bphi)) bphi=bphi/rpsim + else + call equinum_psi(rpsim,zpsim,psinv=bphi,dpsidr=bz,dpsidz=br) + if (present(bphi)) then + psin=bphi + call equinum_fpol(psin,fpol) + bphi=fpol/rpsim + end if + end if + if (present(br)) br=-br/rpsim + if (present(bz)) bz= bz/rpsim + end subroutine bfield + + + subroutine tor_curr(rpsim,zpsim,ajphi) use const_and_precisions, only : wp_,ccj=>mu0inv use gray_params, only : iequil @@ -1067,7 +1067,7 @@ contains r1=zeroc(1) r2=zeroc(2) end if - end subroutine psi_raxis + end subroutine psi_raxis @@ -1082,4 +1082,167 @@ contains call tor_curr(r2,zmaxis,ajphi) end subroutine tor_curr_psi + + + subroutine points_ox(rz,zz,rf,zf,psinvf,info) + use const_and_precisions, only : comp_eps + use minpack, only : hybrj1 + implicit none +! local constants + integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 +! arguments + real(wp_), intent(in) :: rz,zz + real(wp_), intent(out) :: rf,zf,psinvf + integer, intent(out) :: info +! local variables + real(wp_) :: tol + real(wp_), dimension(n) :: xvec,fvec + real(wp_), dimension(lwa) :: wa + real(wp_), dimension(ldfjac,n) :: fjac + + xvec(1)=rz + xvec(2)=zz + tol = sqrt(comp_eps) + call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) + if(info.gt.1) then + print'(a,i2,a,2f8.4)',' info subr points_ox =',info, & + ' O/X coord.',xvec + end if + rf=xvec(1) + zf=xvec(2) + call equinum_psi(rf,zf,psinvf) + end subroutine points_ox + + + + subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) + implicit none +! arguments + integer, intent(in) :: n,iflag,ldfjac + real(wp_), dimension(n), intent(in) :: x + real(wp_), dimension(n), intent(inout) :: fvec + real(wp_), dimension(ldfjac,n), intent(inout) :: fjac +! local variables + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz + + select case(iflag) + case(1) + call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz) + fvec(1) = dpsidr/psia + fvec(2) = dpsidz/psia + case(2) + call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, & + ddpsidrz=ddpsidrz) + fjac(1,1) = ddpsidrr/psia + fjac(1,2) = ddpsidrz/psia + fjac(2,1) = ddpsidrz/psia + fjac(2,2) = ddpsidzz/psia + case default + print*,'iflag undefined' + end select + end subroutine fcnox + + + + subroutine points_tgo(rz,zz,rf,zf,psin0,info) + use const_and_precisions, only : comp_eps + use minpack, only : hybrj1mv + implicit none +! local constants + integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 +! arguments + real(wp_), intent(in) :: rz,zz,psin0 + real(wp_), intent(out) :: rf,zf + integer, intent(out) :: info +! local variables + real(wp_) :: tol + real(wp_), dimension(n) :: xvec,fvec,f0 + real(wp_), dimension(lwa) :: wa + real(wp_), dimension(ldfjac,n) :: fjac + + xvec(1)=rz + xvec(2)=zz + f0(1)=psin0 + f0(2)=0.0_wp_ + tol = sqrt(comp_eps) + call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa) + if(info.gt.1) then + print'(a,i2,a,5f8.4)',' info subr points_tgo =',info, & + ' R,z coord.',xvec,rz,zz,psin0 + end if + rf=xvec(1) + zf=xvec(2) + end + + + + subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag) + use const_and_precisions, only : wp_ + implicit none +! arguments + integer, intent(in) :: n,ldfjac,iflag + real(wp_), dimension(n), intent(in) :: x,f0 + real(wp_), dimension(n), intent(inout) :: fvec + real(wp_), dimension(ldfjac,n), intent(inout) :: fjac +! internal variables + real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz + + select case(iflag) + case(1) + call equinum_psi(x(1),x(2),psinv,dpsidr) + fvec(1) = psinv-f0(1) + fvec(2) = dpsidr/psia-f0(2) + case(2) + call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, & + ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) + fjac(1,1) = dpsidr/psia + fjac(1,2) = dpsidz/psia + fjac(2,1) = ddpsidrr/psia + fjac(2,2) = ddpsidrz/psia + case default + print*,'iflag undefined' + end select + end subroutine fcntgo + + + + subroutine unset_eqspl + implicit none + if(allocated(tr)) deallocate(tr) + if(allocated(tz)) deallocate(tz) + if(allocated(tfp)) deallocate(tfp) + if(allocated(cfp)) deallocate(cfp) + if(allocated(cceq)) deallocate(cceq) + if(allocated(cceq01)) deallocate(cceq01) + if(allocated(cceq10)) deallocate(cceq10) + if(allocated(cceq02)) deallocate(cceq02) + if(allocated(cceq20)) deallocate(cceq20) + if(allocated(cceq11)) deallocate(cceq11) + nsr=0 + nsz=0 + nsf=0 + end subroutine unset_eqspl + + + + subroutine unset_q + implicit none + + if(allocated(psinr)) deallocate(psinr) + if(allocated(cq)) deallocate(cq) + nq=0 + end subroutine unset_q + + + + subroutine unset_rhospl + implicit none + + if(allocated(rhopr)) deallocate(rhopr) + if(allocated(rhotr)) deallocate(rhotr) + if(allocated(crhop)) deallocate(crhop) + if(allocated(crhot)) deallocate(crhot) + nrho=0 + end subroutine unset_rhospl + end module equilibrium diff --git a/src/gray-externals.f90 b/src/gray-externals.f90 index 640fa8c..4efafae 100644 --- a/src/gray-externals.f90 +++ b/src/gray-externals.f90 @@ -414,478 +414,3 @@ !99 format(24(1x,e16.8e3)) !111 format(3i5,20(1x,e16.8e3)) ! end - - - - subroutine prfile - implicit none - write(4,*)' #sst R z phi psi rhot ne Te Btot '// & - 'Nperp Npl ki alpha tau Pt dIds nh iohkw index_rt ddr' - write(8,*) ' #istep j k xt yt zt rt' - write(9,*) ' #istep j k xt yt zt rt' - write(17,*) ' #sst Dr_Nr1 Di_Nr1' - write(33,*) ' #i jk sst x y R z psi tauv Npl alpha index_rt' - write(12,*) ' #sst w1 w2' - write(7,*)'#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav '// & - 'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt '// & - 'Jphimx dPdVmx drhotj drhotp' - write(48,*) '#rhop rhot Jphi Jcdb dPdV Icdins Pins' -! write(66,*) "# psipol0 chipol0 powrfl" - end subroutine prfile - - - - subroutine print_prof - use const_and_precisions, only : wp_ - use equilibrium, only : psinr,nq,fq,frhotor,tor_curr_psi - use coreprofiles, only : density, temp - implicit none -! local constants - real(wp_), parameter :: eps=1.e-4_wp_ -! local variables - integer :: i - real(wp_) :: psin,rhop,rhot,ajphi,te,qq - real(wp_) :: dens,ddens - - write(55,*) ' #psi rhot ne Te q Jphi' - do i=1,nq - psin=psinr(i) - rhop=sqrt(psin) - - call density(psin,dens,ddens) - te=temp(psin) - qq=fq(psin) - rhot=frhotor(rhop) - call tor_curr_psi(max(eps,psin),ajphi) - write(55,"(12(1x,e12.5))") psin,rhot,dens,te,qq,ajphi*1.e-6_wp_ - end do - end subroutine print_prof - - subroutine print_prof_an - use const_and_precisions, only : wp_ - use coreprofiles, only : density, temp - use equilibrium, only : frhotor - implicit none -! local constants - integer, parameter :: nst=51 -! local variables - integer :: i - real(wp_) :: psin,rhop,rhot,te - real(wp_) :: dens,ddens - - write(55,*) ' #psi rhot ne Te' - do i=1,nst - psin=dble(i-1)/dble(nst-1) - rhop=sqrt(psin) - rhot=frhotor(rhop) - call density(psin,dens,ddens) - te=temp(psin) - write(55,"(12(1x,e12.5))") psin,rhot,dens,te - end do - end subroutine print_prof_an - - - - 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,contours_psi - use utils, only : locate, intlin - implicit none -! arguments - integer, intent(in) :: nq - real(wp_), dimension(nq), intent(in) :: psinq,qpsi - real(wp_) :: qval -! local variables - integer :: ncnt,i1,ipr - real(wp_) :: rup,zup,rlw,zlw,rhot,psival - real(wp_), dimension(npoints) :: rcn,zcn - - ncnt=(npoints-1)/2 - -! locate psi surface for q=qval - call locate(abs(qpsi),nq,qval,i1) - if (i1>0.and.i1 h .or. & - ah > 0.0_wp_ .and. a(jm) <= h) then - ix=ix+1 - lx(ix)=-j - end if - if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & - ah > 0.0_wp_ .and. a(j-1) <= h) then - ix=ix+1 - lx(ix)=j - end if - end do - end do - - do jm=nr,mxr,nrqmax - j = jm + nrqmax - ah=a(j)-h - if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & - ah > 0.0_wp_ .and. a(j-1) <= h) then - ix=ix+1 - lx(ix)=j - end if - if (ah <= 0.0_wp_ .and. a(jm) > h .or. & - ah > 0.0_wp_ .and. a(jm) <= h) then - ix=ix+1 - lx(ix)=-j - end if - end do - - do jm=1,mxr,nrqmax - j = jm + nrqmax - if (a(j) <= h .and. a(jm) > h .or. & - a(j) > h .and. a(jm) <= h) then - ix=ix+1 - lx(ix) =-j - end if - end do - - do j=2,nr - if (a(j) <= h .and. a(j-1) > h .or. & - a(j) > h .and. a(j-1) <= h) then - ix=ix+1 - lx(ix)=j - end if - end do - - if(ix<=0) return - -bb: do - in=ix - jx=lx(in) - jfor=0 - lda=1 - ldb=2 - - do - if(jx<0) then - jabs=-jx - jnb = jabs - nrqmax - else - jabs=jx - jnb=jabs-1 - end if - - adn=a(jabs)-a(jnb) - if(adn/=0) px=(a(jabs)-h)/adn - kx = (jabs - 1) / nrqmax - ikx = jabs - nrqmax * kx - 1 - - if(jx<0) then - x = drgrd * ikx - y = dzgrd * (kx - px) - else - x = drgrd * (ikx - px) - y = dzgrd * kx - end if - - icount = icount + 1 - rcon(icount) = x + rqgrid(1) - zcon(icount) = y + zqgrid(1) - mpl= icount - itm = 1 - ja(1,1) = jabs + nrqmax - j=1 - - if(jx<=0) then - ja(1,1) = -jabs-1 - j=2 - end if - - ja(2,1) = -ja(1,1) - ja(3,1) = -jx + 1 - nrqmax - ja(3,2) = -jx - ja(j,2) = jabs - nrqmax - k= 3-j - ja(k,2) = 1-jabs - - if (kx<=0 .or. ikx<=0) then - lda=1 - ldb=lda - else if (ikx + 1 - nr >= 0 .and. jx <= 0) then - lda=2 - ldb=lda - else if(jfor/=0) then - lda=2 - do i=1,3 - if(jfor==ja(i,2)) then - lda=1 - exit - end if - end do - ldb=lda - end if - - flag1=.false. - aa: do k=1,3 - do l=lda,ldb - do i=1,ix - if(lx(i)==ja(k,l)) then - itm=itm+1 - inext= i - if(jfor/=0) exit aa - if(itm .gt. 3) then - flag1=.true. - exit aa - end if - end if - end do - end do - end do aa - - if(.not.flag1) then - lx(in)=0 - if(itm .eq. 1) exit - end if - - jfor=jx - jx=lx(inext) - in = inext - end do - - do - if(lx(ix)/=0) then - if(mpl>=4) then - ncon = ncon + 1 - npts(ncon) = icount - iclast - iclast = icount - end if - exit - end if - ix= ix-1 - if(ix<=0) exit bb - end do - - end do bb - - if(mpl >= 4) then - ncon = ncon + 1 - npts(ncon) = icount - iclast - iclast = icount - end if - end subroutine cniteq - - - - logical function inside_plasma(rrm,zzm) - use const_and_precisions, only : wp_, zero, one - use gray_params, only : iequil - use coreprofiles, only : psdbnd - use equilibrium, only : zbinf,zbsup,equinum_psi,equian - implicit none -! arguments - real(wp_), intent(in) :: rrm,zzm -! local variables - real(wp_) :: psinv - - if(iequil.eq.1) then - call equian(rrm,zzm,psinv) - else - call equinum_psi(rrm,zzm,psinv) - end if - - inside_plasma = (psinv >= zero .and. psinv < psdbnd) .and. & - (psinv >= one .or. (zzm >= zbinf .and. zzm <= zbsup)) - end function inside_plasma - - - - subroutine vacuum_rt(xv0,anv0,xvend,dstvac,ivac) - use const_and_precisions, only : wp_ - use reflections, only : inters_linewall,inside,rlim,zlim,nlim - use beamdata, only : dst - implicit none -! arguments - real(wp_), dimension(3), intent(in) :: xv0,anv0 - real(wp_), dimension(3), intent(out) :: xvend - real(wp_), intent(out) :: dstvac - integer, intent(out) :: ivac -! local variables - integer :: i - real(wp_) :: st,rrm,zzm,smax - real(wp_), dimension(3) :: walln - logical :: plfound -! common/external functions/variables - logical, external :: inside_plasma - -! ivac=1 plasma hit before wall reflection -! ivac=2 wall hit before plasma -! ivac=-1 vessel (and thus plasma) never crossed - - call inters_linewall(xv0/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), & - nlim,smax,walln) - smax=smax*1.0e2_wp_ - rrm=1.0e-2_wp_*sqrt(xv0(1)**2+xv0(2)**2) - zzm=1.0e-2_wp_*xv0(3) - if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then - ! first wall interface is outside-inside - if (dot_product(walln,walln) 1) call expinit @@ -111,16 +103,17 @@ contains ! ======= set environment END ====== ! ======= pre-proc prints BEGIN ====== + call print_headers +! print psi surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout + call print_surfq((/1.5_wp_,2.0_wp_/)) +! print + print*,' ' + print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0 + print'(a,4f8.3)','x00, y00, z00 = ',xv0 ! print Btot=Bres ! print ne, Te, q, Jphi versus psi, rhop, rhot - if (iequil<2) 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 - call prfile + call print_bres(bres) + call print_prof ! ======= pre-proc prints END ====== ! ======= main loop BEGIN ====== @@ -233,42 +226,48 @@ contains ! ======= main loop END ====== ! ======= post-proc BEGIN ====== -! print all ray positions in local reference system - if(nray > 1) call print_projxyzt(st,yw,1) ! print final results on screen write(*,*) write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs write(*,'(a,f9.4)') 'I_tot (kA) = ',icd*1.0e3_wp_ +! print all ray positions in local reference system + if(nray > 1) call print_projxyzt(st,yw,1) ! compute power and current density profiles for all rays 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) - - call postproc_profiles(pabs,icd,rhot_tab,dpdv,jphi, & - rhotpav,drhotpav,rhotjava,drhotjava) - ! print power and current density profiles - do i=1,nnd - write(48,'(7(1x,e16.8e3))') rhop_tab(i),rhot_tab(i), & - jphi(i),jcd(i),dpdv(i),currins(i),pins(i) - end do + call print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt) +! compute profiles width + call postproc_profiles(pabs,icd,rhot_tab,dpdv,jphi, & + rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, & + rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) +! print 0D results + call print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, & + drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, & + st,psipol,chipol,index_rt) + + ! ======= post-proc END ====== ! ======= free memory BEGIN ====== - call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) ! call unset_eqspl ! call unset_q ! call unset_rhospl ! call unset_prfspl - call dealloc_pec - deallocate(jphi,pins,currins) +! call unset_lim +! call dealloc_surfvec +! call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & +! tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) +! call dealloc_pec +! deallocate(jphi,pins,currins) ! ======= free memory END ====== - end subroutine gray + end subroutine gray_main + subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) @@ -300,6 +299,7 @@ contains end subroutine vectinit + subroutine ic_gb(xv0c,anv0c,ak0,wcsi,weta,rcicsi,rcieta,phiw,phir, & ywrk0,ypwrk0,xc0,du10,gri,ggri) ! beam tracing initial conditions igrad=1 @@ -487,7 +487,7 @@ contains dy0t = dcsiw*snphiw + detaw*csphiw x0t = uj(j)*dx0t y0t = uj(j)*dy0t - z0t = -half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t + z0t = -(half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t) dx0 = x0t*csps + snps*(y0t*csth + z0t*snth) dy0 = -x0t*snps + csps*(y0t*csth + z0t*snth) @@ -595,10 +595,13 @@ contains ddr = anx**2 + any**2 + anz**2 - an20 ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt) + call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,ak0,zero,zero,zero, & + zero,zero,zero,zero,zero,0,0,1,ddr,ddi) ! st=0, index_rt=1, Btot=0, psin=-1 end do - write(17,'(3(1x,e16.8e3))') zero,ddr,ddi end subroutine ic_gb + + subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr) ! Runge-Kutta integrator use const_and_precisions, only : wp_ @@ -631,6 +634,8 @@ contains y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4) end subroutine rkstep + + subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery) ! Compute right-hand side terms of the ray equations (dery) ! used in R-K integrator @@ -658,6 +663,7 @@ contains end subroutine rhs + subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr) ! Compute right-hand side terms of the ray equations (dery) @@ -696,6 +702,7 @@ contains end subroutine ywppla_upd + subroutine gradi_upd(ywrk,ak0,xc,du1,gri,ggri) use const_and_precisions, only : wp_,zero,half use beamdata, only : nray,nrayr,nrayth,twodr2 @@ -840,6 +847,8 @@ contains end subroutine gradi_upd + + subroutine solg0(dxv1,dxv2,dxv3,dgg) ! solution of the linear system of 3 eqs : dgg . dxv = dff ! input vectors : dxv1, dxv2, dxv3, dff @@ -895,6 +904,7 @@ contains end subroutine solg3 + subroutine plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv, & xg,yg,derxg,deryg,ajphi) use const_and_precisions, only : wp_,zero,pi,ccj=>mu0inv @@ -1048,6 +1058,7 @@ contains end subroutine plas_deriv + subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & dery,anpl,anpr,ddr,ddi,dersdst,derdnm) use const_and_precisions, only : wp_,zero,one,half,two @@ -1173,6 +1184,7 @@ contains end subroutine disp_deriv + subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,ierr) use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_ @@ -1325,12 +1337,519 @@ contains end do end subroutine set_pol + + +! logical function inside_plasma(rrm,zzm) +! use const_and_precisions, only : wp_, zero, one +! use gray_params, only : iequil +! use equilibrium, only : equian,equinum_psi,zbinf,zbsup +! use coreprofiles, only : psdbnd +! implicit none +! ! arguments +! real(wp_), intent(in) :: rrm,zzm +! ! local variables +! real(wp_) :: psinv +! +! if(iequil.eq.1) then +! call equian(rrm,zzm,psinv) +! else +! call equinum_psi(rrm,zzm,psinv) +! end if +! +! inside_plasma = (psinv >= zero .and. psinv < psdbnd) .and. & +! (psinv >= one .or. (zzm >= zbinf .and. zzm <= zbsup)) +! end function inside_plasma +! +! +! +! subroutine vacuum_rt(xv0,anv0,xvend,dstvac,ivac) +! use const_and_precisions, only : wp_ +! use beamdata, only : dst +! use limiter, only : rlim,zlim,nlim +! implicit none +! ! arguments +! real(wp_), dimension(3), intent(in) :: xv0,anv0 +! real(wp_), dimension(3), intent(out) :: xvend +! real(wp_), intent(out) :: dstvac +! integer, intent(out) :: ivac +! ! local variables +! integer :: i +! real(wp_) :: st,rrm,zzm,smax +! real(wp_), dimension(3) :: walln +! logical :: plfound +! +! ! ivac=1 plasma hit before wall reflection +! ! ivac=2 wall hit before plasma +! ! ivac=-1 vessel (and thus plasma) never crossed +! +! call inters_linewall(xv0/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), & +! nlim,smax,walln) +! smax=smax*1.0e2_wp_ +! rrm=1.0e-2_wp_*sqrt(xv0(1)**2+xv0(2)**2) +! zzm=1.0e-2_wp_*xv0(3) +! if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then +! ! first wall interface is outside-inside +! if (dot_product(walln,walln) h .or. & + ah > 0.0_wp_ .and. a(jm) <= h) then + ix=ix+1 + lx(ix)=-j + end if + if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & + ah > 0.0_wp_ .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + end do + end do + + do jm=nr,mxr,nrqmax + j = jm + nrqmax + ah=a(j)-h + if (ah <= 0.0_wp_ .and. a(j-1) > h .or. & + ah > 0.0_wp_ .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + if (ah <= 0.0_wp_ .and. a(jm) > h .or. & + ah > 0.0_wp_ .and. a(jm) <= h) then + ix=ix+1 + lx(ix)=-j + end if + end do + + do jm=1,mxr,nrqmax + j = jm + nrqmax + if (a(j) <= h .and. a(jm) > h .or. & + a(j) > h .and. a(jm) <= h) then + ix=ix+1 + lx(ix) =-j + end if + end do + + do j=2,nr + if (a(j) <= h .and. a(j-1) > h .or. & + a(j) > h .and. a(j-1) <= h) then + ix=ix+1 + lx(ix)=j + end if + end do + + if(ix<=0) return + +bb: do + in=ix + jx=lx(in) + jfor=0 + lda=1 + ldb=2 + + do + if(jx<0) then + jabs=-jx + jnb = jabs - nrqmax + else + jabs=jx + jnb=jabs-1 + end if + + adn=a(jabs)-a(jnb) + if(adn/=0) px=(a(jabs)-h)/adn + kx = (jabs - 1) / nrqmax + ikx = jabs - nrqmax * kx - 1 + + if(jx<0) then + x = drgrd * ikx + y = dzgrd * (kx - px) + else + x = drgrd * (ikx - px) + y = dzgrd * kx + end if + + icount = icount + 1 + rcon(icount) = x + rqgrid(1) + zcon(icount) = y + zqgrid(1) + mpl= icount + itm = 1 + ja(1,1) = jabs + nrqmax + j=1 + + if(jx<=0) then + ja(1,1) = -jabs-1 + j=2 + end if + + ja(2,1) = -ja(1,1) + ja(3,1) = -jx + 1 - nrqmax + ja(3,2) = -jx + ja(j,2) = jabs - nrqmax + k= 3-j + ja(k,2) = 1-jabs + + if (kx<=0 .or. ikx<=0) then + lda=1 + ldb=lda + else if (ikx + 1 - nr >= 0 .and. jx <= 0) then + lda=2 + ldb=lda + else if(jfor/=0) then + lda=2 + do i=1,3 + if(jfor==ja(i,2)) then + lda=1 + exit + end if + end do + ldb=lda + end if + + flag1=.false. + aa: do k=1,3 + do l=lda,ldb + do i=1,ix + if(lx(i)==ja(k,l)) then + itm=itm+1 + inext= i + if(jfor/=0) exit aa + if(itm .gt. 3) then + flag1=.true. + exit aa + end if + end if + end do + end do + end do aa + + if(.not.flag1) then + lx(in)=0 + if(itm .eq. 1) exit + end if + + jfor=jx + jx=lx(inext) + in = inext + end do + + do + if(lx(ix)/=0) then + if(mpl>=4) then + ncon = ncon + 1 + npts(ncon) = icount - iclast + iclast = icount + end if + exit + end if + ix= ix-1 + if(ix<=0) exit bb + end do + + end do bb + + if(mpl >= 4) then + ncon = ncon + 1 + npts(ncon) = icount - iclast + iclast = icount + end if + end subroutine cniteq + + + + subroutine print_headers + use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm + implicit none + write(uprj0,*) ' #sst j k xt yt zt rt' + write(uprj0+1,*) ' #sst j k xt yt zt rt' + write(uwbm,*) ' #sst w1 w2' + write(udisp,*) ' #sst Dr_Nr Di_Nr' + write(ucenr,*) ' #sst R z phi psin rhot ne Te Btot Nperp Npl ki '// & + 'alpha tau Pt dIds nhmax iohkw index_rt ddr' + write(uoutr,*) ' #i k sst x y R z psin tau Npl alpha index_rt' + write(upec,*) ' #rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt' + write(usumm,*) ' #Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // & + 'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt ' // & + 'Jphimx dPdVmx drhotj drhotp' + end subroutine print_headers + + + + subroutine print_prof + use const_and_precisions, only : wp_ + use equilibrium, only : psinr,nq,fq,frhotor,tor_curr_psi + use coreprofiles, only : density, temp + use units, only : uprfin + implicit none +! local constants + real(wp_), parameter :: eps=1.e-4_wp_ +! local variables + integer :: i + real(wp_) :: psin,rhot,ajphi,dens,ddens + + write(uprfin,*) ' #psi rhot ne Te q Jphi' + do i=1,nq + psin=psinr(i) + rhot=frhotor(sqrt(psin)) + call density(psin,dens,ddens) + call tor_curr_psi(max(eps,psin),ajphi) + + write(uprfin,"(12(1x,e12.5))") psin,rhot,dens,temp(psin),fq(psin),ajphi*1.e-6_wp_ + end do + end subroutine print_prof + + + + subroutine print_bres(bres) + use const_and_precisions, only : wp_ + use gray_params, only : iequil + use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, nq + use units, only : ubres + implicit none +! arguments + real(wp_) :: bres +! local constants + integer, parameter :: icmx=2002 +! local variables + integer :: j,k,n,nconts,nctot + integer, dimension(10) :: ncpts + real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb + real(wp_), dimension(icmx) :: rrcb,zzcb + real(wp_) :: rv(nq), zv(nq) + real(wp_), dimension(nq,nq) :: btotal + + + dr = (rmxm-rmnm)/(nq-1) + dz = (zmxm-zmnm)/(nq-1) + do j=1,nq + rv(j) = rmnm + dr*(j-1) + zv(j) = zmnm + dz*(j-1) + end do + +! Btotal on psi grid + btmx=-1.0e30_wp_ + btmn=1.0e30_wp_ + do k=1,nq + zzk=zv(k) + do j=1,nq + rrj=rv(j) + call bfield(rrj,zzk,bbphi,bbr,bbz) + btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2) + if(btotal(j,k).ge.btmx) btmx=btotal(j,k) + if(btotal(j,k).le.btmn) btmn=btotal(j,k) + enddo + enddo + +! compute Btot=Bres/n with n=1,5 + write(ubres,*)'#i Btot R z' + 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,nq,nq,bbb,nconts,ncpts,nctot,rrcb,zzcb) + do j=1,nctot + write(ubres,'(i6,12(1x,e12.5))') j,bbb,rrcb(j),zzcb(j) + end do + end if + write(ubres,*) + end do + end subroutine print_bres + + + + subroutine print_surfq(qval) + use const_and_precisions, only : wp_, one + use equilibrium, only : psinr,nq,fq,frhotor,rmaxis,zmaxis, & + zbsup,zbinf + use magsurf_data, only : contours_psi,npoints,print_contour + use utils, only : locate, intlin + implicit none +! arguments + real(wp_), dimension(:), intent(in) :: qval +! local variables + integer :: ncnt,i1,i + real(wp_) :: rup,zup,rlw,zlw,rhot,psival + real(wp_), dimension(npoints) :: rcn,zcn + real(wp_), dimension(nq) :: qpsi + +! build q profile on psin grid + do i=1,nq + qpsi(i) = fq(psinr(i)) + end do +! locate psi surface for q=qval + print* + do i=1,size(qval) + call locate(abs(qpsi),nq,qval(i),i1) !!!! check for non monotonous q profile + if (i1>0.and.i1=rtimx .and. jkv(1)==nrayr) rtimx = rti + if(rti<=rtimn .and. jkv(1)==nrayr) rtimn = rti + end do + write(uprj,*) + write(uwbm,'(3(1x,e16.8e3))') st,rtimn,rtimx + end subroutine print_projxyzt + + + subroutine print_output(i,jk,st,qj,xv,psinv,btot,ak0,anpl,anpr,anprim, & 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 - use beamdata, only : nray,nrayth + use beamdata, only : nray,nrayth,jkray1 + use units, only : ucenr,uoutr,udisp implicit none ! arguments integer, intent(in) :: i,jk,nhf,iokhawa,index_rt @@ -1347,8 +1866,7 @@ contains zzm=xv(3)*1.0e-2_wp_ rrm=sqrt(xxm**2 + yym**2) -! central ray only begin -! print dIds in A/m/W, ki in m^-1 +! print central ray trajectory. dIds in A/m/W, ki in m^-1 if(jk.eq.1) then phideg=atan2(yym,xxm)/degree if(psinv>=zero .and. psinv<=one) then @@ -1360,23 +1878,59 @@ contains pt=exp(-tau) 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,didsn,dble(nhf),dble(iokhawa), & - dble(index_rt),ddr + write(ucenr,'(16(1x,e16.8e3),3i5,1x,e16.8e3)') stm,rrm,zzm,phideg, & + psinv,rhot,dens,tekev,btot,anpr,anpl,akim,alpha,tau,pt,didsn, & + nhf,iokhawa,index_rt,ddr end if -! central ray only end ! print conservation of dispersion relation - if(jk==nray) write(17,'(30(1x,e16.8e3))') st,ddr,ddi + if(jk==nray) write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi ! print outer trajectories if(mod(i,istpl0)==0) then - k = jk + nrayth - nray - if(k>0) then - write(33,'(2i5,16(1x,e16.8e3))') i,k,stm,xxm,yym,rrm,zzm, & - psinv,tau,anpl,alpha,dble(index_rt) + k = jk - jkray1 + 1 + if(k>0 .and. k<=nrayth) then + write(uoutr,'(2i5,9(1x,e16.8e3),i5)') i,k,stm,xxm,yym,rrm,zzm, & + psinv,tau,anpl,alpha,index_rt end if end if end subroutine print_output + + + subroutine print_pec(rhop_tab,rhot_tab,jphi,jcd,dpdv,currins,pins,index_rt) + use const_and_precisions, only : wp_ + use units, only : upec + implicit none +! arguments + real(wp_), dimension(:), intent(in) :: rhop_tab,rhot_tab,jphi,jcd,dpdv, & + currins,pins + integer, intent(in) :: index_rt +! local variables + integer :: i + + do i=1,size(rhop_tab) + write(upec,'(7(1x,e16.8e3),i5)') rhop_tab(i),rhot_tab(i), & + jphi(i),jcd(i),dpdv(i),currins(i),pins(i),index_rt + end do + end subroutine print_pec + + + + subroutine print_finals(pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, & + drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, & + stmx,psipol,chipol,index_rt) + use const_and_precisions, only : wp_ + use units, only : usumm + implicit none + real(wp_), intent(in) :: pabs,icd,dpdvp,jphip,rhotpav,rhotjava,drhotpav, & + drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx,ratjbmx, & + stmx,psipol,chipol + integer, intent(in) :: index_rt + + write(usumm,'(15(1x,e12.5),i5,4(1x,e12.5))') icd,pabs,jphip,dpdvp, & + rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, & + stmx,psipol,chipol,index_rt,jphimx,dpdvmx,drhotj,drhotp + end subroutine print_finals + end module graycore diff --git a/src/limiter.f90 b/src/limiter.f90 new file mode 100644 index 0000000..6e5f05c --- /dev/null +++ b/src/limiter.f90 @@ -0,0 +1,33 @@ +module limiter + use const_and_precisions, only : wp_ + ! === 1D array limiter Rlim_i, Zlim_i + integer, public, save :: nlim + real(wp_), save :: rwallm + real(wp_), dimension(:), allocatable, save :: rlim,zlim + +contains + + 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 + + + + subroutine unset_lim + use const_and_precisions, only : zero + implicit none + if(allocated(rlim)) deallocate(rlim) + if(allocated(zlim)) deallocate(zlim) + nlim=0 + rwallm=zero + end subroutine unset_lim + +end module limiter \ No newline at end of file diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index eaf80b5..ea10095 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -23,7 +23,7 @@ module magsurf_data contains - subroutine alloc_surf_anal(ierr) + subroutine alloc_cnt(ierr) implicit none integer, intent(out) :: ierr @@ -32,18 +32,16 @@ contains return end if - call dealloc_surf_anal - allocate(psicon(npsi),rcon(npsi,npoints), & - zcon(npsi,npoints),stat=ierr) - if (ierr/=0) call dealloc_surf_anal - end subroutine alloc_surf_anal + call dealloc_cnt + allocate(psicon(npsi),rcon(npoints,npsi),zcon(npoints,npsi)) + end subroutine alloc_cnt - subroutine dealloc_surf_anal + subroutine dealloc_cnt implicit none if(allocated(psicon)) deallocate(psicon) if(allocated(rcon)) deallocate(rcon) if(allocated(zcon)) deallocate(zcon) - end subroutine dealloc_surf_anal + end subroutine dealloc_cnt subroutine alloc_surfvec(ierr) @@ -56,21 +54,19 @@ contains end if call dealloc_surfvec - allocate(psicon(npsi),rcon(npsi,npoints),zcon(npsi,npoints),pstab(npsi), & + call alloc_cnt(ierr) + allocate(pstab(npsi), & rhot_eq(npsi),rhotqv(npsi),bav(npsi),bmxpsi(npsi),bmnpsi(npsi),varea(npsi), & vvol(npsi),vcurrp(npsi),vajphiav(npsi),qqv(npsi),ffc(npsi),vratja(npsi), & vratjb(npsi),rpstab(npsi),rri(npsi),rbav(npsi),cdadrhot(npsi,4), & cdvdrhot(npsi,4),cbmx(npsi,4),cbmn(npsi,4),crbav(npsi,4),cvol(npsi,4), & crri(npsi,4),carea(npsi,4),cfc(npsi,4),crhotq(npsi,4),cratjpl(npsi,4), & - cratja(npsi,4),cratjb(npsi,4),stat=ierr) - if (ierr/=0) call dealloc_surf_anal + cratja(npsi,4),cratjb(npsi,4)) end subroutine alloc_surfvec subroutine dealloc_surfvec implicit none - if(allocated(psicon)) deallocate(psicon) - if(allocated(rcon)) deallocate(rcon) - if(allocated(zcon)) deallocate(zcon) + call dealloc_cnt if(allocated(pstab)) deallocate(pstab) if(allocated(rhot_eq)) deallocate(rhot_eq) if(allocated(rhotqv)) deallocate(rhotqv) @@ -104,99 +100,6 @@ contains end subroutine dealloc_surfvec - - subroutine contours_psi(h,rup,zup,rlw,zlw,rcn,zcn,ipr) - use const_and_precisions, only : wp_,pi - use equilibrium, only : psiant,psinop,nsr,nsz,cc=>cceq,tr,tz,kspl, & - points_tgo - use dierckx, only : profil,sproota - use reflections, only : rwallm - implicit none -! local constants - integer, parameter :: mest=4 -! arguments - integer, intent(in) :: ipr - real(wp_), intent(in) :: h - real(wp_), intent(inout) :: rup,zup,rlw,zlw - real(wp_), dimension(npoints), intent(out) :: rcn,zcn -! local variables - integer :: np,info,ic,ier,ii,iopt,m - real(wp_) :: ra,rb,za,zb,th,zc,val - real(wp_), dimension(mest) :: zeroc - real(wp_), dimension(nsr) :: czc - - np=(npoints-1)/2 - - ra=rup - rb=rlw - za=zup - zb=zlw - call points_tgo(ra,za,rup,zup,h,info) - call points_tgo(rb,zb,rlw,zlw,h,info) - - th=pi/dble(np) - rcn(1)=rlw - zcn(1)=zlw - rcn(npoints)=rlw - zcn(npoints)=zlw - rcn(np+1)=rup - zcn(np+1)=zup - do ic=2,np - zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_ - iopt=1 - call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) - if(ier.gt.0) print*,' profil =',ier - val=h*psiant+psinop - call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) - if (zeroc(1).gt.rwallm) then - rcn(ic)=zeroc(1) - rcn(npoints+1-ic)=zeroc(2) - else - rcn(ic)=zeroc(2) - rcn(npoints+1-ic)=zeroc(3) - end if - zcn(ic)=zc - zcn(npoints+1-ic)=zc - end do - - if (ipr.gt.0) then - do ii=1,npoints - write(71,'(i6,12(1x,e12.5))') ii,h,rcn(ii),zcn(ii) - end do - write(71,*) - write(71,*) - end if - end subroutine contours_psi - - - - subroutine contours_psi_an(h,rcn,zcn,ipr) - use const_and_precisions, only : wp_,pi - use equilibrium, only : frhotor,aminor,rmaxis,zmaxis - implicit none -! arguments - integer :: ipr - real(wp_) :: h - real(wp_), dimension(npoints) :: rcn,zcn -! local variables - integer :: np,ic - real(wp_) :: rn,th - - np=(npoints-1)/2 - th=pi/dble(np) - rn=frhotor(sqrt(h)) - - do ic=1,npoints - zcn(ic)=zmaxis+aminor*rn*sin(th*(ic-1)) - rcn(ic)=rmaxis+aminor*rn*cos(th*(ic-1)) - - if (ipr.gt.0) write(71,'(i6,12(1x,e12.5))') ic,h,rcn(ic),zcn(ic) - end do - if (ipr.gt.0) write(71,*) - end subroutine contours_psi_an - - - subroutine flux_average use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use gray_params, only : iequil @@ -212,7 +115,7 @@ 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 + integer :: ier,ierr,l,jp,ipr,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,r2iav,currp, & @@ -232,7 +135,6 @@ contains real(wp_) :: fpolv,ddpsidrr,ddpsidzz npsi=nnintp - ninpr=(npsi-1)/10 npoints = 2*ncnt+1 call alloc_surfvec(ierr) @@ -245,8 +147,6 @@ contains ! computation of flux surface averaged quantities - write(71,*)' #i psin R z' - dlam=1.0_wp_/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam @@ -280,8 +180,8 @@ contains ajphiav=-ccj*(ddpsidrr+ddpsidzz)/rmaxis psicon(1)=0.0_wp_ - rcon(1,:)=rmaxis - zcon(1,:)=zmaxis + rcon(:,1)=rmaxis + zcon(:,1)=zmaxis pstab(1)=0.0_wp_ rpstab(1)=0.0_wp_ vcurrp(1)=0.0_wp_ @@ -309,21 +209,14 @@ contains do jp=2,npsi height=dble(jp-1)/dble(npsi-1) if(jp.eq.npsi) height=0.9999_wp_ - ipr=0 - jpr=mod(jp,ninpr) - if(jpr.eq.1) ipr=1 rhopjp=height psinjp=height*height rhotjp=frhotor(rhopjp) psicon(jp)=height - if(iequil<2) then - call contours_psi_an(psinjp,rctemp,zctemp,ipr) - else - call contours_psi(psinjp,rup,zup,rlw,zlw,rctemp,zctemp,ipr) - end if - rcon(jp,:) = rctemp - zcon(jp,:) = zctemp + call contours_psi(psinjp,rctemp,zctemp,rup,zup,rlw,zlw) + rcon(:,jp) = rctemp + zcon(:,jp) = zctemp r2iav=0.0_wp_ anorm=0.0_wp_ @@ -390,7 +283,7 @@ contains ! computation maximum/minimum B values on given flux surface if(btot.le.bmmn) bmmn=btot if(btot.ge.bmmx) bmmx=btot - end do + end do ! bav= [T] , b2av= [T^2] , rbav=/b_min ! anorm = int d l_p/B_p = dV/dpsi/(2pi) @@ -429,8 +322,8 @@ contains vratjb(jp)=ratio_cdbtor qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) qqv(jp)=qq - dadrhotv(jp)=phitedge*frhotor(rhopjp)/fq(psinjp)*dadpsi/pi - dvdrhotv(jp)=phitedge*frhotor(rhopjp)/fq(psinjp)*dvdpsi/pi + dadrhotv(jp)=phitedge*rhotjp/fq(psinjp)*dadpsi/pi + dvdrhotv(jp)=phitedge*rhotjp/fq(psinjp)*dvdpsi/pi ! computation of fraction of circulating/trapped fraction fc, ft ! and of function H(lambda,rhop) @@ -474,19 +367,8 @@ contains dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l) end do end do - - write(56,*)' #rhop rhot || |Bmx| |Bmn| Area Vol |I_pl| fc ratJa ratJb' - - do jp=1,npsi - if(jp.eq.npsi) then - rpstab(jp)=1.0_wp_ - pstab(jp)=1.0_wp_ - end if - rhotjp=frhotor(rpstab(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 + rpstab(npsi)=1.0_wp_ + pstab(npsi)=1.0_wp_ ! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs ! used for computations of dP/dV and J_cd @@ -525,9 +407,19 @@ contains wrk,lwrk,iwrk,kwrk,ier) njpt=njp nlmt=nlm + + + do jp=1,npsi + call print_fluxav(rpstab(jp),frhotor(rpstab(jp)),bav(jp),bmxpsi(jp), & + bmnpsi(jp),varea(jp),vvol(jp),vcurrp(jp),vajphiav(jp), & + ffc(jp),vratja(jp),vratjb(jp)) + end do + + ninpr=(npsi-1)/10 + do jp=ninpr+1,npsi,ninpr + call print_contour(psicon(jp),rcon(:,jp),zcon(:,jp)) + end do - 99 format(20(1x,e12.5)) - end subroutine flux_average @@ -569,4 +461,115 @@ contains end subroutine fluxval + + + subroutine contours_psi(h,rcn,zcn,rup,zup,rlw,zlw) + use const_and_precisions, only : wp_,pi + use gray_params, only : iequil + use dierckx, only : profil,sproota + use equilibrium, only : rmaxis,zmaxis,aminor,frhotor,tr,nsr,tz,nsz,cceq, & + kspl,psiant,psinop,points_tgo + use limiter, only : rwallm + implicit none +! local constants + integer, parameter :: mest=4 +! arguments + real(wp_), intent(in) :: h + real(wp_), dimension(:), intent(out) :: rcn,zcn + real(wp_), intent(inout) :: rup,zup,rlw,zlw +! local variables + integer :: npoints,np,info,ic,ier,ii,iopt,m + real(wp_) :: ra,rb,za,zb,rn,th,zc,val + real(wp_), dimension(mest) :: zeroc + real(wp_), dimension(nsr) :: czc + + npoints=size(rcn) + np=(npoints-1)/2 + th=pi/dble(np) + + if (iequil<2) then + + rn=frhotor(sqrt(h)) + do ic=1,npoints + zcn(ic)=zmaxis+aminor*rn*sin(th*(ic-1)) + rcn(ic)=rmaxis+aminor*rn*cos(th*(ic-1)) + end do + + else + + ra=rup + rb=rlw + za=zup + zb=zlw + call points_tgo(ra,za,rup,zup,h,info) + call points_tgo(rb,zb,rlw,zlw,h,info) + + rcn(1)=rlw + zcn(1)=zlw + rcn(npoints)=rlw + zcn(npoints)=zlw + rcn(np+1)=rup + zcn(np+1)=zup + do ic=2,np + zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_ + iopt=1 + call profil(iopt,tr,nsr,tz,nsz,cceq,kspl,kspl,zc,nsr,czc,ier) + if(ier.gt.0) print*,' profil =',ier + val=h*psiant+psinop + call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) + if (zeroc(1).gt.rwallm) then + rcn(ic)=zeroc(1) + rcn(npoints+1-ic)=zeroc(2) + else + rcn(ic)=zeroc(2) + rcn(npoints+1-ic)=zeroc(3) + end if + zcn(ic)=zc + zcn(npoints+1-ic)=zc + end do + + end if + end subroutine contours_psi + + + + subroutine print_contour(psin,rc,zc) + use const_and_precisions, only : wp_, comp_tiny + use units, only : ucnt + implicit none +! arguments + real(wp_), intent(in) :: psin + real(wp_), dimension(:), intent(in) :: rc, zc +! local variables + integer :: npoints,ic + + if (psin < comp_tiny) then + write(ucnt,*)' #i psin R z' + else + npoints=size(rc) + do ic=1,npoints + write(ucnt,'(i6,12(1x,e12.5))') ic,psin,rc(ic),zc(ic) + end do + write(ucnt,*) + write(ucnt,*) + end if + end subroutine print_contour + + + + subroutine print_fluxav(psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, & + ffc,ratja,ratjb) + use const_and_precisions, only : wp_, comp_tiny + use units, only : uflx + implicit none +! arguments + real(wp_), intent(in) :: psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, & + ffc,ratja,ratjb + + if (psin < comp_tiny) & + write(uflx,*)' #rhop rhot || |Bmx| |Bmn| Area Vol |I_pl| fc ratJa ratJb' + write(uflx,'(20(1x,e12.5))') psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, & + ffc,ratja,ratjb + end subroutine print_fluxav + end module magsurf_data diff --git a/src/main.f90 b/src/main.f90 index 0b857ba..ccfc85a 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,6 +1,6 @@ -program gray_main +program main_std use const_and_precisions, only : wp_,one - use graycore, only : gray + use graycore, only : gray_main use gray_params, only : read_params, antctrl_type,eqparam_type, & prfparam_type,outparam_type,rtrparam_type,hcdparam_type use beams, only : read_beam0, read_beam1, read_beam2 @@ -48,8 +48,9 @@ program gray_main 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)) +! ??? analytical only? change for numerical! +! 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) @@ -108,8 +109,8 @@ program gray_main ! ========================= 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, & + call gray_main(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) ! ======================================================================== @@ -129,4 +130,4 @@ program gray_main 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 +end program main_std \ No newline at end of file diff --git a/src/pec.f90 b/src/pec.f90 index cee1a67..05c0482 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -232,8 +232,9 @@ contains end subroutine pec_tab - subroutine postproc_profiles(pabs,currt,rhot_tab,dpdv,ajphiv, & - rhotpav,drhotpav,rhotjava,drhotjava) + subroutine postproc_profiles(pabs,currt,rhot_tab,dpdv,ajphiv, & + rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, & + rhotp, drhotp, rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx) ! radial average values over power and current density profile use const_and_precisions, only : pi use gray_params, only : nnd @@ -243,14 +244,14 @@ contains real(wp_), intent(in) :: pabs,currt real(wp_), dimension(nnd), intent(in) :: rhot_tab real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv - real(wp_), intent(out) :: rhotpav,rhotjava - real(wp_), intent(out) :: drhotpav,drhotjava - real(wp_) :: rhopjava,rhoppav - real(wp_) :: dpdvp,dpdvmx,rhotp,drhotp - real(wp_) :: ajphip,ajmxfi,rhotjfi,drhotjfi - real(wp_) :: ratjamx,ratjbmx,ratjplmx + real(wp_), intent(out) :: rhotpav,drhotpav,dpdvp + real(wp_), intent(out) :: rhotjava,drhotjava,ajphip + + real(wp_), intent(out) :: rhotp,drhotp,dpdvmx + real(wp_), intent(out) :: rhotjfi,drhotjfi,ajmxfi + real(wp_), intent(out) :: ratjamx,ratjbmx - real(wp_) :: sccsa + real(wp_) :: sccsa,ratjplmx,rhopjava,rhoppav real(wp_) :: rhotjav,rhot2pav,rhot2java,dvdrhotav,dadrhotava rhotpav=zero diff --git a/src/reflections.f90 b/src/reflections.f90 index 1a3391e..2409110 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -2,15 +2,10 @@ module reflections use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge, zero, one implicit none - ! === 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,range2rect,set_lim + public :: wall_refl,range2rect contains @@ -29,6 +24,8 @@ subroutine reflect(ki,nsurf,ko) end if end subroutine reflect + + subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) implicit none real(wp_), intent(in), dimension(3) :: xv,kv @@ -90,6 +87,8 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) if (dot_product(normw,kv)>zero) normw=-normw end subroutine inters_linewall + + subroutine linecone_coord(xv,kv,rs,zs,s,t,n) use utils, only : bubble implicit none @@ -161,6 +160,8 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n) end if end subroutine linecone_coord + + subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr) implicit none real(wp_), dimension(2), intent(in) :: xa,ya,xb,yb @@ -183,6 +184,8 @@ subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr) end if end subroutine interssegm_coord + + function interssegm(xa,ya,xb,yb) implicit none real(wp_), dimension(2), intent(in) :: xa,ya,xb,yb @@ -195,6 +198,21 @@ function interssegm(xa,ya,xb,yb) t>=zero .and. t<=one) interssegm = .true. end function interssegm + + + 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 + + + function inside(xc,yc,n,x,y) use utils, only : locatef, locate_unord, intlinf, bubble implicit none @@ -221,28 +239,10 @@ function inside(xc,yc,n,x,y) inside=(mod(locatef(xint,nj,x),2)==1) end function inside -subroutine alloc_lim(ier) - implicit none - integer, intent(out) :: ier - - if(nlim.lt.0) then - ier = -1 - return - end if - - call dealloc_lim - allocate(rlim(nlim),zlim(nlim), & - stat=ier) - if (ier/=0) call dealloc_lim -end subroutine alloc_lim -subroutine dealloc_lim - implicit none - if(allocated(rlim)) deallocate(rlim) - if(allocated(zlim)) deallocate(zlim) -end subroutine dealloc_lim subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl) + use limiter, only : rlim,zlim,nlim implicit none ! arguments integer :: irfl @@ -316,28 +316,5 @@ subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl) eztr=dot_product(vv3,evrfl) 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 diff --git a/src/units.f90 b/src/units.f90 new file mode 100644 index 0000000..d2db0a0 --- /dev/null +++ b/src/units.f90 @@ -0,0 +1,13 @@ +module units +! STANDARD + integer, parameter :: uprm = 2, ubeam = 97, uprf = 98, ueq = 99 + integer, parameter :: uprfin = 55, uflx = 56, ubres = 70, ucnt = 71 + integer, parameter :: uprj0 = 8, uwbm = 12, ucenr = 4, uoutr = 33 + integer, parameter :: udisp = 17, upec = 48, usumm = 7 +! JINTRAC +! integer, parameter :: uprm =602, ubeam =603, uprf =644, ueq =644 +! integer, parameter :: uprfin =645, uflx =646, ubres =630, ucnt =631 +! integer, parameter :: uprj0 =608, uwbm =612, ucenr =604, uoutr =633 +! integer, parameter :: udisp =617, upec =648, usumm =607, uhead =638 + +end module units \ No newline at end of file diff --git a/srcjetto/gray.f b/srcjetto/gray.f new file mode 100644 index 0000000..84cce7a --- /dev/null +++ b/srcjetto/gray.f @@ -0,0 +1,152 @@ +! Fortran 77 interface to JETTO + subroutine gray(ijetto, mr, mz, mrd, r, z, psin, psiax, psibnd, + . rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, + . qsf, nbeam, powin, alphin, betain, dpdv, jcd, pec, icd, ierr) + use units, only : ucenr,usumm,uprj0,uprj0+1,uwbm,udisp,ubres, + . ucnt,uoutr,ueq,uprfin,uflx,upec + use gray_params, only : read_params + use beams, only : read_beam2 + use graycore, only : gray_main + implicit none +! input arguments + integer ijetto, mr, mz, nbnd, nrho, nbeam + real*8 r(mr), z(mz), psin(mrd,mz) + real*8 psiax, psibnd, rax, zax + real*8 rbnd(nbnd), zbnd(nbnd) + real*8 psijet(nrho), f(nrho), qsf(nrho), te(nrho), dne(nrho) + real*8 zeff(nrho) + real*8 powin(nbeam), alphin(nbeam), betain(nbeam) +! output arguments + real*8 dpdv(nrho), jcd(nrho), pec, icd +! gray_main output arguments + real*8 dpdvloop(nrho), jcdloop(nrho), pecloop, icdloop + integer ierr +! local variables + real*8 rlim(5),zlim(5) + logical firstcall=.true. + save firstcall + +! === input arguments ================================================== +! +! ijetto Equilibrium source (1 EFIT, 2 ESCO) +! If IJETTO=2, then PSIN values are valid only inside +! plasma boudary (PSIN=0 outside) +! mr Size of flux map grid in R direction +! mz Size of flux map grid in Z direction +! mrd Leading dimension of the psin(:,:) array, mrd>mr +! r R coordinates of flux map grid points [m] +! z Z coordinates of flux map grid points [m] +! psin Normalised poloidal flux psin=(psi-psiax)/(psibnd-psiax) +! on the (R, Z) grid. +! psiax Poloidal flux on axis [Wb rad-1] +! psibnd Poloidal flux on boundary [Wb rad-1] +! rax R coordinate of magnetic axis [m] +! zax Z coordinate of magnetic axis [m] +! nbnd Number of points in plasma boundary contour +! rbnd R coordinates of plasma boundary contour [m] +! zbnd Z coordinates of plasma boundary contour [m] +! +! nrho Number of points in JETTO rho grid - +! psijet Normalised poloidal flux on JETTO radial grid +! f Poloidal current stream function f=Bphi*R on JETTO +! radial grid [T m] +! te Electron temperature on JETTO radial grid [eV] +! dne Electron density on JETTO radial grid [m-3] +! zeff Effective nuclear charge Zeff on JETTO radial grid +! qsf Safety factor on JETTO radial grid +! +! nbeam Total number of injected beams +! powin Input ECRH power array [W] (powin(i) =< 0 means i-th beam is unused) +! alphin Beams poloidal injection angles array [rad] +! betain Beams toroidal injection angles array [rad] +! +! === output arguments ================================================= +! +! dpdv Absorbed EC power density on JETTO radial grid [W m-3] +! jcd EC driven flux averaged current density on JETTO +! radial grid [A m-2] +! pec Total absorbed EC power [W] +! icd Total EC driven current [A] +! ierr Return code. IERR>0 on error +! ierr = 90-93: error computing integrals for current drive +! ierr = 94: absorption coefficient alpha < 0 +! ierr = 97: parallel comp. refract. idx N//>0.99 (warning) +! ierr = 98: parallel comp. refract. idx N//>1.05 +! +! === Note ============================================================= +! +! JETTO coordinate system assumes toroidal angle increasing CW +! in GRAY toroidal angle increases CCW --> adapt signs on input data +! +! f is passed as -f +! qsf is passed as -qsf +! +! jcd is returned as -jcd +! icd is returned as -icd +! +! ====================================================================== + +! if first call read parameters from external file + if (firstcall) then + call read_params('gray.data',rtrp,hcdp,antp,eqp,rwallm, + . prfp,outp,uprm) + antp%filenm='graybeam.data' + eqp%filenm='JETTO' + eqp%iequil=ijetto+1 + prfp%filenm='JETTO' + firstcall=.false. + end if + +! set output variables to 0 + do i=1,nrho + dpdv(i) = 0.d0 + jcd(i) = 0.d0 + end do + pec = 0.d0 + icd = 0.d0 + +! loop over beams with power>0 + do j=1,nbeam + if (powin(j).gt.0.0d0) cycle + +! read j-th beam properties from file +! and adjust alpha/beta if out of the allowed range + alpha0=alphin(j) + beta0=betain(j) + p0mw=powin(j)*1.d-6 + call read_beam2(antp%filenm,j,alpha0,beta0,fghz,antp%iox, + . x0,y0,z0,w1,w2,ri1,ri2,phiw,phir,ubeam) + psipol0=antp%psi + chipol0=antp%chi + iox0=antp%iox + +! set simple limiter + r0m=sqrt(x0**2+y0**2)*0.01d0 + call range2rect(rwallm,max(r0m,rv(mr)),zv(1),zv(mz),rlim,zlim) + +! call main subroutine for the j-th beam + subroutine gray_main(r,z,psin(1:mr,:),psibnd-psiax, + . psijet,-f,-qsf,rax,rax,zax,rbnd,zbnd,eqp, + . psijet,te,dne,zeff,prfp,rlim,zlim, + . p0mw,fghz,alpha0,beta0,(/x0,y0,z0/), + . w1,w2,ri1,ri2,phiw,phir,iox0,psipol0,chipol0, + . dpdvloop,jcdloop,pecpool,icdloop,outp,rtrp,hcdp,ierr) + +! add contribution of j-th beam to the total +! adapting output data to JETTO convention on toroidal angle + do i=1,nrho + dpdv(i) = dpdv(i) + dpdvloop(i) + jcd(i) = jcd(i) - jcdloop(i) + end do + pec = pec + pecloop + icd = icd - icdloop + +! end of loop over beams with power>0 + end do + +! close output (debug) files + close(ucenr,usumm,uprj0,uprj0+1,uwbm,udisp,ubres,ucnt,uoutr, + . ueq,uprfin,uflx,upec) + + return + end