From 693af2a763335da6ddd8604712afadecd05538e1 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Wed, 15 Dec 2021 02:30:55 +0100 Subject: [PATCH] merge graysum and gray Note that graysum, now as sum_mode, still doesn't support multipass. --- Makefile | 6 +- src/graycore.f90 | 165 ++++ src/main-sum.f90 | 195 ----- src/main.f90 | 246 +++--- src/sumcore.f90 | 1907 ---------------------------------------------- 5 files changed, 327 insertions(+), 2192 deletions(-) delete mode 100644 src/main-sum.f90 delete mode 100644 src/sumcore.f90 diff --git a/Makefile b/Makefile index 9bc63a4..9f9849c 100644 --- a/Makefile +++ b/Makefile @@ -38,7 +38,7 @@ MODULES = $(filter-out $(OBJDIR)/main%,$(OBJECTS)) # Build outputs GRAY = $(BINDIR)/gray LIBGRAY = $(LIBDIR)/libgray -BINARIES = $(GRAY) $(GRAY)sum +BINARIES = $(GRAY) LIBRARIES = $(LIBGRAY).so ## @@ -89,10 +89,6 @@ install: $(BINARIES) $(LIBRARIES) doc/gray.1 $(GRAY): $(shell ./depend $(OBJDIR)/main.o) | $(BINDIR) $(LD) $(LDFLAGS) -o $@ $^ -# sum variant -$(GRAY)sum: $(shell ./depend $(OBJDIR)/main-sum.o) | $(BINDIR) - $(LD) $(LDFLAGS) -o $@ $^ - # GRAY shared library $(LIBGRAY).so: $(MODULES) | $(LIBDIR) $(LD) -shared $(LDFLAGS) -o $@ $^ diff --git a/src/graycore.f90 b/src/graycore.f90 index 5348f6f..77a9d5b 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -603,6 +603,171 @@ contains + subroutine sum_profiles(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd, & + eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, & + p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, & + psipol0,chipol0, jphi,jcd,dpdv,currins,pins,pabs,icd, & + jphip,dpdvp, & + rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, & + outp,rtrp,hcdp,ierr, rhout) + use const_and_precisions, only : zero, one, degree + use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl + use dispersion, only : expinit + use gray_params, only : eqparam_type, prfparam_type, outparam_type, & + rtrparam_type, hcdparam_type, antctrl_type, set_codepar, print_params, & + iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl + use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff + use beamdata, only : pweight, rayi2jk + use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, & + zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q + use errcodes, only : check_err, print_errn, print_errhcd + use magsurf_data, only : flux_average, dealloc_surfvec + use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst + use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & + rhop_tab, rhot_tab + use limiter, only : set_lim, unset_lim + use utils, only : vmaxmin + implicit none +! arguments + type(eqparam_type), intent(in) :: eqp + type(prfparam_type), intent(in) :: prfp + type(outparam_type), intent(in) :: outp + type(rtrparam_type), intent(in) :: rtrp + type(hcdparam_type), intent(in) :: hcdp + + real(wp_), dimension(:), intent(in) :: psrad, terad, derad, zfc + real(wp_), dimension(:), intent(in) :: rv, zv, psinr, fpol, qpsi + real(wp_), dimension(:), intent(in) :: rbnd, zbnd, rlim, zlim + real(wp_), dimension(:,:), intent(in) :: psin + real(wp_), intent(in) :: psia, rvac, rax, zax + integer, intent(in) :: iox0 + real(wp_), intent(in) :: p0, fghz, psipol0, chipol0 + real(wp_), intent(in) :: alpha0,beta0, w1,w2, ri1,ri2, phiw,phir + real(wp_), dimension(3), intent(in) :: xv0 + + real(wp_), intent(in) :: pabs,icd + real(wp_), dimension(:), intent(in) :: jphi,jcd,dpdv,currins,pins + real(wp_), intent(out) :: jphip,dpdvp, & + rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx + + real(wp_), dimension(:), intent(in), optional :: rhout + + integer, intent(out) :: ierr +! local variables + real(wp_), parameter :: taucr = 12._wp_ + + real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre + real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0 + real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx + real(wp_) :: drhotp,drhotj,dpdvmx,jphimx + + real(wp_), dimension(3) :: xv,anv0,anv,bv + real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null() + real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null() + integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1 + logical :: ins_pl, somein, allout + + real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null() + real(wp_), dimension(:), pointer :: tau0=>null(),alphaabs0=>null(),dids0=>null(), & + ccci0=>null() + real(wp_), dimension(:), pointer :: p0jk=>null() + complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null() + integer, dimension(:), pointer :: iiv=>null() + +! parameters log in file headers + character(len=headw), dimension(headl) :: strheader + type(antctrl_type) :: antp + real(wp_) :: rwall + +! ======== set environment BEGIN ======== + call set_codepar(eqp,prfp,outp,rtrp,hcdp) + + call set_lim(rlim,zlim) + + if(iequil<2) then + call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3)) + else + call set_eqspl(rv,zv,psin, psia, psinr,fpol, qpsi, eqp%ssplps,eqp%ssplf, & + rvac, rax,zax, rbnd,zbnd, eqp%ixp) + ! qpsi used for rho_pol/rho_tor mapping (initializes fq,frhotor,frhopol) + end if +! compute flux surface averaged quantities + call flux_average ! requires frhotor for dadrhot,dvdrhot + + if(iprof==0) then + call set_prfan(terad,derad,zfc) + else + call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd) + end if + + call xgygcoeff(fghz,ak0,bres,xgcn) + call launchangles2n(alpha0,beta0,xv0,anv0) + call init_btr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) + + if(iwarm > 1) call expinit + +! ======= set environment END ====== + +! ======= pre-proc prints BEGIN ====== + antp%alpha=alpha0 + antp%beta=beta0 + antp%power=p0 + antp%psi=psipol0 + antp%chi=chipol0 + antp%iox=iox0 +!!!!! missing values + antp%ibeam=0 + antp%filenm='' + rwall=0._wp_ + call print_params(rtrp,hcdp,antp,eqp,rwall,prfp,outp,strheader) + call print_headers(strheader, 0) +! 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 + call print_bres(bres) + call print_prof + call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2), & + sin(beta0*degree)) +! ======= pre-proc prints END ====== + +! ======= post-proc BEGIN ====== + +! compute power and current density profiles for all rays + call pec_init(ipec,rhout) + nnd=size(rhop_tab) +! print power and current density profiles + 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,p0,zero,zero) ! cpl1=cpl2=0 for a single pass + +! ======= post-proc END ====== + +! ======= free memory BEGIN ====== + call unset_eqspl + call unset_q + call unset_rhospl + call unset_prfspl + 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 +! ======= free memory END ====== + end subroutine sum_profiles + + subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) use const_and_precisions, only : wp_, zero implicit none diff --git a/src/main-sum.f90 b/src/main-sum.f90 deleted file mode 100644 index 73c7ff6..0000000 --- a/src/main-sum.f90 +++ /dev/null @@ -1,195 +0,0 @@ -program main_sum - use const_and_precisions, only : wp_,one - use sumcore, only : sum_profiles - 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 - use equilibrium, only : read_equil_an,read_eqdsk,change_cocos,eq_scal, & - set_rhospl,setqphi_num,frhopolv - use coreprofiles, only : read_profiles_an,read_profiles,tene_scal - use reflections, only : range2rect - implicit none - type(antctrl_type) :: antp - type(eqparam_type) :: eqp - type(prfparam_type) :: prfp - type(outparam_type) :: outp - type(rtrparam_type) :: rtrp - type(hcdparam_type) :: hcdp - - real(wp_), dimension(:), allocatable :: psrad, terad, derad, zfc - real(wp_), dimension(:), allocatable :: rv, zv, psinr, fpol, qpsi - real(wp_), dimension(:), allocatable :: rbnd, zbnd, rlim, zlim - real(wp_), dimension(:,:), allocatable :: psin - real(wp_) :: psia, rvac, rax, zax - integer :: iox0 - real(wp_) :: p0mw, fghz, psipol0, chipol0 - real(wp_) :: alpha0, beta0, x0, y0, z0, w1, w2, ri1, ri2, phiw, phir - - real(wp_) :: pec,icd - - integer :: ierr - real(wp_), dimension(:), allocatable :: xrad, rhot - real(wp_), dimension(:), allocatable :: jphi,jcd,dpdv,currins,pins,rtin,rpin - real(wp_) :: rwallm, rmxm, r0m, z0m, dzmx - - integer :: i,j,k,n,ngam,irt - character(len=255) :: fn,dumstr - real(wp_), dimension(5) :: f48v - real(wp_) :: gam,alp,bet, jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx - - -! ======= read parameters BEGIN ======= - call read_params('gray_params.data',rtrp,hcdp,antp,eqp,rwallm,prfp,outp) -! ======= read parameters END ======= - -! ======= read input data BEGIN ======= -!------------ equilibrium ------------ - if(eqp%iequil<2) then - call read_equil_an(eqp%filenm, rtrp%ipass, rv, zv, fpol, qpsi, rlim, zlim) -! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0) - psia = sign(one,qpsi(2)*fpol(1)) - else - call read_eqdsk(eqp%filenm, rv,zv,psin, psia, psinr,fpol,qpsi, rvac, & - rax,zax, rbnd,zbnd, rlim,zlim, eqp%ipsinorm,eqp%idesc,eqp%ifreefmt) - call change_cocos(psia, fpol, qpsi, eqp%icocos, 3) - end if -! re-scale B/I and/or force signs. If sgn=0 on input, set to fpol/-psia signs on output - call eq_scal(psia, fpol, eqp%sgni, eqp%sgnb, eqp%factb) -! ??? 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) - else - call read_profiles(prfp%filenm, xrad, terad, derad, zfc) - allocate(psrad(size(xrad))) - if(prfp%irho==0) then ! xrad==rhot - allocate(rhot(size(psinr))) - call setqphi_num(psinr,qpsi,psia,rhot) - call set_rhospl(sqrt(psinr),rhot) - deallocate(rhot) - psrad=frhopolv(xrad)**2 - else if(prfp%irho == 1) then ! xrad==rhop - psrad=xrad**2 - else - psrad=xrad - end if - deallocate(xrad) - end if -! re-scale input data - call tene_scal(terad,derad,prfp%factte,prfp%factne,eqp%factb,prfp%iscal, & - prfp%iprof) -!------------- antenna -------------- -! interpolate beam table if antctrl%ibeam>0 - select case (antp%ibeam) - case (2) -! to be completed: now 1st beamd always selected, iox read from table - call read_beam2(antp%filenm,1,antp%alpha,antp%beta,fghz,antp%iox,x0,y0,z0, & - w1,w2,ri1,ri2,phiw,phir) - case (1) - call read_beam1(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0, & - w1,w2,ri1,ri2,phiw,phir) - case default - call read_beam0(antp%filenm,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir) - end select - alpha0=antp%alpha - beta0=antp%beta - p0mw=antp%power - psipol0=antp%psi - chipol0=antp%chi - iox0=antp%iox -!--------------- wall --------------- -! set simple limiter if not read from EQDSK -! need to clean up... - r0m=sqrt(x0**2+y0**2)*0.01_wp_ - dzmx=abs(rtrp%ipass)*rtrp%dst*rtrp%nstep*0.01_wp_ - z0m=z0*0.01_wp_ - if (.not.allocated(rlim).or.rtrp%ipass<0) then - if (allocated(rlim)) deallocate(rlim) - if (allocated(zlim)) deallocate(zlim) - allocate(rlim(5)) - allocate(zlim(5)) - if (rtrp%ipass<0) rtrp%ipass = -rtrp%ipass - if(eqp%iequil<2) then - rmxm=(rv(1)+rv(2))*0.01_wp_ - else - rmxm=rv(size(rv)) - end if - call range2rect(rwallm,max(r0m,rmxm),z0m-dzmx,z0m+dzmx,rlim,zlim) - end if -! ======= read input data END ======= - -! ========================= MAIN SUBROUTINE CALL ========================= - allocate(jphi(outp%nrho),jcd(outp%nrho),dpdv(outp%nrho), & - currins(outp%nrho),pins(outp%nrho),rtin(outp%nrho),rpin(outp%nrho)) - - open(100,file='filelist.txt',action='read',status='old') - read(100,*) n,ngam - do i=1,n - read(100,*) fn - open(100+i,file=fn,action='read',status='old') - do j=1,22 - read(100+i,*) dumstr - end do - end do - close(100) - - open(100+n+1,file='f48sum.txt',action='write',status='unknown') - open(100+n+2,file='f7sum.txt',action='write',status='unknown') - - do k=1,ngam - jphi=0.0_wp_ - jcd=0.0_wp_ - dpdv=0.0_wp_ - currins=0.0_wp_ - pins=0.0_wp_ - do j=1,outp%nrho - do i=1,n - read(100+i,*) gam,alp,bet,rpin(j),rtin(j),f48v(1:5),irt - jphi(j)=jphi(j)+f48v(1) - jcd(j)=jcd(j)+f48v(2) - dpdv(j)=dpdv(j)+f48v(3) - currins(j)=currins(j)+f48v(4) - pins(j)=pins(j)+f48v(5) - end do - write(100+n+1,'(10(1x,e16.8e3),i5)') gam,alp,bet,rpin(j),rtin(j), & - jphi(j),jcd(j),dpdv(j),currins(j),pins(j),irt - end do - pec=pins(outp%nrho) - icd=currins(outp%nrho) - write(100+n+1,*) - call sum_profiles(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, jphi,jcd,dpdv,currins,pins,pec,icd, & - jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, & - outp,rtrp,hcdp,ierr) - write(100+n+2,'(15(1x,e12.5),i5,4(1x,e12.5))') gam,alp,bet,icd,pec, & - jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx - end do - do i=1,n+2 - close(100+i) - end do -! ======================================================================== - -! ======= control prints BEGIN ======= - if(ierr/=0) print*,' IERR = ', ierr - print*,' ' - print*,'Pabs (MW) = ', pec - print*,'Icd (kA) = ', icd*1.0e3_wp_ -! ======= control prints END ======= - -! ======= free memory BEGIN ======= - if(allocated(psrad)) deallocate(psrad) - if(allocated(terad)) deallocate(terad, derad, zfc) - if(allocated(rv)) deallocate(rv, zv, fpol, qpsi) - if(allocated(psin)) deallocate(psin, psinr) - if(allocated(rbnd)) deallocate(rbnd,zbnd) - if(allocated(rlim)) deallocate(rlim,zlim) - if(allocated(dpdv)) deallocate(jphi,jcd,dpdv,currins,pins,rtin,rpin) -! ======= free memory END ====== -end program main_sum \ No newline at end of file diff --git a/src/main.f90 b/src/main.f90 index 8dc5c7e..1807d3e 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,16 +1,18 @@ program main_std use const_and_precisions, only : wp_,one - 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 - use equilibrium, only : read_equil_an,read_eqdsk,change_cocos,eq_scal, & - set_rhospl,setqphi_num,frhopolv - use coreprofiles, only : read_profiles_an,read_profiles,tene_scal - use reflections, only : range2rect + use graycore, only : gray_main, sum_profiles + 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 + use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, & + eq_scal, set_rhospl, setqphi_num, frhopolv + use coreprofiles, only : read_profiles_an, read_profiles, tene_scal + use reflections, only : range2rect + implicit none - type(antctrl_type) :: antp - type(eqparam_type) :: eqp + type(antctrl_type) :: antp + type(eqparam_type) :: eqp type(prfparam_type) :: prfp type(outparam_type) :: outp type(rtrparam_type) :: rtrp @@ -21,9 +23,8 @@ program main_std real(wp_), dimension(:), allocatable :: rbnd, zbnd, rlim, zlim real(wp_), dimension(:,:), allocatable :: psin real(wp_) :: psia, rvac, rax, zax - integer :: iox0 - real(wp_) :: p0mw, fghz, psipol0, chipol0 - real(wp_) :: alpha0, beta0, x0, y0, z0, w1, w2, ri1, ri2, phiw, phir + real(wp_) :: fghz + real(wp_) :: x0, y0, z0, w1, w2, ri1, ri2, phiw, phir real(wp_) :: pec,icd @@ -31,109 +32,184 @@ program main_std real(wp_), dimension(:), allocatable :: xrad, rhot, dpdv, jcd real(wp_) :: rwallm, rmxm, r0m, z0m, dzmx -! ======= read parameters BEGIN ======= - call read_params('gray_params.data',rtrp,hcdp,antp,eqp,rwallm,prfp,outp) -! ======= read parameters END ======= + logical :: sum_mode = .FALSE. -! ======= read input data BEGIN ======= -!------------ equilibrium ------------ + ! ------- sum mode variables ------- + real(wp_), dimension(:), allocatable :: jphi, currins, pins, rtin, rpin + integer :: i,j,k,n,ngam,irt + character(len=255) :: fn,dumstr + real(wp_), dimension(5) :: f48v + real(wp_) :: gam,alp,bet, jphip,dpdvp, & + rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx + ! ---------------------------------- + + call read_params('gray_params.data', rtrp, hcdp, antp, eqp, rwallm, prfp, outp) + + ! ======= read input data BEGIN ======= + + !------------ equilibrium ------------ if(eqp%iequil<2) then call read_equil_an(eqp%filenm, rtrp%ipass, rv, zv, fpol, qpsi, rlim, zlim) -! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0) + ! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0) psia = sign(one,qpsi(2)*fpol(1)) else - call read_eqdsk(eqp%filenm, rv,zv,psin, psia, psinr,fpol,qpsi, rvac, & - rax,zax, rbnd,zbnd, rlim,zlim, eqp%ipsinorm,eqp%idesc,eqp%ifreefmt) + call read_eqdsk(eqp%filenm, rv,zv,psin, psia, psinr,fpol,qpsi, rvac, & + rax,zax, rbnd,zbnd, rlim,zlim, & + eqp%ipsinorm,eqp%idesc,eqp%ifreefmt) call change_cocos(psia, fpol, qpsi, eqp%icocos, 3) end if -! re-scale B/I and/or force signs. If sgn=0 on input, set to fpol/-psia signs on output + ! 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) -! ??? 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 + ! ??? 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) else call read_profiles(prfp%filenm, xrad, terad, derad, zfc) allocate(psrad(size(xrad))) - if(prfp%irho==0) then ! xrad==rhot + if(prfp%irho == 0) then ! xrad==rhot allocate(rhot(size(psinr))) call setqphi_num(psinr,qpsi,psia,rhot) call set_rhospl(sqrt(psinr),rhot) deallocate(rhot) - psrad=frhopolv(xrad)**2 + psrad = frhopolv(xrad)**2 else if(prfp%irho == 1) then ! xrad==rhop - psrad=xrad**2 + psrad = xrad**2 else - psrad=xrad + psrad = xrad end if deallocate(xrad) end if -! re-scale input data - call tene_scal(terad,derad,prfp%factte,prfp%factne,eqp%factb,prfp%iscal, & - prfp%iprof) -!------------- antenna -------------- -! interpolate beam table if antctrl%ibeam>0 + ! re-scale input data + call tene_scal(terad, derad, prfp%factte, prfp%factne, & + eqp%factb, prfp%iscal, prfp%iprof) + ! ---------------------------------- + + !------------- antenna -------------- + ! interpolate beam table if antctrl%ibeam>0 select case (antp%ibeam) - case (2) -! to be completed: now 1st beamd always selected, iox read from table - call read_beam2(antp%filenm,1,antp%alpha,antp%beta,fghz,antp%iox,x0,y0,z0, & - w1,w2,ri1,ri2,phiw,phir) - case (1) - call read_beam1(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0, & - w1,w2,ri1,ri2,phiw,phir) - case default - call read_beam0(antp%filenm,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir) + case (2) + ! to be completed: now 1st beamd always selected, iox read from table + call read_beam2(antp%filenm, 1, antp%alpha, antp%beta, & + fghz, antp%iox, x0, y0, z0, & + w1, w2, ri1, ri2, phiw, phir) + case (1) + call read_beam1(antp%filenm, antp%alpha, antp%beta, & + fghz, x0, y0, z0, & + w1, w2, ri1, ri2, phiw, phir) + case default + call read_beam0(antp%filenm, & + fghz, x0, y0, z0, & + w1, w2, ri1, ri2, phiw, phir) end select - alpha0=antp%alpha - beta0=antp%beta - p0mw=antp%power - psipol0=antp%psi - chipol0=antp%chi - iox0=antp%iox -!--------------- wall --------------- -! set simple limiter if not read from EQDSK -! need to clean up... - r0m=sqrt(x0**2+y0**2)*0.01_wp_ - dzmx=abs(rtrp%ipass)*rtrp%dst*rtrp%nstep*0.01_wp_ - z0m=z0*0.01_wp_ - if (.not.allocated(rlim).or.rtrp%ipass<0) then - if (allocated(rlim)) deallocate(rlim) - if (allocated(zlim)) deallocate(zlim) - allocate(rlim(5)) - allocate(zlim(5)) - if (rtrp%ipass<0) rtrp%ipass = -rtrp%ipass - if(eqp%iequil<2) then - rmxm=(rv(1)+rv(2))*0.01_wp_ - else - rmxm=rv(size(rv)) + ! ---------------------------------- + + !--------------- wall --------------- + ! set simple limiter if not read from EQDSK + ! need to clean up... + r0m=sqrt(x0**2+y0**2)*0.01_wp_ + dzmx=abs(rtrp%ipass)*rtrp%dst*rtrp%nstep*0.01_wp_ + z0m=z0*0.01_wp_ + if (.not.allocated(rlim).or.rtrp%ipass<0) then + if (allocated(rlim)) deallocate(rlim) + if (allocated(zlim)) deallocate(zlim) + allocate(rlim(5)) + allocate(zlim(5)) + if (rtrp%ipass<0) rtrp%ipass = -rtrp%ipass + if(eqp%iequil<2) then + rmxm=(rv(1)+rv(2))*0.01_wp_ + else + rmxm=rv(size(rv)) + end if + call range2rect(rwallm,max(r0m,rmxm),z0m-dzmx,z0m+dzmx,rlim,zlim) end if - call range2rect(rwallm,max(r0m,rmxm),z0m-dzmx,z0m+dzmx,rlim,zlim) + ! ---------------------------------- + + ! ======= read input data END ======= + + ! ========================= MAIN SUBROUTINE CALL ========================= + allocate(dpdv(outp%nrho), jcd(outp%nrho)) + if (sum_mode) then + allocate(jphi(outp%nrho), currins(outp%nrho), & + pins(outp%nrho), rtin(outp%nrho), rpin(outp%nrho)) + + open(100,file='filelist.txt',action='read',status='old') + read(100,*) n,ngam + do i=1,n + read(100,*) fn + open(100+i,file=fn,action='read',status='old') + do j=1,22 + read(100+i,*) dumstr + end do + end do + close(100) + + open(100+n+1,file='f48sum.txt',action='write',status='unknown') + open(100+n+2,file='f7sum.txt',action='write',status='unknown') + + do k=1,ngam + jphi=0.0_wp_ + jcd=0.0_wp_ + dpdv=0.0_wp_ + currins=0.0_wp_ + pins=0.0_wp_ + do j=1,outp%nrho + do i=1,n + read(100+i,*) gam,alp,bet,rpin(j),rtin(j),f48v(1:5),irt + jphi(j)=jphi(j)+f48v(1) + jcd(j)=jcd(j)+f48v(2) + dpdv(j)=dpdv(j)+f48v(3) + currins(j)=currins(j)+f48v(4) + pins(j)=pins(j)+f48v(5) + end do + write(100+n+1,'(10(1x,e16.8e3),i5)') gam,alp,bet,rpin(j),rtin(j), & + jphi(j),jcd(j),dpdv(j),currins(j),pins(j),irt + end do + pec=pins(outp%nrho) + icd=currins(outp%nrho) + write(100+n+1,*) + call sum_profiles(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & + psrad,terad,derad,zfc,prfp, rlim,zlim, & + antp%beta,fghz,antp%alpha,antp%beta, & + (/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,antp%iox, & + antp%psi, antp%chi, jphi,jcd,dpdv,currins,pins,pec,icd, & + jphip,dpdvp, rhotj,rhotjava,rhotp,rhotpav, & + drhotjava,drhotpav, ratjamx,ratjbmx, outp,rtrp,hcdp,ierr) + write(100+n+2,'(15(1x,e12.5),i5,4(1x,e12.5))') gam,alp,bet,icd,pec, & + jphip,dpdvp, & + rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx + end do + do i=1,n+2 + close(100+i) + end do + else + call gray_main(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, & + psrad,terad,derad,zfc,prfp, rlim,zlim, & + antp%power,fghz,antp%alpha,antp%beta, & + (/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,antp%iox, & + antp%psi,antp%chi, dpdv,jcd,pec,icd, outp,rtrp,hcdp,ierr) end if -! ======= read input data END ======= + ! ======================================================================== -! ========================= MAIN SUBROUTINE CALL ========================= - allocate(dpdv(outp%nrho),jcd(outp%nrho)) - 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) -! ======================================================================== - -! ======= control prints BEGIN ======= + ! ======= control prints BEGIN ======= if(ierr/=0) print*,' IERR = ', ierr print*,' ' - print*,'Pabs (MW), Icd (kA) = ', pec,icd*1.0e3_wp_ -! ======= control prints END ======= + print*,'Pabs (MW) = ', pec + print*,'Icd (kA) = ', icd*1.0e3_wp_ + ! ======= control prints END ======= -! ======= free memory BEGIN ======= + ! ======= free memory BEGIN ======= if(allocated(psrad)) deallocate(psrad) if(allocated(terad)) deallocate(terad, derad, zfc) if(allocated(rv)) deallocate(rv, zv, fpol, qpsi) if(allocated(psin)) deallocate(psin, psinr) - if(allocated(rbnd)) deallocate(rbnd,zbnd) - if(allocated(rlim)) deallocate(rlim,zlim) + if(allocated(rbnd)) deallocate(rbnd, zbnd) + if(allocated(rlim)) deallocate(rlim, zlim) if(allocated(dpdv)) deallocate(dpdv, jcd) -! ======= free memory END ====== -end program main_std \ No newline at end of file + if(sum_mode) deallocate(jphi, currins, pins, rtin, rpin) + ! ======= free memory END ====== +end program main_std diff --git a/src/sumcore.f90 b/src/sumcore.f90 deleted file mode 100644 index f58ac57..0000000 --- a/src/sumcore.f90 +++ /dev/null @@ -1,1907 +0,0 @@ -module sumcore - use const_and_precisions, only : wp_ - implicit none - -contains - - subroutine sum_profiles(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd, & - eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, & - p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, & - psipol0,chipol0, jphi,jcd,dpdv,currins,pins,pabs,icd, & - jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, & - outp,rtrp,hcdp,ierr, rhout) - use const_and_precisions, only : zero, one, degree - use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl - use dispersion, only : expinit - use gray_params, only : eqparam_type, prfparam_type, outparam_type, & - rtrparam_type, hcdparam_type, antctrl_type, set_codepar, print_params, & - iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl - use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff - use beamdata, only : pweight, rayi2jk - use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, & - zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q - use errcodes, only : check_err, print_errn, print_errhcd - use magsurf_data, only : flux_average, dealloc_surfvec - use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst - use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & - rhop_tab, rhot_tab - use limiter, only : set_lim, unset_lim - use utils, only : vmaxmin - implicit none -! arguments - type(eqparam_type), intent(in) :: eqp - type(prfparam_type), intent(in) :: prfp - type(outparam_type), intent(in) :: outp - type(rtrparam_type), intent(in) :: rtrp - type(hcdparam_type), intent(in) :: hcdp - - real(wp_), dimension(:), intent(in) :: psrad, terad, derad, zfc - real(wp_), dimension(:), intent(in) :: rv, zv, psinr, fpol, qpsi - real(wp_), dimension(:), intent(in) :: rbnd, zbnd, rlim, zlim - real(wp_), dimension(:,:), intent(in) :: psin - real(wp_), intent(in) :: psia, rvac, rax, zax - integer, intent(in) :: iox0 - real(wp_), intent(in) :: p0, fghz, psipol0, chipol0 - real(wp_), intent(in) :: alpha0,beta0, w1,w2, ri1,ri2, phiw,phir - real(wp_), dimension(3), intent(in) :: xv0 - - real(wp_), intent(in) :: pabs,icd - real(wp_), dimension(:), intent(in) :: jphi,jcd,dpdv,currins,pins - real(wp_), intent(out) :: jphip,dpdvp, & - rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx - - real(wp_), dimension(:), intent(in), optional :: rhout - - integer, intent(out) :: ierr -! local variables - real(wp_), parameter :: taucr = 12._wp_ - - real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre - real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0 - real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx - real(wp_) :: drhotp,drhotj,dpdvmx,jphimx - - real(wp_), dimension(3) :: xv,anv0,anv,bv - real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null() - real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null() - integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1 - logical :: ins_pl, somein, allout - - real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null() - real(wp_), dimension(:), pointer :: tau0=>null(),alphaabs0=>null(),dids0=>null(), & - ccci0=>null() - real(wp_), dimension(:), pointer :: p0jk=>null() - complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null() - integer, dimension(:), pointer :: iiv=>null() - -! parameters log in file headers - character(len=headw), dimension(headl) :: strheader - type(antctrl_type) :: antp - real(wp_) :: rwall - -! ======= set environment BEGIN ====== - call set_codepar(eqp,prfp,outp,rtrp,hcdp) - - call set_lim(rlim,zlim) - - if(iequil<2) then - call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3)) - else - call set_eqspl(rv,zv,psin, psia, psinr,fpol, qpsi, eqp%ssplps,eqp%ssplf, & - rvac, rax,zax, rbnd,zbnd, eqp%ixp) - ! qpsi used for rho_pol/rho_tor mapping (initializes fq,frhotor,frhopol) - end if -! compute flux surface averaged quantities - call flux_average ! requires frhotor for dadrhot,dvdrhot - - if(iprof==0) then - call set_prfan(terad,derad,zfc) - else - call set_prfspl(psrad, terad, derad, zfc, prfp%sspld, prfp%psnbnd) - end if - - call xgygcoeff(fghz,ak0,bres,xgcn) - call launchangles2n(alpha0,beta0,xv0,anv0) - call init_btr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - - if(iwarm > 1) call expinit - -! ======= set environment END ====== - -! ======= pre-proc prints BEGIN ====== - antp%alpha=alpha0 - antp%beta=beta0 - antp%power=p0 - antp%psi=psipol0 - antp%chi=chipol0 - antp%iox=iox0 -!!!!! missing values - antp%ibeam=0 - antp%filenm='' - rwall=0._wp_ - call print_params(rtrp,hcdp,antp,eqp,rwall,prfp,outp,strheader) - call print_headers(strheader) -! 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 - call print_bres(bres) - call print_prof - call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2), & - sin(beta0*degree)) -! ======= pre-proc prints END ====== - -! ======= post-proc BEGIN ====== - -! compute power and current density profiles for all rays - call pec_init(ipec,rhout) - nnd=size(rhop_tab) -! print power and current density profiles - 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 unset_eqspl - call unset_q - call unset_rhospl - call unset_prfspl - 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 -! ======= free memory END ====== - end subroutine sum_profiles - - - - subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) - use const_and_precisions, only : wp_, zero - implicit none -! arguments - real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci - real(wp_), dimension(:), intent(out) :: tau0,alphaabs0,dids0,ccci0 - integer, dimension(:), intent(out) :: iiv -!! common/external functions/variables -! integer :: jclosest -! real(wp_), dimension(3) :: anwcl,xwcl -! -! common/refln/anwcl,xwcl,jclosest -! -! jclosest=nrayr+1 -! anwcl(1:3)=0.0_wp_ -! xwcl(1:3)=0.0_wp_ - - psjki = zero - ppabs = zero - ccci = zero - tau0 = zero - alphaabs0 = zero - dids0 = zero - ccci0 = zero - iiv = 1 - - 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 -! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! - use const_and_precisions, only : wp_,izero,zero,one,pi,half,two,degree,ui=>im - use math, only : catand - use gray_params, only : idst - use beamdata, only : nray,nrayr,nrayth,rwmax - implicit none -! arguments - real(wp_), dimension(3), intent(in) :: xv0c,anv0c - real(wp_), intent(in) :: ak0 - real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir - real(wp_), dimension(6,nray), intent(out) :: ywrk0,ypwrk0 - real(wp_), dimension(3,nray), intent(out) :: gri - real(wp_), dimension(3,3,nray), intent(out) :: ggri - real(wp_), dimension(3,nrayth,nrayr), intent(out) :: xc0,du10 - -! local variables - integer :: j,k,jk - real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw,alfak - real(wp_) :: wwcsi,wweta,sk,sw,dk,dw,rci1,ww1,rci2,ww2,wwxx,wwyy,wwxy - real(wp_) :: rcixx,rciyy,rcixy,dwwxx,dwwyy,dwwxy,d2wwxx,d2wwyy,d2wwxy - real(wp_) :: drcixx,drciyy,drcixy,dr,da,ddfu,dcsiw,detaw,dx0t,dy0t - real(wp_) :: x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0,gxt,gyt,gzt,gr2 - real(wp_) :: gxxt,gyyt,gzzt,gxyt,gxzt,gyzt,dgr2xt,dgr2yt,dgr2zt - real(wp_) :: dgr2x,dgr2y,dgr2z,pppx,pppy,denpp,ppx,ppy - real(wp_) :: anzt,anxt,anyt,anx,any,anz,an20,an0 - real(wp_) :: du1tx,du1ty,du1tz,denom,ddr,ddi - real(wp_), dimension(nrayr) :: uj - real(wp_), dimension(nrayth) :: sna,csa - complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1,dqi2 - complex(wp_) :: dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy - - csth=anv0c(3) - snth=sqrt(one-csth**2) - if(snth > zero) then - csps=anv0c(2)/snth - snps=anv0c(1)/snth - else - csps=one - snps=zero - end if - -! Gaussian beam: exp[-ik0 zt] exp[-i k0/2 S(xt,yt,zt)] -! xt,yt,zt, cartesian coordinate system with zt along the beamline and xt in the z = 0 plane -! S(xt,yt,zt) = S_real +i S_imag = Qxx(zt) xt^2 + Qyy(zt) yt^2 + 2 Qxy(zt) xt yt -! (csiw, etaw) and (csiR, etaR) intensity and phase ellipse, rotated by angle phiw and phiR -! S(xt,yt,zt) = csiR^2 / Rccsi +etaR^2 /Rceta - i (csiw^2 Wcsi +etaw^2 Weta) -! Rccsi,eta curvature radius at the launching point -! Wcsi,eta =2/(k0 wcsi,eta^2) with wcsi,eta^2 beam size at the launching point - - phiwrad = phiw*degree - phirrad = phir*degree - csphiw = cos(phiwrad) - snphiw = sin(phiwrad) -! csphir = cos(phirrad) -! snphir = sin(phirrad) - - wwcsi = two/(ak0*wcsi**2) - wweta = two/(ak0*weta**2) - - if(phir/=phiw) then - sk = rcicsi + rcieta - sw = wwcsi + wweta - dk = rcicsi - rcieta - dw = wwcsi - wweta - ts = -(dk*sin(2*phirrad) - ui*dw*sin(2*phiwrad)) - tc = (dk*cos(2*phirrad) - ui*dw*cos(2*phiwrad)) - phic = half*catand(ts/tc) - ddd = dk*cos(2*(phirrad+phic)) - ui*dw*cos(2*(phiwrad+phic)) - sss = sk - ui*sw - qi1 = half*(sss + ddd) - qi2 = half*(sss - ddd) - rci1 = dble(qi1) - rci2 = dble(qi2) - ww1 = -dimag(qi1) - ww2 = -dimag(qi2) - else - rci1 = rcicsi - rci2 = rcieta - ww1 = wwcsi - ww2 = wweta - phic = -phiwrad - qi1 = rci1 - ui*ww1 - qi2 = rci2 - ui*ww2 - end if - -! w01=sqrt(2.0_wp_/(ak0*ww1)) -! d01=-rci1/(rci1**2+ww1**2) -! w02=sqrt(2.0_wp_/(ak0*ww2)) -! d02=-rci2/(rci2**2+ww2**2) - - qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2 - qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2 - qqxy = -(qi1 - qi2)*sin(phic)*cos(phic) - wwxx = -dimag(qqxx) - wwyy = -dimag(qqyy) - wwxy = -dimag(qqxy) - rcixx = dble(qqxx) - rciyy = dble(qqyy) - rcixy = dble(qqxy) - - dqi1 = -qi1**2 - dqi2 = -qi2**2 - d2qi1 = 2*qi1**3 - d2qi2 = 2*qi2**3 - dqqxx = dqi1*cos(phic)**2 + dqi2*sin(phic)**2 - dqqyy = dqi1*sin(phic)**2 + dqi2*cos(phic)**2 - dqqxy = -(dqi1 - dqi2)*sin(phic)*cos(phic) - d2qqxx = d2qi1*cos(phic)**2 + d2qi2*sin(phic)**2 - d2qqyy = d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2 - d2qqxy = -(d2qi1 - d2qi2)*sin(phic)*cos(phic) - - dwwxx = -dimag(dqqxx) - dwwyy = -dimag(dqqyy) - dwwxy = -dimag(dqqxy) - d2wwxx = -dimag(d2qqxx) - d2wwyy = -dimag(d2qqyy) - d2wwxy = -dimag(d2qqxy) - drcixx = dble(dqqxx) - drciyy = dble(dqqyy) - drcixy = dble(dqqxy) - - if(nrayr > 1) then - dr = rwmax/dble(nrayr-1) - else - dr = one - end if - ddfu = two*dr**2/ak0 ! twodr2 = 2*dr**2 = 2*rwmax/dble(nrayr-1) - do j = 1, nrayr - uj(j) = dble(j-1) - end do - - da=2*pi/dble(nrayth) - do k=1,nrayth - alfak = (k-1)*da - sna(k) = sin(alfak) - csa(k) = cos(alfak) - end do - -! central ray - jk=1 - gri(:,1) = zero - ggri(:,:,1) = zero - - ywrk0(1:3,1) = xv0c - ywrk0(4:6,1) = anv0c - ypwrk0(1:3,1) = anv0c - ypwrk0(4:6,1) = zero - - do k=1,nrayth - dcsiw = dr*csa(k)*wcsi - detaw = dr*sna(k)*weta - dx0t = dcsiw*csphiw - detaw*snphiw - dy0t = dcsiw*snphiw + detaw*csphiw - du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu - du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu - - xc0(:,k,1) = xv0c - du10(1,k,1) = du1tx*csps + snps*du1ty*csth - du10(2,k,1) = -du1tx*snps + csps*du1ty*csth - du10(3,k,1) = -du1ty*snth - end do - ddr = zero - ddi = zero - -! loop on rays jk>1 - j=2 - k=0 - do jk=2,nray - k=k+1 - if(k > nrayth) then - j=j+1 - k=1 - end if - -! csiw=u*dcsiw -! etaw=u*detaw -! csir=x0t*csphir+y0t*snphir -! etar=-x0t*snphir+y0t*csphir - dcsiw = dr*csa(k)*wcsi - detaw = dr*sna(k)*weta - dx0t = dcsiw*csphiw - detaw*snphiw - dy0t = dcsiw*snphiw + detaw*csphiw - x0t = uj(j)*dx0t - y0t = uj(j)*dy0t - 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) - dz0 = z0t*csth - y0t*snth - x0 = xv0c(1) + dx0 - y0 = xv0c(2) + dy0 - z0 = xv0c(3) + dz0 - - gxt = x0t*wwxx + y0t*wwxy - gyt = x0t*wwxy + y0t*wwyy - gzt = half*(x0t**2*dwwxx + y0t**2*dwwyy ) + x0t*y0t*dwwxy - gr2 = gxt*gxt + gyt*gyt + gzt*gzt - gxxt = wwxx - gyyt = wwyy - gzzt = half*(x0t**2*d2wwxx + y0t**2*d2wwyy) + x0t*y0t*d2wwxy - gxyt = wwxy - gxzt = x0t*dwwxx + y0t*dwwxy - gyzt = x0t*dwwxy + y0t*dwwyy - dgr2xt = 2*(gxt*gxxt + gyt*gxyt + gzt*gxzt) - dgr2yt = 2*(gxt*gxyt + gyt*gyyt + gzt*gyzt) - dgr2zt = 2*(gxt*gxzt + gyt*gyzt + gzt*gzzt) - dgr2x = dgr2xt*csps + snps*(dgr2yt*csth + dgr2zt*snth) - dgr2y = -dgr2xt*snps + csps*(dgr2yt*csth + dgr2zt*snth) - dgr2z = dgr2zt*csth - dgr2yt*snth - - gri(1,jk) = gxt*csps + snps*(gyt*csth + gzt*snth) - gri(2,jk) = -gxt*snps + csps*(gyt*csth + gzt*snth) - gri(3,jk) = gzt*csth - gyt*snth - ggri(1,1,jk) = gxxt*csps**2 & - + snps**2 *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) & - +2*snps*csps*(gxyt*csth + gxzt*snth) - ggri(2,1,jk) = csps*snps & - *(-gxxt+csth**2*gyyt + snth**2*gzzt + 2*csth*snth*gyzt) & - +(csps**2 - snps**2)*(snth*gxzt + csth*gxyt) - ggri(3,1,jk) = csth*snth*snps*(gzzt - gyyt) + (csth**2 - snth**2) & - *snps*gyzt + csps*(csth*gxzt - snth*gxyt) - ggri(1,2,jk) = ggri(2,1,jk) - ggri(2,2,jk) = gxxt*snps**2 & - + csps**2 *(gyyt*csth**2 + gzzt*snth**2 + 2*snth*csth*gyzt) & - -2*snps*csps*(gxyt*csth + gxzt*snth) - ggri(3,2,jk) = csth*snth*csps*(gzzt - gyyt) + (csth**2-snth**2) & - *csps*gyzt + snps*(snth*gxyt - csth*gxzt) - ggri(1,3,jk) = ggri(3,1,jk) - ggri(2,3,jk) = ggri(3,2,jk) - ggri(3,3,jk) = gzzt*csth**2 + gyyt*snth**2 - 2*csth*snth*gyzt - - du1tx = (dx0t*wwxx + dy0t*wwxy)/ddfu - du1ty = (dx0t*wwxy + dy0t*wwyy)/ddfu - du1tz = half*uj(j)*(dx0t**2*dwwxx + dy0t**2*dwwyy + 2*dx0t*dy0t*dwwxy)/ddfu - - du10(1,k,j) = du1tx*csps + snps*(du1ty*csth + du1tz*snth) - du10(2,k,j) = -du1tx*snps + csps*(du1ty*csth + du1tz*snth) - du10(3,k,j) = du1tz*csth - du1ty*snth - - pppx = x0t*rcixx + y0t*rcixy - pppy = x0t*rcixy + y0t*rciyy - denpp = pppx*gxt + pppy*gyt - if (denpp/=zero) then - ppx = -pppx*gzt/denpp - ppy = -pppy*gzt/denpp - else - ppx = zero - ppy = zero - end if - - anzt = sqrt((one + gr2)/(one + ppx**2 + ppy**2)) - anxt = ppx*anzt - anyt = ppy*anzt - - anx = anxt*csps + snps*(anyt*csth + anzt*snth) - any =-anxt*snps + csps*(anyt*csth + anzt*snth) - anz = anzt*csth - anyt*snth - - an20 = one + gr2 - an0 = sqrt(an20) - - xc0(1,k,j) = x0 - xc0(2,k,j) = y0 - xc0(3,k,j) = z0 - - ywrk0(1,jk) = x0 - ywrk0(2,jk) = y0 - ywrk0(3,jk) = z0 - ywrk0(4,jk) = anx - ywrk0(5,jk) = any - ywrk0(6,jk) = anz - - select case(idst) - case(1) -! integration variable: c*t - denom = one - case(2) -! integration variable: Sr - denom = an20 - case default ! idst=0 -! integration variable: s - denom = an0 - end select - ypwrk0(1,jk) = anx/denom - ypwrk0(2,jk) = any/denom - ypwrk0(3,jk) = anz/denom - ypwrk0(4,jk) = dgr2x/(2*denom) - ypwrk0(5,jk) = dgr2y/(2*denom) - ypwrk0(6,jk) = dgr2z/(2*denom) - - 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,(/zero,zero,zero/), & - ak0,zero,zero,(/zero,zero,zero/),zero,zero,zero,zero,zero,zero, & - 0,0,1,ddr,ddi) ! st=0, index_rt=1, B=0, N=0, psin=-1 - end do - end subroutine ic_gb - - - - subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr) -! Runge-Kutta integrator - use const_and_precisions, only : wp_ -! use gray_params, only : igrad - use beamdata, only : h,hh,h6 - implicit none - real(wp_), intent(in) :: sox,bres,xgcn - real(wp_), dimension(6), intent(inout) :: y - real(wp_), dimension(6), intent(in) :: yp - real(wp_), dimension(3), intent(in) :: dgr - real(wp_), dimension(3,3), intent(in) :: ddgr - - real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4 - real(wp_) :: gr2 - real(wp_), dimension(3) :: dgr2 - -! if(igrad.eq.1) then - gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2 - dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3)) -! end if - fk1 = yp - - yy = y + fk1*hh - call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2) - yy = y + fk2*hh - call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk3) - yy = y + fk3*h - call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk4) - - 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 - use const_and_precisions, only : wp_ - implicit none -! arguments - real(wp_), dimension(6), intent(in) :: y - real(wp_), intent(in) :: sox,bres,xgcn,gr2 - real(wp_), dimension(3), intent(in) :: dgr2,dgr - real(wp_), dimension(3,3), intent(in) :: ddgr - real(wp_), dimension(6), intent(out) :: dery -! local variables - real(wp_) :: psinv,dens,btot,xg,yg,anpl,anpr,ajphi - real(wp_) :: ddr,ddi,dersdst,derdnm - real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg - real(wp_), dimension(3,3) :: derbv - - xv = y(1:3) - call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg, & - ajphi) - - anv = y(4:6) - call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & - dery,anpl,anpr,ddr,ddi,dersdst,derdnm) - end subroutine rhs - - - - subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & - bv,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr) -! Compute right-hand side terms of the ray equations (dery) -! used after full R-K step and grad(S_I) update - use errcodes, only : pnpl - implicit none -! arguments - real(wp_), dimension(3), intent(in) :: xv,anv - real(wp_), dimension(3), intent(in) :: dgr - real(wp_), dimension(3,3), intent(in) :: ddgr - real(wp_), intent(in) :: sox,bres,xgcn - real(wp_), dimension(6), intent(out) :: dery - real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr - real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm - real(wp_), dimension(3), intent(out) :: bv - integer, intent(out) :: ierr -! local variables - real(wp_) :: gr2,ajphi - real(wp_), dimension(3) :: dgr2,derxg,deryg - real(wp_), dimension(3,3) :: derbv - real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ - - gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2 - dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3)) - call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg,ajphi) - call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & - dery,anpl,anpr,ddr,ddi,dersdst,derdnm) - - ierr=0 - if( abs(anpl) > anplth1) then - if(abs(anpl) > anplth2) then - ierr=ibset(ierr,pnpl+1) - else - ierr=ibset(ierr,pnpl) - end if - end if - 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 - implicit none - real(wp_), intent(in) :: ak0 - real(wp_), dimension(6,nray), intent(in) :: ywrk - real(wp_), dimension(3,nrayth,nrayr), intent(inout) :: xc,du1 - real(wp_), dimension(3,nray), intent(out) :: gri - real(wp_), dimension(3,3,nray), intent(out) :: ggri -! local variables - real(wp_), dimension(3,nrayth,nrayr) :: xco,du1o - integer :: jk,j,jm,jp,k,km,kp - real(wp_) :: ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz - real(wp_) :: dfuu,dffiu,gx,gxx,gxy,gxz,gy,gyy,gyz,gz,gzz - real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgu - real(wp_), dimension(3,3) :: dgg,dff - -! update position and du1 vectors - xco = xc - du1o = du1 - - jk = 1 - do j=1,nrayr - do k=1,nrayth - if(j>1) jk=jk+1 - xc(1:3,k,j)=ywrk(1:3,jk) - end do - end do - -! compute grad u1 for central ray - j = 1 - jp = 2 - do k=1,nrayth - if(k == 1) then - km = nrayth - else - km = k-1 - end if - if(k == nrayth) then - kp = 1 - else - kp = k+1 - end if - dxv1 = xc(:,k ,jp) - xc(:,k ,j) - dxv2 = xc(:,kp,jp) - xc(:,km,jp) - dxv3 = xc(:,k ,j) - xco(:,k ,j) - call solg0(dxv1,dxv2,dxv3,dgu) - du1(:,k,j) = dgu - end do - gri(:,1) = zero - -! compute grad u1 and grad(S_I) for all the other rays - dfuu=twodr2/ak0 ! twodr2 = 2*dr**2 = 2*(rwmax/(nrayr-1))**2 - jm=1 - j=2 - k=0 - dffiu = dfuu - do jk=2,nray - k=k+1 - if(k > nrayth) then - jm = j - j = j+1 - k = 1 - dffiu = dfuu*jm - end if - kp = k+1 - km = k-1 - if (k == 1) then - km=nrayth - else if (k == nrayth) then - kp=1 - end if - dxv1 = xc(:,k ,j) - xc(:,k ,jm) - dxv2 = xc(:,kp,j) - xc(:,km,j) - dxv3 = xc(:,k ,j) - xco(:,k ,j) - call solg0(dxv1,dxv2,dxv3,dgu) - du1(:,k,j) = dgu - gri(:,jk) = dgu(:)*dffiu - end do - -! compute derivatives of grad u and grad(S_I) for rays jk>1 - ggri(:,:,1) = zero - jm=1 - j=2 - k=0 - dffiu = dfuu - do jk=2,nray - k=k+1 - if(k > nrayth) then - jm=j - j=j+1 - k=1 - dffiu = dfuu*jm - end if - kp=k+1 - km=k-1 - if (k == 1) then - km=nrayth - else if (k == nrayth) then - kp=1 - end if - dxv1 = xc(:,k ,j) - xc(:,k ,jm) - dxv2 = xc(:,kp,j) - xc(:,km,j) - dxv3 = xc(:,k ,j) - xco(:,k ,j) - dff(:,1) = du1(:,k ,j) - du1(:,k ,jm) - dff(:,2) = du1(:,kp,j) - du1(:,km,j) - dff(:,3) = du1(:,k ,j) - du1o(:,k ,j) - call solg3(dxv1,dxv2,dxv3,dff,dgg) - -! derivatives of u - ux = du1(1,k,j) - uy = du1(2,k,j) - uz = du1(3,k,j) - uxx = dgg(1,1) - uyy = dgg(2,2) - uzz = dgg(3,3) - uxy = (dgg(1,2) + dgg(2,1))*half - uxz = (dgg(1,3) + dgg(3,1))*half - uyz = (dgg(2,3) + dgg(3,2))*half - -! derivatives of S_I and Grad(S_I) - gx = ux*dffiu - gy = uy*dffiu - gz = uz*dffiu - gxx = dfuu*ux*ux + dffiu*uxx - gyy = dfuu*uy*uy + dffiu*uyy - gzz = dfuu*uz*uz + dffiu*uzz - gxy = dfuu*ux*uy + dffiu*uxy - gxz = dfuu*ux*uz + dffiu*uxz - gyz = dfuu*uy*uz + dffiu*uyz - - ggri(1,1,jk)=gxx - ggri(2,1,jk)=gxy - ggri(3,1,jk)=gxz - ggri(1,2,jk)=gxy - ggri(2,2,jk)=gyy - ggri(3,2,jk)=gyz - ggri(1,3,jk)=gxz - ggri(2,3,jk)=gyz - ggri(3,3,jk)=gzz - end do - - 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 -! output vector : dgg -! dff=(1,0,0) - use const_and_precisions, only : wp_ - implicit none -! arguments - real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3 - real(wp_), dimension(3), intent(out) :: dgg -! local variables - real(wp_) :: denom,aa1,aa2,aa3 - - aa1 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3)) - aa2 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3)) - aa3 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3)) - - denom = dxv1(1)*aa1 - dxv2(1)*aa2 + dxv3(1)*aa3 - - dgg(1) = aa1/denom - dgg(2) = -(dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3))/denom - dgg(3) = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2))/denom - end subroutine solg0 - - subroutine solg3(dxv1,dxv2,dxv3,dff,dgg) -! rhs "matrix" dff, result in dgg - use const_and_precisions, only : wp_ - implicit none -! arguments - real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3 - real(wp_), dimension(3,3), intent(in) :: dff - real(wp_), dimension(3,3), intent(out) :: dgg -! local variables - real(wp_) denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 - - a11 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3)) - a21 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3)) - a31 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3)) - - a12 = (dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3)) - a22 = (dxv1(1)*dxv3(3) - dxv3(1)*dxv1(3)) - a32 = (dxv1(1)*dxv2(3) - dxv2(1)*dxv1(3)) - - a13 = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2)) - a23 = (dxv1(1)*dxv3(2) - dxv3(1)*dxv1(2)) - a33 = (dxv1(1)*dxv2(2) - dxv2(1)*dxv1(2)) - - denom = dxv1(1)*a11 - dxv2(1)*a21 + dxv3(1)*a31 - - dgg(:,1) = ( dff(:,1)*a11 - dff(:,2)*a21 + dff(:,3)*a31)/denom - dgg(:,2) = (-dff(:,1)*a12 + dff(:,2)*a22 - dff(:,3)*a32)/denom - dgg(:,3) = ( dff(:,1)*a13 - dff(:,2)*a23 + dff(:,3)*a33)/denom - 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 - use gray_params, only : iequil - use equilibrium, only : psia,equinum_fpol,equinum_psi,equian,sgnbphi - use coreprofiles, only : density - implicit none -! arguments - real(wp_), dimension(3), intent(in) :: xv - real(wp_), intent(in) :: xgcn,bres - real(wp_), intent(out) :: psinv,dens,btot,xg,yg - real(wp_), dimension(3), intent(out) :: bv,derxg,deryg - real(wp_), dimension(3,3), intent(out) :: derbv -! local variables - integer :: jv - real(wp_) :: xx,yy,zz - real(wp_) :: b2tot,csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm - real(wp_), dimension(3) :: dbtot,bvc - real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv - real(wp_) :: brr,bphi,bzz,ajphi,dxgdpsi - real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddenspsin - - xg = zero - yg = 99._wp_ - psinv = -1._wp_ - dens = zero - btot = zero - ajphi = zero - derxg = zero - deryg = zero - bv = zero - derbv = zero - - if(iequil==0) return - - dbtot = zero - dbv = zero - dbvcdc = zero - dbvcdc = zero - dbvdc = zero - - xx = xv(1) - yy = xv(2) - zz = xv(3) - -! cylindrical coordinates - rr2 = xx**2 + yy**2 - rr = sqrt(rr2) - csphi = xx/rr - snphi = yy/rr - - bv(1) = -snphi*sgnbphi - bv(2) = csphi*sgnbphi - -! convert from cm to meters - zzm = 1.0e-2_wp_*zz - rrm = 1.0e-2_wp_*rr - - if(iequil==1) then - call equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, & - ddpsidrr,ddpsidzz,ddpsidrz) - else - call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz) - call equinum_fpol(psinv,fpolv,dfpv) - end if - -! compute yg and derivative - if(psinv < zero) then - bphi = fpolv/rrm - btot = abs(bphi) - yg = btot/bres - return - end if - -! compute xg and derivative - call density(psinv,dens,ddenspsin) - xg = xgcn*dens - dxgdpsi = xgcn*ddenspsin/psia - -! B = f(psi)/R e_phi+ grad psi x e_phi/R - bphi = fpolv/rrm - brr =-dpsidz/rrm - bzz = dpsidr/rrm - -! bvc(i) = B_i in cylindrical coordinates - bvc(1) = brr - bvc(2) = bphi - bvc(3) = bzz - -! bv(i) = B_i in cartesian coordinates - bv(1)=bvc(1)*csphi - bvc(2)*snphi - bv(2)=bvc(1)*snphi + bvc(2)*csphi - bv(3)=bvc(3) - -! dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) - dbvcdc(1,1) = -ddpsidrz/rrm - brr/rrm - dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm - dbvcdc(3,1) = ddpsidrr/rrm - bzz/rrm - dbvcdc(1,3) = -ddpsidzz/rrm - dbvcdc(2,3) = dfpv*dpsidz/rrm - dbvcdc(3,3) = ddpsidrz/rrm - -! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) - dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi - dbvdc(2,1) = dbvcdc(1,1)*snphi + dbvcdc(2,1)*csphi - dbvdc(3,1) = dbvcdc(3,1) - dbvdc(1,2) = -bv(2) - dbvdc(2,2) = bv(1) - dbvdc(3,2) = dbvcdc(3,2) - dbvdc(1,3) = dbvcdc(1,3)*csphi - dbvcdc(2,3)*snphi - dbvdc(2,3) = dbvcdc(1,3)*snphi + dbvcdc(2,3)*csphi - dbvdc(3,3) = dbvcdc(3,3) - - drrdx = csphi - drrdy = snphi - dphidx = -snphi/rrm - dphidy = csphi/rrm - -! dbv(iv,jv) = d Bcart(iv) / dxvcart(jv) - dbv(:,1) = drrdx*dbvdc(:,1) + dphidx*dbvdc(:,2) - dbv(:,2) = drrdy*dbvdc(:,1) + dphidy*dbvdc(:,2) - dbv(:,3) = dbvdc(:,3) - -! B magnitude and derivatives - b2tot = bv(1)**2 + bv(2)**2 + bv(3)**2 - btot = sqrt(b2tot) - - dbtot = (bv(1)*dbv(1,:) + bv(2)*dbv(2,:) + bv(3)*dbv(3,:))/btot - - yg = btot/Bres - -! convert spatial derivatives from dummy/m -> dummy/cm -! to be used in rhs - -! bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z - deryg = 1.0e-2_wp_*dbtot/Bres - bv = bv/btot - do jv=1,3 - derbv(:,jv) = 1.0e-2_wp_*(dbv(:,jv) - bv(:)*dbtot(jv))/btot - end do - - derxg(1) = 1.0e-2_wp_*drrdx*dpsidr*dxgdpsi - derxg(2) = 1.0e-2_wp_*drrdy*dpsidr*dxgdpsi - derxg(3) = 1.0e-2_wp_*dpsidz *dxgdpsi - -! current density computation in Ampere/m^2, ccj==1/mu_0 - ajphi = ccj*(dbvcdc(1,3) - dbvcdc(3,1)) -! ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3)) -! ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2)) - - 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 - use gray_params, only : idst,igrad - implicit none -! arguments - real(wp_), intent(in) :: xg,yg,gr2,sox - real(wp_), intent(out) :: anpl,anpr,ddr,ddi,derdnm,dersdst - real(wp_), dimension(3), intent(in) :: anv,bv,derxg,deryg - real(wp_), dimension(3), intent(in) :: dgr2,dgr - real(wp_), dimension(3,3), intent(in) :: ddgr,derbv - real(wp_), dimension(6), intent(out) :: dery -! local variables - integer :: iv - real(wp_) :: yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl,an2,an2s - real(wp_) :: dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel - real(wp_) :: derdom,dfdiadnpl,dfdiadxg,dfdiadyg,fdia,bdotgr !,vgm - real(wp_), dimension(3) :: derdxv,danpldxv,derdnv,dbgr !,vgv - - an2 = anv(1)*anv(1) + anv(2)*anv(2) + anv(3)*anv(3) - anpl = anv(1)*bv(1) + anv(2)*bv(2) + anv(3)*bv(3) - - anpl2 = anpl**2 - dnl = one - anpl2 - anpr2 = max(an2-anpl2,zero) - anpr = sqrt(anpr2) - yg2 = yg**2 - - an2s = one - dan2sdxg = zero - dan2sdyg = zero - dan2sdnpl = zero - del = zero - fdia = zero - dfdiadnpl = zero - dfdiadxg = zero - dfdiadyg = zero - - duh = one - xg - yg2 - if(xg > zero) then - del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2) - an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh - - dan2sdxg = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2 & - + sox*xg*anpl2/(del*duh) - one - dan2sdyg = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2 & - + two*sox*xg*(one - xg)*anpl2/(yg*del*duh) - dan2sdnpl = - xg*yg2*anpl/duh & - - sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh) - - if(igrad > 0) then - ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) & - - yg2*dnl**3)/yg2/del**3 - fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh - derdel = two*(one - xg)*anpl2*(one + 3.0_wp_*anpl2**2) & - - dnl**2*(one + 3.0_wp_*anpl2)*yg2 - derdel = 4.0_wp_*derdel/(yg*del)**5 - ddelnpl2y = two*(one - xg)*derdel - ddelnpl2x = yg*derdel - dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2) & - /(yg2*del**5) - dfdiadxg = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2) & - *ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2) - dfdiadyg = - two*yg*xg*(one - xg)/duh**2 & - - sox*xg*yg*(two*(one - xg)*ddelnpl2 & - + yg*duh*ddelnpl2y)/(two*duh**2) - end if - end if - - bdotgr = bv(1)*dgr(1) + bv(2)*dgr(2) + bv(3)*dgr(3) - do iv=1,3 - dbgr(iv) = dgr(1)*derbv(1,iv) + bv(1)*ddgr(1,iv) & - + dgr(2)*derbv(2,iv) + bv(2)*ddgr(2,iv) & - + dgr(3)*derbv(3,iv) + bv(3)*ddgr(3,iv) - danpldxv(iv) = anv(1)*derbv(1,iv) + anv(2)*derbv(2,iv) + anv(3)*derbv(3,iv) - end do - - derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl + & - igrad*dgr2) & - + fdia*bdotgr*dbgr + half*bdotgr**2 & - *(derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl) - derdnv = two*anv + (half*bdotgr**2*dfdiadnpl - dan2sdnpl)*bv - - derdnm = sqrt(derdnv(1)**2 + derdnv(2)**2 + derdnv(3)**2) - - derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl & - + two*igrad*gr2 - bdotgr**2*(fdia + xg*dfdiadxg & - + half*yg*dfdiadyg & - + half*anpl*dfdiadnpl) - - if (idst == 0) then -! integration variable: s - denom = derdnm - else if (idst == 1) then -! integration variable: c*t - denom = -derdom - else -! integration variable: Sr - denom = anv(1)*derdnv(1) + anv(2)*derdnv(2) + anv(3)*derdnv(3) - end if - -! coefficient for integration in s -! ds/dst, where st is the integration variable - dersdst = derdnm/denom - -! rhs vector - dery(1:3) = derdnv(:)/denom - dery(4:6) = -derdxv(:)/denom - -! vgv : ~ group velocity -! vgm=0 -! do iv=1,3 -! vgv(iv)=-derdnv(iv)/derdom -! vgm=vgm+vgv(iv)**2 -! end do -! vgm=sqrt(vgm) - -! ddr : dispersion relation (real part) -! ddi : dispersion relation (imaginary part) - ddr = an2 - an2s - igrad*(gr2 - half*bdotgr**2*fdia) - ddi = derdnv(1)*dgr(1) + derdnv(2)*dgr(2) + derdnv(3)*dgr(3) - - 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_ - use gray_params, only : iwarm,ilarm,ieccd,imx - use coreprofiles, only : fzeff - use equilibrium, only : sgnbphi - use dispersion, only : harmnumber, warmdisp - use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl - use errcodes, only : palph - use magsurf_data, only : fluxval - implicit none -! arguments - real(wp_),intent(in) ::psinv,ak0,bres - real(wp_),intent(in) :: xg,yg,tekev,dens,anpl,anpr,derdnm,sox - real(wp_),intent(out) :: anprre,anprim,alpha,didp - integer, intent(out) :: nhmin,nhmax,iokhawa - integer, intent(out) :: ierr -! local constants - real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ -! local variables - real(wp_) :: rbavi,rrii,rhop - integer :: lrm,ithn,ierrcd - real(wp_) :: amu,ratiovgr,rbn,rbx - real(wp_) :: zeff,cst2,bmxi,bmni,fci - real(wp_), dimension(:), pointer :: eccdpar=>null() - real(wp_) :: effjcd,effjcdav,akim,btot - complex(wp_) :: ex,ey,ez - - alpha=zero - anprim=zero - anprre=zero - didp=zero - nhmin=0 - nhmax=0 - iokhawa=0 - ierr=0 - - if(tekev>zero) then -! absorption computation - amu=mc2/tekev - call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) - if(nhmin.gt.0) then - lrm=max(ilarm,nhmax) - call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim, & - iwarm,imx,ex,ey,ez) - akim=ak0*anprim - ratiovgr=2.0_wp_*anpr/derdnm!*vgm - alpha=2.0_wp_*akim*ratiovgr - if(alpha: effjcdav [A m/W ] - if(ieccd>0) then -! current drive computation - zeff=fzeff(psinv) - ithn=1 - if(lrm>nhmin) ithn=2 - rhop=sqrt(psinv) - call fluxval(rhop,rri=rrii,rbav=rbavi,bmn=bmni,bmx=bmxi,fc=fci) - btot=yg*bres - rbn=btot/bmni - rbx=btot/bmxi - - select case(ieccd) - case(1) -! cohen model - call setcdcoeff(zeff,rbn,rbx,cst2,eccdpar) - call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierrcd) - case(2) -! no trapping - call setcdcoeff(zeff,cst2,eccdpar) - call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierrcd) - case default -! neoclassical model - call setcdcoeff(zeff,rbx,fci,amu,rhop,cst2,eccdpar) - call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd) - end select - ierr=ierr+ierrcd - if(associated(eccdpar)) deallocate(eccdpar) - effjcdav=rbavi*effjcd - didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii) - end if - end if - end if - end subroutine alpha_effj - - - - subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0) - use const_and_precisions, only : wp_,degree,zero,one,half,im - use beamdata, only : nray,nrayth - use equilibrium, only : bfield - use gray_params, only : ipol - use polarization, only : pol_limit, polellipse, stokes_ce, stokes_ell - implicit none -! arguments - real(wp_), dimension(6,nray), intent(in) :: ywrk0 - real(wp_), intent(in) :: sox,bres - real(wp_), intent(inout) :: psipol0, chipol0 - complex(wp_), dimension(nray), intent(out) :: ext0, eyt0 -! local variables - integer :: j,k,jk - real(wp_), dimension(3) :: xmv, anv, bv - real(wp_) :: rm, csphi, snphi, bphi, br, bz, qq, uu, vv, deltapol - - j=1 - k=0 - do jk=1,nray - k=k+1 - if(jk == 2 .or. k > nrayth) then - j=j+1 - k=1 - end if - - if(ipol == 0) then - xmv=ywrk0(1:3,jk)*0.01_wp_ ! convert from cm to m - anv=ywrk0(4:6,jk) - rm=sqrt(xmv(1)**2+xmv(2)**2) - csphi=xmv(1)/rm - snphi=xmv(2)/rm - call bfield(rm,xmv(3),bphi,br,bz) -! bv(i) = B_i in cartesian coordinates - bv(1)=br*csphi-bphi*snphi - bv(2)=br*snphi+bphi*csphi - bv(3)=bz - call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk)) - - if (jk == 1) then - call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv) - call polellipse(qq,uu,vv,psipol0,chipol0) - psipol0=psipol0/degree ! convert from rad to degree - chipol0=chipol0/degree - end if - else - call stokes_ell(chipol0*degree,psipol0*degree,qq,uu,vv) - if(qq**2 < one) then - deltapol=asin(vv/sqrt(one - qq**2)) - ext0(jk)= sqrt(half*(one + qq)) - eyt0(jk)= sqrt(half*(one - qq))*exp(-im*deltapol) - else - ext0(jk)= one - eyt0(jk)= zero - end if - endif - 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(strheader) - use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm - implicit none -! arguments - character(len=*), dimension(:), intent(in) :: strheader -! local variables - integer :: i,l - - l=size(strheader) - do i=1,l - write(uprj0,'(1x,a)') strheader(i) - write(uprj0+1,'(1x,a)') strheader(i) - write(uwbm,'(1x,a)') strheader(i) - write(udisp,'(1x,a)') strheader(i) - write(ucenr,'(1x,a)') strheader(i) - write(uoutr,'(1x,a)') strheader(i) - write(upec,'(1x,a)') strheader(i) - write(usumm,'(1x,a)') strheader(i) - end do - write(uprj0,'(1x,a)') '#sst j k xt yt zt rt' - write(uprj0+1,'(1x,a)') '#sst j k xt yt zt rt' - write(uwbm,'(1x,a)') '#sst w1 w2' - write(udisp,'(1x,a)') '#sst Dr_Nr Di_Nr' - write(ucenr,'(1x,a)') '#sst R z phi psin rhot ne Te Btot Bx By Bx Nperp Npl '// & - 'Nx Ny Nz ki alpha tau Pt dIds nhmax iohkw index_rt ddr' - write(uoutr,'(1x,a)') '#i k sst x y R z psin tau Npl alpha index_rt' - write(upec,'(1x,a)') '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt' - write(usumm,'(1x,a)') '#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 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_maps(bres,xgcn,r0,anpl0) - use const_and_precisions, only : wp_ - use gray_params, only : iequil - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, & - equinum_fpol, nq - use coreprofiles, only : density, temp - use units, only : umaps - implicit none -! arguments - real(wp_), intent(in) :: bres,xgcn,r0,anpl0 -! local variables - integer :: j,k - real(wp_) :: dr,dz,zk,rj,bphi,br,bz,btot,psin,ne,dne,te,xg,yg,anpl - real(wp_), dimension(nq) :: r, z - - - dr = (rmxm-rmnm)/(nq-1) - dz = (zmxm-zmnm)/(nq-1) - do j=1,nq - r(j) = rmnm + dr*(j-1) - z(j) = zmnm + dz*(j-1) - end do - - write(umaps,*)'#R z psin Br Bphi Bz Btot ne Te X Y Npl' - do j=1,nq - rj=r(j) - anpl=anpl0*r0/rj - do k=1,nq - zk=z(k) - if (iequil < 2) then - call equian(rj,zk,psinv=psin,fpolv=bphi,dpsidr=bz,dpsidz=br) - else - call equinum_psi(rj,zk,psinv=psin,dpsidr=bz,dpsidz=br) - call equinum_fpol(psin,fpolv=bphi) - end if - br = -br/rj - bphi = bphi/rj - bz = bz/rj - btot = sqrt(br**2+bphi**2+bz**2) - yg = btot/bres - te = temp(psin) - call density(psin,ne,dne) - xg = xgcn*ne - write(umaps,'(12(x,e12.5))') rj,zk,psin,br,bphi,bz,btot,ne,te,xg,yg,anpl - enddo - write(umaps,*) - enddo - - end subroutine print_maps - - - - - 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 :: 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,bv,ak0,anpl,anpr,anv, & - 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,jkray1 - use units, only : ucenr,uoutr,udisp - implicit none -! arguments - integer, intent(in) :: i,jk,nhf,iokhawa,index_rt - real(wp_), dimension(3), intent(in) :: xv,bv,anv - real(wp_), intent(in) :: st,qj,psinv,btot,ak0,anpl,anpr,anprim - real(wp_), intent(in) :: dens,tekev,alpha,tau,dids,ddr,ddi -! local variables - real(wp_) :: stm,xxm,yym,zzm,rrm,phideg,rhot,akim,pt,didsn - integer :: k - - stm=st*1.0e-2_wp_ - xxm=xv(1)*1.0e-2_wp_ - yym=xv(2)*1.0e-2_wp_ - zzm=xv(3)*1.0e-2_wp_ - rrm=sqrt(xxm**2 + yym**2) - -! 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 - rhot=frhotor(sqrt(psinv)) - else - rhot=1.0_wp_ - end if - akim=anprim*ak0*1.0e2_wp_ - pt=exp(-tau) - didsn=dids*1.0e2_wp_/qj - - write(ucenr,'(22(1x,e16.8e3),3i5,1x,e16.8e3)') stm,rrm,zzm,phideg, & - psinv,rhot,dens,tekev,btot,bv,anpr,anpl,anv,akim,alpha,tau,pt,didsn, & - nhf,iokhawa,index_rt,ddr - end if - -! print conservation of dispersion relation - if(jk==nray) write(udisp,'(30(1x,e16.8e3))') st,ddr,ddi - -! print outer trajectories - if(mod(i,istpl0)==0) then - 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 sumcore