From 27f1793f14e5c95c0e1a067cf3aa6c04fcad61ce Mon Sep 17 00:00:00 2001 From: Daniele Micheletti Date: Tue, 26 Mar 2019 14:21:22 +0000 Subject: [PATCH] trunk: added multiple passes calculation; added multipass module; graycore/gray_main: subroutine rewritten; equilibrium/read_equil_an: limiter coordinates read from equilibrium file for ipass>1; err_codes/check_err: istop=1 only for alpha<0; pec/pec_tab: fixed index assignements in loops; polarization/pol_limit: fixed ext,eyt normalization --- Makefile | 10 +- src/equilibrium.f90 | 25 +- src/errcodes.f90 | 12 +- src/graycore.f90 | 680 +++++++++++++++++++++++++++++++------------ src/main-sum.f90 | 4 +- src/main.f90 | 4 +- src/multipass.f90 | 280 ++++++++++++++++++ src/pec.f90 | 19 +- src/polarization.f90 | 3 + 9 files changed, 832 insertions(+), 205 deletions(-) create mode 100755 src/multipass.f90 diff --git a/Makefile b/Makefile index 3aac4ad..9146b39 100644 --- a/Makefile +++ b/Makefile @@ -5,8 +5,8 @@ EXE=gray MAINOBJ=main.o OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \ dierckx.o dispersion.o eccd.o eierf.o errcodes.o graycore.o \ - gray_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 + gray_params.o equilibrium.o limiter.o magsurf_data.o math.o minpack.o multipass.o \ + numint.o pec.o polarization.o quadpack.o reflections.o simplespline.o units.o utils.o # Alternative search paths vpath %.f90 src @@ -15,7 +15,7 @@ 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 +# FFLAGS=-O0 -Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check DIRECTIVES = -DREVISION="'$(shell svnversion src)'" @@ -30,7 +30,7 @@ main.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \ graycore.o gray_params.o reflections.o graycore.o: const_and_precisions.o beamdata.o beams.o coreprofiles.o \ dispersion.o eccd.o equilibrium.o errcodes.o gray_params.o \ - pec.o polarization.o limiter.o units.o utils.o + pec.o polarization.o limiter.o units.o utils.o reflections.o multipass.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 @@ -49,6 +49,8 @@ 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 +multipass.o: const_and_precisions.o beamdata.o equilibrium.o gray_params.o \ + polarization.o reflections.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 diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 8e96d80..d02f8cc 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -13,7 +13,7 @@ module equilibrium ! === 2d spline psi(r,z), normalization and derivatives ========== integer, save :: nsr, nsz - real(wp_), save :: psia, psiant, psinop + real(wp_), save :: psia, psiant, psinop, psnbd real(wp_), dimension(:), allocatable, save :: tr,tz real(wp_), dimension(:), allocatable, save :: cceq, cceq01, cceq10, & cceq20, cceq02, cceq11 @@ -151,15 +151,16 @@ contains - subroutine read_equil_an(filenm,rv,zv,fpol,q,unit) + subroutine read_equil_an(filenm,ipass,rv,zv,fpol,q,rlim,zlim,unit) use utils, only : get_free_unit implicit none ! arguments character(len=*), intent(in) :: filenm + integer, intent(in) :: ipass integer, optional, intent(in) :: unit - real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q + real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q,rlim,zlim ! local variables - integer :: u + integer :: i, u, nlim real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen if (present(unit)) then @@ -191,6 +192,22 @@ contains q(1)=q0 q(2)=qa q(3)=alq + + if(ipass.ge.2) then + + ! get size of boundary and limiter arrays and allocate them + read (u,*) nlim + if (allocated(rlim)) deallocate(rlim) + if (allocated(zlim)) deallocate(zlim) + + ! store boundary and limiter data + if(nlim>0) then + allocate(rlim(nlim),zlim(nlim)) + read(u,*) (rlim(i),zlim(i),i=1,nlim) + end if + + end if + close(u) end subroutine read_equil_an diff --git a/src/errcodes.f90 b/src/errcodes.f90 index b50c583..a1a7f7f 100644 --- a/src/errcodes.f90 +++ b/src/errcodes.f90 @@ -13,12 +13,16 @@ contains ! arguments integer, intent(in) :: ierr integer, intent(inout) :: istop - if(ibits(ierr,pnpl, lnpl )>1 .or. & ! N// too large - ibits(ierr,palph,lalph)==1) then ! alpha < 0 - istop = 1 +! if(ibits(ierr,pnpl, lnpl )>1 .or. & ! N// too large +! ibits(ierr,palph,lalph)==1) then ! alpha < 0 +! istop = 1 ! else ! istop = 0 - end if +! end if + +! if(ibits(ierr,pnpl, lnpl )>1) istop = 1 ! N// too large + if(ibits(ierr,palph,lalph)==1) istop = 1 ! alpha < 0 + end subroutine check_err diff --git a/src/graycore.f90 b/src/graycore.f90 index 1e3139e..7884d95 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -8,16 +8,16 @@ contains eqp,psrad,terad,derad,zfc,prfp, rlim,zlim, & p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0, & psipol0,chipol0, dpdv,jcd,pabs,icd, outp,rtrp,hcdp,ierr, rhout) - use const_and_precisions, only : zero, one, degree + use const_and_precisions, only : zero, one, degree, comp_tiny 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 + iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl, ipass 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 + zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q, psnbd 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 @@ -25,6 +25,10 @@ contains rhop_tab, rhot_tab use limiter, only : set_lim, unset_lim use utils, only : vmaxmin + use reflections, only : inside + use multipass, only : alloc_multipass, dealloc_multipass, initbeam, & + initmultipass, turnoffray, plasma_in, plasma_out, wall_out + implicit none ! arguments type(eqparam_type), intent(in) :: eqp @@ -49,38 +53,50 @@ contains 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_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) + character, dimension(2), parameter :: mode=(/'O','X'/) + + real(wp_) :: sox,ak0,bres,xgcn,xg,yg,rrm,zzm,alpha,didp,anpl,anpr,anprim,anprre + real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx + real(wp_) :: pabs_beam,icd_beam,cpl_beam1,cpl_beam2,cpl_cbeam1,cpl_cbeam2 + real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl,cpl0 real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg 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 + integer :: i,j,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt + integer :: ip,ib,iopmin,ipar,iO + integer :: igrad_b,ipol,istop_pass,nbeam_pass,nlim + logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff - 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() - - real(wp_), dimension(:), allocatable :: jphi,pins,currins + real(wp_), dimension(:,:,:), pointer :: yynext=>null(),yypnext=>null() + real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null() + real(wp_), dimension(:,:), pointer :: taus=>null(),stnext=>null(), & + yw0=>null(),ypw0=>null(),cpls=>null() + real(wp_), dimension(:), pointer :: p0ray=>null(),tau0=>null(),alphaabs0=>null(), & + dids0=>null(),ccci0=>null(),tau1=>null(),etau1=>null(),cpl1=>null(),lgcpl1=>null() + real(wp_), dimension(:), pointer :: p0jk=>null() + real(wp_), dimension(:), pointer :: jphi_beam=>null(),pins_beam=>null(), & + currins_beam=>null(), dpdv_beam=>null(),jcd_beam=>null(),stv=>null(), & + psipv=>null(),chipv=>null() + complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null() + integer, dimension(:), pointer :: iiv=>null(),iop=>null(),iow=>null() + logical, dimension(:), pointer :: iwait=>null() + logical, dimension(:,:), pointer :: iroff=>null() ! parameters log in file headers character(len=headw), dimension(headl) :: strheader type(antctrl_type) :: antp real(wp_) :: rwall -! ======= set environment BEGIN ====== +! ======== set environment BEGIN ======== call set_codepar(eqp,prfp,outp,rtrp,hcdp) call set_lim(rlim,zlim) + nlim = size(zlim) if(iequil<2) then call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3)) @@ -104,10 +120,16 @@ contains tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) if(iwarm > 1) call expinit + + call pec_init(ipec,rhout) + nnd=size(rhop_tab) + call alloc_multipass(nnd,iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & + stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & + pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) -! ======= set environment END ====== +! ========= set environment END ========= -! ======= pre-proc prints BEGIN ====== +! ======== pre-proc prints BEGIN ======== antp%alpha=alpha0 antp%beta=beta0 antp%power=p0 @@ -118,8 +140,16 @@ contains antp%ibeam=0 antp%filenm='' rwall=0._wp_ + psnbd=prfp%psnbnd ! plasma boundary + ipol=rtrp%ipol + + pabs=zero ! gray_main output + icd=zero + dpdv=zero + jcd=zero + call print_params(rtrp,hcdp,antp,eqp,rwall,prfp,outp,strheader) - call print_headers(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 @@ -130,161 +160,435 @@ contains ! 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 ====== + call print_maps(bres,xgcn,0.01_wp_*sqrt(xv0(1)**2+xv0(2)**2),sin(beta0*degree)) -! ======= main loop BEGIN ====== - iox=iox0 - sox=-1.0_wp_ - if(iox==2) sox=1.0_wp_ - call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) - call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri) +! ========= pre-proc prints END ========= + +! =========== main loop BEGIN =========== + call initmultipass(ipol,iox0,iroff,yynext,yypnext,yw0,ypw0, & + stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) + + if(ipol.eq.0) then + if(iox0.eq.2) then ! only X mode on 1st pass + cpl0 = (/zero,one/) + else ! only O mode on 1st pass + cpl0 = (/one,zero/) + end if + end if + + sox=one ! mode inverted for each beam + iox=2 ! start with O: sox=-1, iox=1 psipol=psipol0 chipol=chipol0 - call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) call pweight(p0,p0jk) - st0 = zero - if(nray>1) call print_projxyzt(st0,yw,0) ! iproj=0 ==> nfilp=8 - - somein = .false. ! becomes true if at least part of the beam enters the plasma -! beam/ray propagation - do i=1,nstep - -! advance one step with "frozen" grad(S_I) - st=i*dst+st0 - do jk=1,nray - call rkstep(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk)) - end do - -! update position and grad - if(igrad == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) - - allout = .true. ! becomes false if at least part of the beam is inside the plsama - ierr = 0 - istop = 0 - do jk=1,nray -! compute derivatives with updated gradient and local plasma values - xv = yw(1:3,jk) - anv = yw(4:6,jk) - call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), & - psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn) -! update global error code and print message - if (ierrn/=0) then - ierr = ior(ierr,ierrn) - call print_errn(ierrn,i,anpl) - end if - - zzm = xv(3)*0.01_wp_ - ins_pl = (psinv>=zero .and. psinv=zbinf .and. zzm<=zbsup) -! test if the beam is completely out of the plsama - allout = allout .and. .not.ins_pl -! test if at least part of the beam has entered the plsama - somein = somein .or. ins_pl - -! compute ECRH&CD - if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. tau0(jk)<=taucr) then - tekev=temp(psinv) - call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & - sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd) - if (ierrhcd/=0) then - ierr = ior(ierr,ierrhcd) - call print_errhcd(ierrhcd,i,anprre,anprim,alpha) - end if - else - tekev=zero - alpha=zero - didp=zero - anprim=zero - anprre=anpr - nharm=0 - nhf=0 - iokhawa=0 - end if - if(nharm>0) iiv(jk)=i - - psjki(jk,i) = psinv -! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s) - tau=tau0(jk)+0.5_wp_*(alphaabs0(jk)+alpha)*dersdst*dst - pow=p0jk(jk)*exp(-tau) !*exp(-tau1v(jk)) - ppabs(jk,i)=p0jk(jk)-pow - - dids=didp*pow*alpha - ccci(jk,i)=ccci0(jk)+0.5_wp_*(dids0(jk)+dids)*dersdst*dst - - tau0(jk)=tau - alphaabs0(jk)=alpha - dids0(jk)=dids - ccci0(jk)=ccci(jk,i) - - call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,bv,ak0,anpl,anpr, & - anv,anprim,dens,tekev,alpha,tau,dids,nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) - - end do - -! print ray positions for j=nrayr in local reference system - if (mod(i,istpr0) == 0) then - if(nray > 1) call print_projxyzt(st,yw,0) + nbeam_pass=1 ! max n of beam per pass + index_rt=0 ! global beam index: 1,O 2,X 1st pass + ! / \ / \ + do ip=1,ipass ! 3,O 4,X 5,O 6,X 2nd pass + + pabs_pass = zero + icd_pass = zero + istop_pass = 0 ! stop flag for current pass + nbeam_pass = 2*nbeam_pass ! max n of beams in current pass + + if(ip.gt.1) then + du1 = zero ! igrad=0 from 2nd pass + gri = zero + ggri = zero + if(ip.eq.ipass) cpl = (/zero,zero/) ! no successive passes end if + +! =========== beam loop BEGIN =========== + do ib=1,nbeam_pass + + sox = -sox ! invert mode + iox = 3-iox ! O-mode at ip=1,ib=1 + + index_rt = index_rt +1 + iO = 2*index_rt +1 ! * index_rt of O-mode derived ray (iX=iO+1) -! check for any error code and stop if necessary - call check_err(ierr,istop) -! test whether further trajectory integration is unnecessary - call vmaxmin(tau0,nray,taumn,taumx) - if ((taumn > taucr) .or. (somein .and. allout)) istop = 1 + call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, & + pins_beam,currins_beam,dpdv_beam,jcd_beam) + + if(iboff) then ! no propagation for current beam + istop_pass = istop_pass +1 ! * +1 non propagating beam + cycle + end if + + call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) + call print_headers((/' '/),index_rt) + + if(ip.eq.1) then ! first pass + igrad_b = igrad + + tau1 = zero + etau1 = one + cpl1 = one + lgcpl1 = zero + p0ray = p0jk ! * initial beam power - if(istop == 1) exit + call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri,index_rt) ! * initial conditions + call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization + + do jk=1,nray + if(iwait(jk)) cycle + zzm = yw(3,jk)*0.01_wp_ + rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ + + if(inside(rlim,zlim,nlim,rrm,zzm)) then ! * start propagation in/outside vessel? + iow(jk) = 1 ! + inside + else + iow(jk) = 0 ! + outside + end if + end do + + else ! successive passes + ipar = (index_rt+1)/2-1 ! * parent beam index + yw = yynext(:,:,ipar) ! * starting coordinates from + ypw = yypnext(:,:,ipar) ! parent beam last step + stv = stnext(:,ipar) ! * starting step from parent beam last step + iow = 1 ! * start propagation inside vessel + + tau1 = taus(:,index_rt) ! * tau from previous passes + etau1 = exp(-tau1) + cpl1 = cpls(:,index_rt) ! * coupling from previous passes + lgcpl1 = -log(cpl1) + p0ray = p0jk * etau1 * cpl1 ! * initial beam power + end if + iop = 0 ! start propagation outside plasma + + if(nray>1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,0) ! iproj=0 ==> nfilp=8 + +! ======= propagation loop BEGIN ======= + do i=1,nstep + + ! advance one step with "frozen" grad(S_I) + do jk=1,nray + if(iwait(jk)) cycle ! jk ray is waiting for next pass + stv(jk) = stv(jk) + dst ! current ray step + call rkstep(sox,bres,xgcn,yw(:,jk),ypw(:,jk),gri(:,jk),ggri(:,:,jk),igrad_b) + end do + ! update position and grad + if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) + + ierr = 0 + istop = 0 ! stop flag for current beam + iopmin = 10 + +! =========== rays loop BEGIN =========== + do jk=1,nray + if(iwait(jk)) cycle ! jk ray is waiting for next pass + + ! compute derivatives with updated gradient and local plasma values + xv = yw(1:3,jk) + anv = yw(4:6,jk) + call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), & + psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm, & + ierrn,igrad_b) + ! update global error code and print message + if(ierrn/=0) then + ierr = ior(ierr,ierrn) + call print_errn(ierrn,i,anpl) + end if + + ! check entrance/exit plasma/wall + zzm = xv(3)*0.01_wp_ + rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_ + + ins_pl = (psinv>=zero .and. psinv=zbinf .and. zzm<=zbsup ! in/out plasma? + ins_wl = inside(rlim,zlim,nlim,rrm,zzm) ! in/out vessel? + ent_pl = (mod(iop(jk),2).eq.0 .and. ins_pl) ! enter plasma + ext_pl = (mod(iop(jk),2).eq.1 .and. .not.ins_pl) ! exit plasma + ent_wl = (mod(iow(jk),2).eq.0 .and. ins_wl) ! enter vessel + ext_wl = (mod(iow(jk),2).eq.1 .and. .not.ins_wl) ! exit vessel + + if(ent_pl) then ! ray enters plasma + call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt) + + if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray never entered in plasma) => continue current pass + + if(ipol.eq.0) then ! + IF single mode propagation + cpl = cpl0 + p0ray(jk) = p0ray(jk)*cpl(iox) + else if(cpl(iox).lt.etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays + call turnoffray(jk,ip+1,2*ib+2-iox,iroff) + iwait(jk) = .true. ! . stop advancement and H&CD computation for current ray + if(cpl(iox).le.comp_tiny) cpl(iox)=etaucr + else ! + ELSE assign coupled power to current ray + p0ray(jk) = p0ray(jk)*cpl(iox) + end if + cpls(jk,index_rt) = cpl(iox) + + if(jk.eq.1) then + write(*,*) + write(*,'("1st pass coupling (central ray, ",a1,"-mode)",f9.4)'), & + mode(iox),cpl(iox) + psipv(index_rt) = psipol ! + polarization angles at plasma boundary for central ray + chipv(index_rt) = chipol + end if + + else if(iop(jk).gt.2) then ! * 2nd entrance on 1st pass / entrance on 2nd+ pass => end of current pass for ray jk + igrad_b = 0 ! + switch to ray-tracing + iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray + + if(ip.lt.ipass) then ! not last pass + yynext(:,jk,index_rt) = yw0(:,jk) ! + copy starting coordinates + yypnext(:,jk,index_rt) = ypw0(:,jk) ! for next pass from last step + stnext(jk,index_rt) = stv(jk) - dst ! + starting step for next pass = last step + + if(cpl(1).lt.etaucr) then ! + low coupled power for O-mode => de-activate derived rays + call turnoffray(jk,ip+1,2*ib-1,iroff) + if(cpl(1).le.comp_tiny) cpl(1)=etaucr + end if + if(cpl(2).lt.etaucr) then ! + low coupled power for X-mode => de-activate derived rays + call turnoffray(jk,ip+1,2*ib,iroff) + if(cpl(2).le.comp_tiny) cpl(2)=etaucr + end if + + taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! + starting tau for next O-mode pass + cpls(jk,iO) = cpl1(jk) * cpl(1) ! + cumulative coupling for next O-mode pass + cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! + cumulative coupling for next X-mode pass + + if(jk.eq.1) then ! + polarization angles at plasma boundary for central ray + psipv(iO:iO+1) = psipol + chipv(iO:iO+1) = chipol + end if + else + cpl = (/zero,zero/) + end if + end if + + else if(ext_pl) then ! ray exits plasma + call plasma_out(jk,xv,anv,bres,sox,iop,ext,eyt) + end if + + if(ent_wl) then ! ray enters vessel + iow(jk)=iow(jk)+1 + + else if(ext_wl) then ! ray exit vessel + call wall_out(jk,ins_pl,xv,anv,bres,sox,psipol,chipol,iow,iop,ext,eyt) + yw(:,jk) = (/xv,anv/) ! * updated coordinates (reflected) + igrad_b = 0 ! * switch to ray-tracing + + call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), & + psinv,dens,btot,bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm, & + ierrn,igrad_b) ! * update derivatives after reflection + if(ierrn/=0) then ! * update global error code and print message + ierr = ior(ierr,ierrn) + call print_errn(ierrn,i,anpl) + end if + + if(jk.eq.1 .and. ip.eq.1) then ! * 1st pass, polarization angles at reflection for central ray + psipv(index_rt) = psipol + chipv(index_rt) = chipol + end if + + if(ins_pl) then ! * plasma-wall overlapping => wall crossing=plasma crossing => end of current pass + iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray + + if(ip.lt.ipass) then ! + not last pass + yynext(:,jk,index_rt) = (/xv,anv/) ! . starting coordinates + yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point + stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection + + call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt) ! . ray re-enters plasma after reflection + + if(cpl(1).lt.etaucr) then ! . low coupled power for O-mode? => de-activate derived rays + call turnoffray(jk,ip+1,2*ib-1,iroff) + if(cpl(1).le.comp_tiny) cpl(1)=etaucr + end if + if(cpl(2).lt.etaucr) then ! . low coupled power for X-mode? => de-activate derived rays + call turnoffray(jk,ip+1,2*ib,iroff) + if(cpl(2).le.comp_tiny) cpl(2)=etaucr + end if + + taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass + cpls(jk,iO) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass + cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass + + if(jk.eq.1) then ! + polarization angles at plasma boundary for central ray + psipv(iO:iO+1) = psipol + chipv(iO:iO+1) = chipol + end if + end if + end if + end if + + iopmin = min(iopmin,iop(jk)) + if(ip.lt.ipass) then ! not last pass + yw0(:,jk) = yw(:,jk) ! * store current coordinates + ypw0(:,jk) = ypw(:,jk) ! in case pass ends on next step + end if + + ! compute ECRH&CD + if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. & + (tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check + tekev=temp(psinv) + call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & + sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd) + if(ierrhcd/=0) then + ierr = ior(ierr,ierrhcd) + call print_errhcd(ierrhcd,i,anprre,anprim,alpha) + end if + else + tekev=zero + alpha=zero + didp=zero + anprim=zero + anprre=anpr + nharm=0 + nhf=0 + iokhawa=0 + end if + if(nharm>0) iiv(jk)=i + + psjki(jk,i) = psinv + ! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s) + tau=tau0(jk)+0.5_wp_*(alphaabs0(jk)+alpha)*dersdst*dst + pow=p0ray(jk)*exp(-tau) !*exp(-tau1v(jk)) + ppabs(jk,i)=p0ray(jk)-pow + dids=didp*pow*alpha + ccci(jk,i)=ccci0(jk)+0.5_wp_*(dids0(jk)+dids)*dersdst*dst + tau0(jk)=tau + alphaabs0(jk)=alpha + dids0(jk)=dids + ccci0(jk)=ccci(jk,i) + + if(iwait(jk)) then + ppabs(jk,i:nstep) = ppabs(jk,i-1) + ccci(jk,i:nstep) = ccci(jk,i-1) + psjki(jk,i:nstep) = psjki(jk,i-1) + else + call print_output(i,jk,stv(jk),p0ray(jk)/etau1(jk),xv,psinv, & + btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, & + nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) !****** p0ray(jk)/etau1(jk) fattore normalizzazione per dids da controllare + end if + + end do +! ============ rays loop END ============ + + if(i==nstep) then ! step limit reached? + do jk=1,nray + if(iop(jk)<3) call turnoffray(jk,ip,ib,iroff) ! * ray hasn't exited+reentered the plasma by last step => stop ray + end do + end if + + ! print ray positions for j=nrayr in local reference system + if(mod(i,istpr0) == 0) then + if(nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,0) + end if + + ! check for any error code and stop if necessary + call check_err(ierr,istop) + ! test whether further trajectory integration is unnecessary + call vmaxmin(tau1+tau0+lgcpl1,nray,taumn,taumx) ! test on tau + coupling + if(taumn > taucr .or. all(iroff(:,index_rt))) istop = 1 ! (residual power~0) or (no ray active) => stop beam + + if(istop == 1) then ! stop propagation for current beam + istop_pass = istop_pass +1 ! * +1 non propagating beam + if(ip.lt.ipass) call turnoffray(0,ip,ib,iroff) ! * de-activate derived beams + exit + else if(all(iwait)) then ! all rays in current beam are waiting for next pass + exit + end if + end do +! ======== propagation loop END ======== + + ! print all ray positions in local reference system + if(nray > 1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,1) + +! =========== post-proc BEGIN =========== + ! compute total absorbed power and driven current for current beam + if(i>nstep) i=nstep + pabs_beam = sum(ppabs(:,i)) + icd_beam = sum(ccci(:,i)) + call vmaxmin(tau0,nray,taumn,taumx) ! taumn,taumx for print + + ! compute power and current density profiles for all rays + call spec(psjki,ppabs,ccci,iiv,pabs_beam,icd_beam,dpdv_beam,jphi_beam,jcd_beam, & + pins_beam,currins_beam) + + pabs_pass(iox) = pabs_pass(iox) + pabs_beam ! 0D results for current pass, sum on O/X mode beams + icd_pass(iox) = icd_pass(iox) + icd_beam + + if(ip.lt.ipass) then + cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1) / sum(p0ray * exp(-tau0)) ! average O-mode coupling + cpl_beam2 = one - cpl_beam1 + cpl_cbeam1 = cpls(1,iO)/cpl1(1) + cpl_cbeam2 = one - cpl_cbeam1 + else + cpl_beam1 = zero + cpl_beam2 = zero + end if + + ! print final results for pass on screen + write(*,*) + write(*,'("End of propagation for beam ",i5," (pass ",i3,", ",a1," mode)")'), & + index_rt,ip,mode(iox) + write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',stv(1) + write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx + write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs_beam + write(*,'(a,f9.4)') 'I_tot (kA) = ',icd_beam*1.0e3_wp_ + if(ip.lt.ipass) then + write(*,'(a,2(f9.4,1x))') 'Coupling (average, O/X):',cpl_beam1,cpl_beam2 ! average coupling for next O/X beams + write(*,'(a,2(f9.4,1x))') 'Coupling (c. beam, O/X):',cpl_cbeam1,cpl_cbeam2 ! central beam coupling for next O/X beams + end if + + call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, & + pins_beam,ip) ! *print power and current density profiles for current beam + + call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam,jphi_beam, & + rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip,rhotp,drhotp,rhotj, & + drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam + + call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav,rhotjava, & + drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx, & + ratjbmx,stv(1),psipv(index_rt),chipv(index_rt),index_rt,sum(p0ray), & + cpl_beam1,cpl_beam2) ! *print 0D results for current beam +! ============ post-proc END ============ + + end do +! ============ beam loop END ============ + +! ======= cumulative prints BEGIN ======= + pabs = pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output] + icd = icd + sum(icd_pass) + + ! print final results for pass on screen + write(*,*) + write(*,'("# End of pass ",i3)'),ip + write(*,'(a,f9.4,f9.4)') '# Pabs_tot (MW) [O,X mode] = ',pabs_pass(1),pabs_pass(2) + write(*,'(a,f9.4,f9.4)') '# I_tot (kA) [O,X mode] = ', & + icd_pass(1)*1.0e3_wp_,icd_pass(2)*1.0e3_wp_ +! ======== cumulative prints END ======== + + if(istop_pass == nbeam_pass) exit ! no active beams + end do +! ============ main loop END ============ -! compute total absorbed power and driven current - if (i>nstep) i=nstep - pabs = sum(ppabs(:,i)) - icd = sum(ccci(:,i)) -! ======= main loop END ====== - -! ======= post-proc BEGIN ====== -! print final results on screen - write(*,*) - write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',st - write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx - write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs - write(*,'(a,f9.4)') 'I_tot (kA) = ',icd*1.0e3_wp_ -! print all ray positions in local reference system - if(nray > 1) call print_projxyzt(st,yw,1) - -! compute power and current density profiles for all rays - call pec_init(ipec,rhout) - nnd=size(rhop_tab) - allocate(jphi(nnd),pins(nnd),currins(nnd)) - call spec(psjki,ppabs,ccci,iiv,pabs,icd,dpdv,jphi,jcd,pins,currins) -! 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 ====== +! print final results on screen + write(*,*) + write(*,'(a)') '## Final results:' + write(*,'(a,f9.4)') '## Pabs_tot (MW) = ',pabs + write(*,'(a,f9.4)') '## I_tot (kA) = ',icd*1.0e3_wp_ + +! ========== 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_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, & + alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) call dealloc_pec - deallocate(jphi,pins,currins) -! ======= free memory END ====== + call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0, & + stnext,stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & + pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) +! =========== free memory END =========== end subroutine gray_main @@ -320,7 +624,7 @@ contains subroutine ic_gb(xv0c,anv0c,ak0,wcsi,weta,rcicsi,rcieta,phiw,phir, & - ywrk0,ypwrk0,xc0,du10,gri,ggri) + ywrk0,ypwrk0,xc0,du10,gri,ggri,index_rt) ! 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 @@ -329,6 +633,7 @@ contains use beamdata, only : nray,nrayr,nrayth,rwmax implicit none ! arguments + integer, intent(in) :: index_rt real(wp_), dimension(3), intent(in) :: xv0c,anv0c real(wp_), intent(in) :: ak0 real(wp_), intent(in) :: wcsi,weta,rcicsi,rcieta,phiw,phir @@ -616,13 +921,13 @@ contains 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,0,1,ddr,ddi,zero,zero,(/zero,zero,zero/)) ! st=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0 + 0,0,0,index_rt,ddr,ddi,zero,zero,(/zero,zero,zero/)) ! st=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0 end do end subroutine ic_gb - subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr) + subroutine rkstep(sox,bres,xgcn,y,yp,dgr,ddgr,igrad) ! Runge-Kutta integrator use const_and_precisions, only : wp_ ! use gray_params, only : igrad @@ -633,6 +938,7 @@ contains real(wp_), dimension(6), intent(in) :: yp real(wp_), dimension(3), intent(in) :: dgr real(wp_), dimension(3,3), intent(in) :: ddgr + integer, intent(in) :: igrad real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4 real(wp_) :: gr2 @@ -645,18 +951,18 @@ contains fk1 = yp yy = y + fk1*hh - call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2) + call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk2,igrad) yy = y + fk2*hh - call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk3) + call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk3,igrad) yy = y + fk3*h - call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk4) + call rhs(sox,bres,xgcn,yy,gr2,dgr2,dgr,ddgr,fk4,igrad) y = y + h6*(fk1 + 2*fk2 + 2*fk3 + fk4) end subroutine rkstep - subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery) + subroutine rhs(sox,bres,xgcn,y,gr2,dgr2,dgr,ddgr,dery,igrad) ! Compute right-hand side terms of the ray equations (dery) ! used in R-K integrator use const_and_precisions, only : wp_ @@ -667,6 +973,7 @@ contains real(wp_), dimension(3), intent(in) :: dgr2,dgr real(wp_), dimension(3,3), intent(in) :: ddgr real(wp_), dimension(6), intent(out) :: dery + integer, intent(in) :: igrad ! local variables real(wp_) :: psinv,dens,btot,xg,yg,anpl,anpr,ajphi real(wp_) :: ddr,ddi,dersdst,derdnm @@ -679,13 +986,13 @@ contains 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) + dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad) end subroutine rhs subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & - bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr) + bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr,igrad) ! 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 @@ -701,6 +1008,7 @@ contains real(wp_), dimension(3), intent(out) :: bv integer, intent(out) :: ierr real(wp_), dimension(3), intent(out) :: derxg + integer, intent(in) :: igrad ! local variables real(wp_) :: gr2,ajphi real(wp_), dimension(3) :: dgr2,deryg @@ -711,7 +1019,7 @@ contains 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) + dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad) ierr=0 if( abs(anpl) > anplth1) then @@ -1082,9 +1390,9 @@ contains subroutine disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & - dery,anpl,anpr,ddr,ddi,dersdst,derdnm) + dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad) use const_and_precisions, only : wp_,zero,one,half,two - use gray_params, only : idst,igrad + use gray_params, only : idst implicit none ! arguments real(wp_), intent(in) :: xg,yg,gr2,sox @@ -1093,6 +1401,7 @@ contains real(wp_), dimension(3), intent(in) :: dgr2,dgr real(wp_), dimension(3,3), intent(in) :: ddgr,derbv real(wp_), dimension(6), intent(out) :: dery + integer, intent(in) :: igrad ! local variables integer :: iv real(wp_) :: yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl,an2,an2s @@ -1662,11 +1971,12 @@ bb: do - subroutine print_headers(strheader) + subroutine print_headers(strheader,index_rt) use units, only : uprj0,uwbm,udisp,ucenr,uoutr,upec,usumm implicit none ! arguments character(len=*), dimension(:), intent(in) :: strheader + integer, intent(in) :: index_rt ! local variables integer :: i,l @@ -1681,17 +1991,20 @@ bb: do write(upec,'(1x,a)') strheader(i) write(usumm,'(1x,a)') strheader(i) end do + + if(index_rt.gt.0) return + 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 '// & + write(ucenr,'(1x,a)') '#sst R z phi psin rhot ne Te Btot Bx By Bz Nperp Npl '// & 'Nx Ny Nz ki alpha tau Pt dIds nhmin nhmax iohkw index_rt ddr Xg Yg dXgdx dXgdy dXgdz' 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' + 'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt' // & + 'Jphimx dPdVmx drhotj drhotp P0 cplO cplX' end subroutine print_headers @@ -1872,13 +2185,13 @@ bb: do - subroutine print_projxyzt(st,ywrk,iproj) + subroutine print_projxyzt(stv,ywrk,iproj) use const_and_precisions, only : wp_, comp_huge, zero, one use beamdata, only : nray, nrayr, nrayth, rayi2jk use units, only : uprj0,uwbm implicit none ! arguments - real(wp_), intent(in) :: st + real(wp_), dimension(:), intent(in) :: stv real(wp_), dimension(:,:), intent(in) :: ywrk integer, intent(in) :: iproj ! local variables @@ -1921,14 +2234,14 @@ bb: do jkv=rayi2jk(jk) if(.not.(iproj==0 .and. jk==1)) & - write(uprj,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') st,jkv,xti,yti,zti,rti + write(uprj,'(1x,e16.8e3,2i5,4(1x,e16.8e3))') stv(jk),jkv,xti,yti,zti,rti if(iproj==1 .and. jkv(2)==nrayth) write(uprj,*) if(rti>=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 + write(uwbm,'(3(1x,e16.8e3))') stv(1),rtimn,rtimx end subroutine print_projxyzt @@ -2006,24 +2319,27 @@ bb: do 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 + write(upec,*) '' 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) + stmx,psipol,chipol,index_rt,p0,cpl1,cpl2) 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 + stmx,psipol,chipol,p0,cpl1,cpl2 integer, intent(in) :: index_rt - write(usumm,'(15(1x,e12.5),i5,4(1x,e12.5))') icd,pabs,jphip,dpdvp, & + 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 + stmx,psipol,chipol,index_rt,jphimx,dpdvmx,drhotj,drhotp,p0, & + cpl1,cpl2 + write(usumm,*) '' end subroutine print_finals end module graycore diff --git a/src/main-sum.f90 b/src/main-sum.f90 index 2bd9407..73c7ff6 100644 --- a/src/main-sum.f90 +++ b/src/main-sum.f90 @@ -46,7 +46,7 @@ program main_sum ! ======= read input data BEGIN ======= !------------ equilibrium ------------ if(eqp%iequil<2) then - call read_equil_an(eqp%filenm, rv, zv, fpol, qpsi) + 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 @@ -104,7 +104,7 @@ program main_sum ! 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_ + 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) diff --git a/src/main.f90 b/src/main.f90 index 7dca69c..8dc5c7e 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -38,7 +38,7 @@ program main_std ! ======= read input data BEGIN ======= !------------ equilibrium ------------ if(eqp%iequil<2) then - call read_equil_an(eqp%filenm, rv, zv, fpol, qpsi) + 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 @@ -96,7 +96,7 @@ program main_std ! 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_ + 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) diff --git a/src/multipass.f90 b/src/multipass.f90 new file mode 100755 index 0000000..e0393a9 --- /dev/null +++ b/src/multipass.f90 @@ -0,0 +1,280 @@ +module multipass + use const_and_precisions, only : wp_, zero, half, one, degree, czero + use beamdata, only : dst, nray + use gray_params, only : ipass + use polarization, only : pol_limit, stokes_ce, polellipse + use reflections, only : wall_refl + use equilibrium, only : bfield + + implicit none + + integer, save :: nbeam_max + integer, save :: nbeam_tot + +contains + +! ------------------------------ + subroutine plasma_in(i,xv,anv,bres,sox,cpl,psipol1,chipol1,iop,ext,eyt) ! ray enters plasma + implicit none + integer, intent(in) :: i ! ray index + real(wp_), dimension(3), intent(in) :: xv,anv + real(wp_), intent(in) :: bres,sox + real(wp_), dimension(2), intent(out) :: cpl ! coupling (O/X) + real(wp_), intent(out) :: psipol1,chipol1 + integer, dimension(:), intent(inout), pointer :: iop + complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt +! + real(wp_) :: rm,csphi,snphi,bphi,br,bz + real(wp_) :: qq1,uu1,vv1,qq,uu,vv,powcpl + real(wp_), dimension(3) :: bv,xmv + complex(wp_) :: ext1,eyt1 +! + iop(i)=iop(i)+1 + + xmv=xv*0.01_wp_ ! convert from cm to m + 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(1)=br*csphi-bphi*snphi + bv(2)=br*snphi+bphi*csphi + bv(3)=bz + + call pol_limit(anv,bv,bres,sox,ext1,eyt1) + call stokes_ce(ext1,eyt1,qq1,uu1,vv1) ! stokes parameter at plasma entrance + call stokes_ce(ext(i),eyt(i),qq,uu,vv) ! stokes parameter at plasma exit/wall reflection + powcpl = half*(one + vv*vv1+uu*uu1+qq*qq1) + + if(sox.eq.-one) then ! incoming mode = O + cpl=(/powcpl,one-powcpl/) + else ! incoming mode = X + cpl=(/one-powcpl,powcpl/) + end if + + if(i.eq.1) then + call polellipse(qq1,uu1,vv1,psipol1,chipol1) + psipol1=psipol1/degree ! convert from rad to degree + chipol1=chipol1/degree + else + psipol1 = zero + chipol1 = zero + end if + end subroutine plasma_in +! ------------------------------ + subroutine plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! ray exits plasma + implicit none + integer, intent(in) :: i ! ray index + real(wp_), dimension(3), intent(in) :: xv,anv + real(wp_), intent(in) :: bres,sox + integer, dimension(:), intent(inout), pointer :: iop + complex(wp_), dimension(:), intent(out), pointer :: ext,eyt +! + real(wp_) :: rm,csphi,snphi,bphi,br,bz + real(wp_), dimension(3) :: bv,xmv +! + iop(i)=iop(i)+1 + + xmv=xv*0.01_wp_ ! convert from cm to m + 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(1)=br*csphi-bphi*snphi + bv(2)=br*snphi+bphi*csphi + bv(3)=bz + call pol_limit(anv,bv,bres,sox,ext(i),eyt(i)) ! polarization at plasma exit + + end subroutine plasma_out +! ------------------------------ +! subroutine wall_in(i) ! ray enters vessel +! implicit none +! integer, intent(in) :: i ! ray index +! +! iow(i)=iow(i)+1 +! end subroutine wall_in +! ------------------------------ + subroutine wall_out(i,ins,xv,anv,bres,sox,psipol1,chipol1,iow,iop,ext,eyt) ! ray exits vessel + implicit none + integer, intent(in) :: i ! ray index + logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) + real(wp_), dimension(3), intent(inout) :: xv,anv + real(wp_), intent(in) :: bres,sox + real(wp_), intent(out) :: psipol1,chipol1 + integer, dimension(:), intent(inout), pointer :: iow,iop + complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt + ! + integer :: irfl + real(wp_) :: qq1,uu1,vv1 + real(wp_), dimension(3) :: xvrfl,anvrfl,walln + complex(wp_) :: ext1,eyt1 + ! + iow(i)=iow(i)+1 + if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) + call wall_refl(xv-dst*anv,anv,ext(i),eyt(i),xvrfl,anvrfl,ext1,eyt1,walln,irfl) + ext(i) = ext1 + eyt(i) = eyt1 + xv = xvrfl + anv = anvrfl + + if(i.eq.1) then ! polarization angles at wall reflection for central ray + call stokes_ce(ext1,eyt1,qq1,uu1,vv1) + call polellipse(qq1,uu1,vv1,psipol1,chipol1) + psipol1=psipol1/degree ! convert from rad to degree + chipol1=chipol1/degree + else + psipol1 = zero + chipol1 = zero + end if + + end subroutine wall_out +! ------------------------------ + subroutine initbeam(i,iroff,iboff,iwait,stv,jphi_beam,pins_beam,currins_beam, & + dpdv_beam,jcd_beam) ! initialization at beam propagation start + implicit none + integer, intent(in) :: i ! beam index + logical, dimension(:,:), intent(in), pointer :: iroff + logical, intent(out) :: iboff + logical, dimension(:), intent(out), pointer :: iwait + real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam, & + currins_beam,dpdv_beam,jcd_beam,stv + + iboff = .false. ! beam status (F = active, T = inactive) + iwait = iroff(:,i) ! copy ray status for current beam from global ray status + if(all(iwait)) then ! no rays active => stop beam + iboff = .true. + write(*,*) + write(*,'("Beam ",i5," inactive")'),i + else if(.not.all(.not.iwait)) then ! no all rays active + write(*,*) + write(*,'("WARNING: not all rays in beam ",i5," are active")'),i + end if + + stv = zero ! starting step + + jphi_beam = zero ! 1D beam profiles + pins_beam = zero + currins_beam = zero + dpdv_beam = zero + jcd_beam = zero + end subroutine initbeam +! ------------------------------ + subroutine initmultipass(i,iox,iroff,yynext,yypnext,yw0,ypw0,stnext,p0ray, & + taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) ! initialization before pass loop + implicit none + integer, intent(in) :: i, iox ! ipol, mode active on 1st pass + logical, dimension(:,:), intent(out), pointer :: iroff + real(wp_), dimension(:), intent(out), pointer :: p0ray,tau1,etau1,cpl1,lgcpl1, & + psipv,chipv + real(wp_), dimension(:,:), intent(out), pointer :: yw0,ypw0,stnext,taus,cpls + real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext +! + iroff = .false. ! global ray status (F = active, T = inactive) + if(i.eq.0) call turnoffray(0,1,3-iox,iroff) ! ipol = 0 => stop other mode (iox=1/2 -> stop ib=2/1) + yynext = zero ! starting beam coordinates (1) + yypnext = zero ! starting beam coordinates (2) + yw0 = zero ! temporary beam coordinates (1) + ypw0 = zero ! temporary beam coordinates (2) + stnext = zero ! starting beam step + p0ray = zero ! starting beam power + taus = zero ! starting beam tau + tau1 = zero ! total beam tau + etau1 = one + cpls = one + cpl1 = one + lgcpl1 = zero + psipv = zero ! psi polarization angle at vacuum-plasma boundary + chipv = zero ! chi polarization angle at vacuum-plasma boundary + end subroutine initmultipass +! ------------------------------ + subroutine turnoffray(jk,ip,ib,iroff) ! turn off ray propagation + implicit none + integer, intent(in) :: jk, ip, ib ! ray (0=all rays), pass, beam indexes + logical, dimension(:,:), intent(out), pointer :: iroff +! + integer :: ipx, i1, i2 +! + if(jk==0) then ! stop all rays + do ipx=ip,ipass + i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1 + i2 = 2**ipx-2+2**(ipx-ip)*ib + iroff(:,i1:i2) = .true. + end do + else ! only stop ray jk + do ipx=ip,ipass + i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1 + i2 = 2**ipx-2+2**(ipx-ip)*ib + iroff(jk,i1:i2) = .true. + end do + end if + + end subroutine turnoffray +! ------------------------------ + subroutine alloc_multipass(dim,iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & + stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & + pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) + implicit none + integer :: dim + logical, dimension(:), intent(out), pointer :: iwait + logical, dimension(:,:), intent(out), pointer :: iroff + integer, dimension(:), intent(out), pointer :: iop,iow + real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam,currins_beam, & + dpdv_beam,jcd_beam,stv,tau1,etau1,cpl1,lgcpl1,p0ray,psipv,chipv + real(wp_), dimension(:,:), intent(out), pointer :: taus,cpls,stnext,yw0,ypw0 + real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext + + call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext,stv, & + p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam,pins_beam,currins_beam, & + dpdv_beam,jcd_beam,psipv,chipv) + + nbeam_max = 2**ipass ! max n of beams active at a time + nbeam_tot = 2**(ipass+1)-2 ! total n of beams + + allocate(iwait(nray),iroff(nray,nbeam_tot),iop(nray),iow(nray), & + yynext(6,nray,nbeam_max-2),yypnext(6,nray,nbeam_max-2), & + yw0(6,nray),ypw0(6,nray),stnext(nray,nbeam_tot),stv(nray), & + p0ray(nray),taus(nray,nbeam_tot),tau1(nray),etau1(nray), & + cpls(nray,nbeam_tot),cpl1(nray),lgcpl1(nray),jphi_beam(dim), & + pins_beam(dim),currins_beam(dim),dpdv_beam(dim),jcd_beam(dim), & + psipv(nbeam_tot),chipv(nbeam_tot)) + + end subroutine alloc_multipass +! ------------------------------ + subroutine dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & + stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & + pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) + implicit none + logical, dimension(:), intent(out), pointer :: iwait + logical, dimension(:,:), intent(out), pointer :: iroff + integer, dimension(:), intent(out), pointer :: iop,iow + real(wp_), dimension(:), intent(out), pointer :: stv,p0ray,tau1,etau1,cpl1,lgcpl1, & + jphi_beam,pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv + real(wp_), dimension(:,:), intent(out), pointer :: yw0,ypw0,stnext,taus,cpls + real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext + + if (associated(iwait)) deallocate(iwait) + if (associated(iroff)) deallocate(iroff) + if (associated(iop)) deallocate(iop) + if (associated(iow)) deallocate(iow) + if (associated(yynext)) deallocate(yynext) + if (associated(yypnext)) deallocate(yypnext) + if (associated(yw0)) deallocate(yw0) + if (associated(ypw0)) deallocate(ypw0) + if (associated(stnext)) deallocate(stnext) + if (associated(stv)) deallocate(stv) + if (associated(p0ray)) deallocate(p0ray) + if (associated(taus)) deallocate(taus) + if (associated(tau1)) deallocate(tau1) + if (associated(etau1)) deallocate(etau1) + if (associated(cpls)) deallocate(cpls) + if (associated(cpl1)) deallocate(cpl1) + if (associated(lgcpl1)) deallocate(lgcpl1) + if (associated(jphi_beam)) deallocate(jphi_beam) + if (associated(pins_beam)) deallocate(pins_beam) + if (associated(currins_beam)) deallocate(currins_beam) + if (associated(dpdv_beam)) deallocate(dpdv_beam) + if (associated(jcd_beam)) deallocate(jcd_beam) + if (associated(psipv)) deallocate(psipv) + if (associated(chipv)) deallocate(chipv) + end subroutine dealloc_multipass + +end module multipass diff --git a/src/pec.f90 b/src/pec.f90 index 5151da9..f0b0e62 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -174,7 +174,7 @@ contains if(ise0 == 0) then if(xxi(i) < one) then ise0 = i - isev(is) = i - 1 + isev(is) = max(i - 1, 1) is = is + 1 end if else @@ -219,9 +219,14 @@ contains if (idecr == -1) indi = ind - 1 rt1 = xtab1(indi) call locatex(xxi,iise,iise0,iise,rt1,itb1) - if(itb1 >= iise0 .and. itb1 < iise) then - call intlin(xxi(itb1), ypt(itb1),xxi(itb1+1), ypt(itb1+1),rt1,ppa2) - call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1),rt1,cci2) + if(itb1 >= iise0 .and. itb1 <= iise) then + if(itb1 == iise) then + ppa2=ypt(itb1) + cci2=yamp(itb1) + else + call intlin(xxi(itb1), ypt(itb1),xxi(itb1+1), ypt(itb1+1),rt1,ppa2) + call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1),rt1,cci2) + end if dppa = ppa2 - ppa1 didst = cci2 - cci1 wdpdv(ind) = wdpdv(ind) + dppa @@ -238,7 +243,7 @@ contains rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, & rhotp, drhotp, rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx) ! radial average values over power and current density profile - use const_and_precisions, only : pi + use const_and_precisions, only : pi, comp_eps use gray_params, only : nnd use equilibrium, only : frhopol use magsurf_data, only : fluxval @@ -277,8 +282,8 @@ contains end if ! factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile - drhotpav = sqrt(8._wp_*(rhot2pav -rhotpav**2)) - drhotjava = sqrt(8._wp_*(rhot2java-rhotjava**2)) + drhotpav = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps))) + drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps))) rhoppav = frhopol(rhotpav) rhopjava = frhopol(rhotjava) diff --git a/src/polarization.f90 b/src/polarization.f90 index 75c0390..4fb65d9 100644 --- a/src/polarization.f90 +++ b/src/polarization.f90 @@ -86,6 +86,9 @@ contains if (den>zero) then ext = (ff*csgam - ui*anpl*sngam)/sqrt(den) eyt = (-ff*sngam - ui*anpl*csgam)/sqrt(den) + den = sqrt(abs(ext)**2+abs(eyt)**2) + ext = ext/den + eyt = eyt/den else ! only for XM (sox=+1) when N//=0 ext = -ui*sngam eyt = -ui*csgam