From d8fc72dd3e9e0f0c11111f592ea832db7ac1f03f Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Wed, 13 Sep 2017 09:31:30 +0000 Subject: [PATCH] re-added routines to sum profiles (removed by mistake in rev 176) --- Makefile.sum | 75 ++ src/main-sum.f90 | 195 +++++ src/sumcore.f90 | 1907 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 2177 insertions(+) create mode 100644 Makefile.sum create mode 100644 src/main-sum.f90 create mode 100644 src/sumcore.f90 diff --git a/Makefile.sum b/Makefile.sum new file mode 100644 index 0000000..c9d8f58 --- /dev/null +++ b/Makefile.sum @@ -0,0 +1,75 @@ +# Executable name +EXE=sumgray + +# Objects list +MAINOBJ=main-sum.o +OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \ + dierckx.o dispersion.o eccd.o eierf.o errcodes.o sumcore.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 +vpath %.f src + +# Fortran compiler name and flags +FC=gfortran +FFLAGS=-O3 +#FFLAGS=-Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check + +DIRECTIVES = -DREVISION="'$(shell svnversion src)'" + +all: $(EXE) + +# Build executable from object files +$(EXE): $(MAINOBJ) $(OTHOBJ) + $(FC) $(FFLAGS) -o $@ $^ + +# Dependencies on modules +main-sum.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \ + sumcore.o gray_params.o reflections.o +sumcore.o: const_and_precisions.o beamdata.o beams.o coreprofiles.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 +coreprofiles.o: const_and_precisions.o dierckx.o gray_params.o simplespline.o \ + utils.o +dierckx.o: const_and_precisions.o +dispersion.o: const_and_precisions.o eierf.o errcodes.o math.o quadpack.o +eccd.o: const_and_precisions.o conical.o dierckx.o errcodes.o magsurf_data.o \ + numint.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 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 units.o utils.o +math.o: const_and_precisions.o +minpack.o: const_and_precisions.o +numint.o: const_and_precisions.o +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 limiter.o utils.o +simplespline.o: const_and_precisions.o +utils.o: const_and_precisions.o + +# General object compilation command +%.o: %.f90 + $(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< + +.PHONY: clean install +# Remove output files +clean: + rm -rf *.o *.mod $(EXE) + +install: + @if [ -f $(EXE) ]; then \ + cp $(EXE) ~/bin/; \ + else \ + echo File $(EXE) does not exist. Run \'make\' first; \ + fi diff --git a/src/main-sum.f90 b/src/main-sum.f90 new file mode 100644 index 0000000..2bd9407 --- /dev/null +++ b/src/main-sum.f90 @@ -0,0 +1,195 @@ +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, rv, zv, fpol, qpsi) +! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0) + psia = sign(one,qpsi(2)*fpol(1)) + else + call read_eqdsk(eqp%filenm, rv,zv,psin, psia, psinr,fpol,qpsi, rvac, & + rax,zax, rbnd,zbnd, rlim,zlim, eqp%ipsinorm,eqp%idesc,eqp%ifreefmt) + call change_cocos(psia, fpol, qpsi, eqp%icocos, 3) + end if +! re-scale B/I and/or force signs. If sgn=0 on input, set to fpol/-psia signs on output + call eq_scal(psia, fpol, eqp%sgni, eqp%sgnb, eqp%factb) +! ??? 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=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/sumcore.f90 b/src/sumcore.f90 new file mode 100644 index 0000000..f58ac57 --- /dev/null +++ b/src/sumcore.f90 @@ -0,0 +1,1907 @@ +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