From f5ab40f54f93c597c5e785e595b5a8873dc00dc1 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Tue, 10 May 2022 12:36:37 +0200 Subject: [PATCH] src/gray_core.f90: cleanup - add some comments - annotate loops - indent comments - remove trailing whitespace - reduce usage of opaque global variables - use Fortran 90 logical operators - use Fortran 2003 array syntax --- src/beamdata.f90 | 3 +- src/gray_core.f90 | 907 ++++++++++++++++++++++++-------------------- src/gray_params.f90 | 2 +- src/main.f90 | 11 +- 4 files changed, 504 insertions(+), 419 deletions(-) diff --git a/src/beamdata.f90 b/src/beamdata.f90 index 6b18c27..ee8c768 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -12,7 +12,7 @@ contains use gray_params, only : raytracing_parameters use const_and_precisions, only : half,two implicit none - type(raytracing_parameters), intent(in) :: rtrparam + type(raytracing_parameters), intent(inout) :: rtrparam real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, & gri,psjki,ppabs,ccci real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri @@ -41,6 +41,7 @@ contains end if nray = (nrayr-1)*nrayth + 1 jkray1 = (jray1-2)*nrayth + 2 + rtrparam%nray = nray if(nrayr>1) then twodr2 = two*(rwmax/(nrayr-1))**2 diff --git a/src/gray_core.f90 b/src/gray_core.f90 index aff88c9..b49337b 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -12,13 +12,12 @@ contains use coreprofiles, only : temp, fzeff use dispersion, only : expinit use gray_params, only : gray_parameters, gray_data, gray_results, & - print_parameters, & - iwarm, ipec, istpr0, igrad, headw, headl, ipass + print_parameters use beams, only : xgygcoeff, launchangles2n use beamdata, only : pweight, rayi2jk 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 beamdata, only : init_btr, dealloc_beam use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & rhop_tab, rhot_tab use utils, only : vmaxmin @@ -31,9 +30,9 @@ contains implicit none ! subroutine arguments - type(gray_parameters), intent(in) :: params - type(gray_data), intent(in) :: data - type(gray_results), intent(out) :: results + type(gray_parameters), intent(inout) :: params + type(gray_data), intent(in) :: data + type(gray_results), intent(out) :: results ! Predefined grid for the output profiles (optional) real(wp_), dimension(:), intent(in), optional :: rhout @@ -43,7 +42,7 @@ contains ! local variables real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) - character, dimension(2), parameter :: mode=(/'O','X'/) + 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 @@ -104,13 +103,13 @@ contains tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv) ! Initialise the dispersion module - if(iwarm > 1) call expinit + if(params%ecrh_cd%iwarm > 1) call expinit ! Initialise the magsurf_data module call flux_average ! requires frhotor for dadrhot,dvdrhot - + ! Initialise the output profiles - call pec_init(ipec, rhout) + call pec_init(params%output%ipec, rhout) nnd = size(rhop_tab) ! number of radial profile points call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & @@ -155,18 +154,18 @@ contains call initmultipass(params%raytracing%ipol, params%antenna%iox, & iroff,yynext,yypnext,yw0,ypw0, & stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) - - if(params%raytracing%ipol .eq. 0) then - if(params%antenna%iox .eq. 2) then ! only X mode on 1st pass - cpl0 = (/zero,one/) + + if(params%raytracing%ipol == 0) then + if(params%antenna%iox == 2) then ! only X mode on 1st pass + cpl0 = [zero, one] else ! only O mode on 1st pass - cpl0 = (/one,zero/) + cpl0 = [one, zero] end if end if sox=one ! mode inverted for each beam iox=2 ! start with O: sox=-1, iox=1 - + psipol = params%antenna%psi chipol = params%antenna%chi call pweight(params%antenna%power, p0jk) @@ -175,7 +174,7 @@ contains index_rt=0 ! global beam index: 1,O 2,X 1st pass ! | | | | call log_debug('pass loop start', mod='gray_core', proc='gray_main') ! 3,O 4,X 5,O 6,X 2nd pass - do ip=1,ipass + main_loop: do ip=1,params%raytracing%ipass write (msg, '("pass: ",g0)') ip call log_info(msg, mod='gray_core', proc='gray_main') @@ -183,21 +182,22 @@ contains 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 + + if(ip > 1) then du1 = zero gri = zero ggri = zero - if(ip.eq.ipass) cpl = (/zero,zero/) ! no successive passes + if(ip == params%raytracing%ipass) cpl = [zero, zero] ! no successive passes end if - -! =========== beam loop BEGIN =========== + + ! =========== beam loop BEGIN =========== call log_debug('beam loop start', mod='gray_core', proc='gray_main') - do ib=1,nbeam_pass - + + beam_loop: 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) @@ -212,29 +212,27 @@ contains call log_info(" beam is off", mod='gray_core', proc='gray_main') cycle end if - + call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) - - if(ip.eq.1) then ! 1st pass - igrad_b = igrad ! * input value, igrad_b=0 from 2nd pass - + + if(ip == 1) then ! 1st pass + igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass + tau1 = zero ! * tau from previous passes etau1 = one cpl1 = one ! * coupling from previous passes lgcpl1 = zero p0ray = p0jk ! * initial beam power - call ic_gb(params%antenna%pos, anv0, ak0, & - params%antenna%w(1),params%antenna%w(2), & - params%antenna%ri(1),params%antenna%ri(2), & - params%antenna%phi(1),params%antenna%phi(2), & - 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 + call ic_gb(params, anv0, ak0, yw, ypw, & ! * initial conditions + xc, du1, gri, ggri, index_rt) + call set_pol(yw, bres, sox, params%raytracing%ipol, & ! * initial polarization + psipol, chipol, ext, eyt) + + do jk=1,params%raytracing%nray 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(data%equilibrium%rlim, data%equilibrium%zlim, & nlim, rrm, zzm)) then ! * start propagation in/outside vessel? iow(jk) = 1 ! + inside @@ -242,7 +240,7 @@ contains iow(jk) = 0 ! + outside end if end do - + else ! 2nd+ passes ipar = (index_rt+1)/2-1 ! * parent beam index yw = yynext(:,:,ipar) ! * starting coordinates from @@ -257,63 +255,64 @@ contains 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 ======= + + if(params%raytracing%nray>1 .and. all(.not.iwait)) & ! iproj=0 ==> nfilp=8 + call print_projxyzt(stv,yw,0) + + ! ======= propagation loop BEGIN ======= call log_debug(' propagation loop start', mod='gray_core', proc='gray_main') - do i=1,nstep - ! advance one step with "frozen" grad(S_I) - do jk=1,nray + propagation_loop: do i=1,params%raytracing%nstep + ! advance one step with "frozen" grad(S_I) + do jk=1,params%raytracing%nray if(iwait(jk)) cycle ! jk ray is waiting for next pass - stv(jk) = stv(jk) + dst ! current ray step + stv(jk) = stv(jk) + params%raytracing%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 + ! update position and grad if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) - + error = 0 istop = 0 ! stop flag for current beam iopmin = 10 - -! =========== rays loop BEGIN =========== - do jk=1,nray + + ! =========== rays loop BEGIN =========== + rays_loop: do jk=1,params%raytracing%nray if(iwait(jk)) cycle ! jk ray is waiting for next pass - - ! compute derivatives with updated gradient and local plasma values + + ! 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 + ! update global error code and print message if(ierrn/=0) then error = ior(error,ierrn) call print_errn(ierrn,i,anpl) end if - - ! check entrance/exit plasma/wall + + ! 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=zero .and. psinv continue current pass - - if(params%raytracing%ipol .eq. 0) then ! + IF single mode propagation + if(iop(jk) == 1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass + + if(params%raytracing%ipol == 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 + else if(cpl(iox) < 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 @@ -321,58 +320,58 @@ contains p0ray(jk) = p0ray(jk)*cpl(iox) end if cpls(jk,index_rt) = cpl(iox) - - if(jk.eq.1) then + + if(jk == 1) then write (msg,'(" 1st pass - central ray (",a1,"-mode) c=",g0.4)') & mode(iox), cpl(iox) call log_info(msg, mod='gray_core', proc='gray_main') 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 + + else if(iop(jk) > 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 + + if(ip < params%raytracing%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 + stnext(jk,index_rt) = stv(jk) - params%raytracing%dst ! . starting step for next pass = last step + + if(cpl(1) < 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 + if(cpl(2) < 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 + + if(jk == 1) then ! . polarization angles at plasma boundary for central ray psipv(iO:iO+1) = psipol chipv(iO:iO+1) = chipol end if else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass - cpl = (/zero,zero/) + 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 ! * out->in - + 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) + 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 @@ -380,55 +379,58 @@ contains error = ior(error,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 + + if(jk == 1 .and. ip == 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+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 + + if(ip < params%raytracing%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 + + if(cpl(1) < 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 + if(cpl(2) < 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 + if(jk == 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 + if(ip < params%raytracing%ipass) then ! not last pass yw0(:,jk) = yw(:,jk) ! * store current coordinates in case ypw0(:,jk) = ypw(:,jk) ! current pass ends on next step end if - - ! compute ECRH&CD if (inside plasma & power available>0 & ray still active) - if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. & + + ! Compute ECRH&CD if (inside plasma & power available>0 & ray still active) + ! Note: power check is τ + τ₀ + log(coupling O/X) < τ_critic + + if(ierrn==0 .and. params%ecrh_cd%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) + call alpha_effj(params%ecrh_cd, psinv, xg, yg, dens, tekev, & + ak0, bres, derdnm, anpl, anpr, sox, anprre, & + anprim, alpha, didp, nharm, nhf, iokhawa, ierrhcd) if(ierrhcd/=0) then error = ior(error,ierrhcd) call print_errhcd(ierrhcd,i,anprre,anprim,alpha) @@ -444,83 +446,90 @@ contains 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) - + + ! Computation of the ray τ, dP/ds, P(s), dI/ds, I(s) + + ! optical depth: τ = ∫α(s)ds using the trapezoid rule + tau = tau0(jk) + 0.5_wp_*(alphaabs0(jk) + alpha) * dersdst * params%raytracing%dst + + pow = p0ray(jk) * exp(-tau) ! residual power: P = P₀exp(-τ) + ppabs(jk,i) = p0ray(jk) - pow ! absorbed power: P_abs = P₀ - P + dids = didp * pow * alpha ! current drive: dI/ds = dI/dP⋅P⋅α + + ! current: I = ∫dI/ds⋅ds using the trapezoid rule + ccci(jk,i) = ccci0(jk) + 0.5_wp_*(dids0(jk) + dids) * dersdst * params%raytracing%dst + + tau0(jk) = tau + alphaabs0(jk) = alpha + dids0(jk) = dids + ccci0(jk) = ccci(jk,i) + if(iwait(jk)) then ! copy values from last pass for inactive ray - ppabs(jk,i:nstep) = ppabs(jk,i-1) - ccci(jk,i:nstep) = ccci(jk,i-1) - psjki(jk,i:nstep) = psjki(jk,i-1) + ppabs(jk,i:params%raytracing%nstep) = ppabs(jk,i-1) + ccci(jk,i:params%raytracing%nstep) = ccci(jk,i-1) + psjki(jk,i:params%raytracing%nstep) = psjki(jk,i-1) else call print_output(i,jk,stv(jk),p0ray(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/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) end if - - end do -! ============ rays loop END ============ - - if(i==nstep) then ! step limit reached? - do jk=1,nray + + end do rays_loop + ! ============ rays loop END ============ + + if(i==params%raytracing%nstep) then ! step limit reached? + do jk=1,params%raytracing%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) + + ! print ray positions for j=nrayr in local reference system + if(mod(i,params%output%istpr) == 0) then + if(params%raytracing%nray > 1 .and. all(.not.iwait)) & + call print_projxyzt(stv,yw,0) end if - - ! check for any error code and stop if necessary + + ! check for any error code and stop if necessary call check_err(error,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 - + ! test whether further trajectory integration is unnecessary + call vmaxmin(tau1+tau0+lgcpl1,params%raytracing%nray,taumn,taumx) ! test on tau + coupling + 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 + if(ip < params%raytracing%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 => do not increase istop_pass exit end if - end do + end do propagation_loop call log_debug(' propagation loop end', mod='gray_core', proc='gray_main') -! ======== propagation loop END ======== + ! ======== 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 + ! print all ray positions in local reference system + if(params%raytracing%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>params%raytracing%nstep) i=params%raytracing%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 vmaxmin(tau0,params%raytracing%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 .and. iopmin.gt.2) then ! not last pass AND at least one ray re-entered plasma - cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1, MASK=iop.gt.2) / & - sum(p0ray * exp(-tau0), MASK=iop.gt.2) ! * average O-mode coupling for next beam (on active rays) + + if(ip < params%raytracing%ipass .and. iopmin > 2) then ! not last pass AND at least one ray re-entered plasma + cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1, MASK=iop > 2) / & + sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays) cpl_beam2 = one - cpl_beam1 ! * average X-mode coupling for next beam - - if(iop(1).gt.2) then ! * central ray O/X-mode coupling for next beam + + if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam cpl_cbeam1 = cpls(1,iO)/cpl1(1) cpl_cbeam2 = one - cpl_cbeam1 end if @@ -543,7 +552,7 @@ contains write(msg, '(3x,a,g0.3," MW")') 'current drive: I=', icd_beam*1.0e3_wp_ call log_info(msg, mod='gray_core', proc='gray_main') - if(ip < ipass) then + if(ip < params%raytracing%ipass) then write (msg,'(3x,a,(g0.4,", ",g0.4))') & ! average coupling for next O/X beams (=0 if no ray re-entered plasma) 'next couplings [O,X mode]: c=', cpl_beam1, cpl_beam2 call log_info(msg, mod='gray_core', proc='gray_main') @@ -555,22 +564,22 @@ contains 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 + ! ============ post-proc END ============ + + end do beam_loop call log_debug('beam loop end', mod='gray_core', proc='gray_main') -! ============ beam loop END ============ - -! ======= cumulative prints BEGIN ======= + ! ============ beam loop END ============ + + ! ======= cumulative prints BEGIN ======= results%pabs = results%pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output] results%icd = results%icd + sum(icd_pass) @@ -584,14 +593,14 @@ contains write(msg, '(" current drive [O,X mode] I=",g0.4,", ",g0.4," kA")') & icd_pass(1)*1.0e3_wp_, icd_pass(2)*1.0e3_wp_ call log_info(msg, mod='gray_core', proc='gray_main') -! ======== cumulative prints END ======== - - if(istop_pass == nbeam_pass) exit ! no active beams - end do - call log_debug('pass loop end', mod='gray_core', proc='gray_main') -! ============ main loop END ============ + ! ======== cumulative prints END ======== -! ========== free memory BEGIN ========== + if(istop_pass == nbeam_pass) exit ! no active beams + end do main_loop + call log_debug('pass loop end', mod='gray_core', proc='gray_main') + ! ============ main loop END ============ + + ! ========== free memory BEGIN ========== call dealloc_surfvec call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, & alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) @@ -599,7 +608,7 @@ contains 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 =========== + ! =========== free memory END =========== end subroutine gray_main @@ -635,26 +644,30 @@ contains - subroutine ic_gb(xv0c,anv0c,ak0,wcsi,weta,rcicsi,rcieta,phiw,phir, & - ywrk0,ypwrk0,xc0,du10,gri,ggri,index_rt) -! beam tracing initial conditions igrad=1 -! !!!!!! check ray tracing initial conditions igrad=0 !!!!!! + subroutine ic_gb(params, anv0c, ak0, 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 : zero,one,pi,half,two,degree,ui=>im - use math, only : catand - use gray_params, only : idst - use beamdata, only : nray,nrayr,nrayth,rwmax + use math, only : catand + use gray_params, only : gray_parameters + 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 - 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 + + ! subroutine arguments + type(gray_parameters), intent(in) :: params + real(wp_), dimension(3), intent(in) :: anv0c + real(wp_), intent(in) :: ak0 + real(wp_), dimension(6, nray), intent(out) :: ywrk0, ypwrk0 + real(wp_), dimension(3, nrayth, nrayr), intent(out) :: xc0, du10 + real(wp_), dimension(3, nray), intent(out) :: gri + real(wp_), dimension(3, 3, nray), intent(out) :: ggri + integer, intent(in) :: index_rt + + ! local variables + real(wp_), dimension(3) :: xv0c + real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir 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 @@ -670,6 +683,15 @@ contains complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1,dqi2 complex(wp_) :: dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy + ! initial beam parameters + xv0c = params%antenna%pos + wcsi = params%antenna%w(1) + weta = params%antenna%w(2) + rcicsi = params%antenna%ri(1) + rcieta = params%antenna%ri(2) + phiw = params%antenna%phi(1) + phir = params%antenna%phi(2) + csth=anv0c(3) snth=sqrt(one-csth**2) if(snth > zero) then @@ -680,20 +702,20 @@ contains 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 + ! 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) + !csphir = cos(phirrad) + !snphir = sin(phirrad) wwcsi = two/(ak0*wcsi**2) wweta = two/(ak0*weta**2) @@ -723,11 +745,11 @@ contains 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) + + !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 @@ -768,7 +790,7 @@ contains 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 + end do da=2*pi/dble(nrayth) do k=1,nrayth @@ -777,7 +799,7 @@ contains csa(k) = cos(alfak) end do -! central ray + ! central ray jk=1 gri(:,1) = zero ggri(:,:,1) = zero @@ -786,7 +808,7 @@ contains 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 @@ -803,7 +825,7 @@ contains ddr = zero ddi = zero -! loop on rays jk>1 + ! loop on rays jk>1 j=2 k=0 do jk=2,nray @@ -813,10 +835,10 @@ contains k=1 end if -! csiw=u*dcsiw -! etaw=u*detaw -! csir=x0t*csphir+y0t*snphir -! etar=-x0t*snphir+y0t*csphir + !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 @@ -848,7 +870,7 @@ contains 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 @@ -911,16 +933,17 @@ contains 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 + ! integration variable + select case(params%raytracing%idst) + case(0) + ! optical path: s + denom = an0 + case(1) + ! "time": c⋅t + denom = one + case(2) + ! real eikonal: S_R + denom = an20 end select ypwrk0(1,jk) = anx/denom ypwrk0(2,jk) = any/denom @@ -931,9 +954,10 @@ contains ddr = anx**2 + any**2 + anz**2 - an20 ddi = 2*(anxt*gxt + anyt*gyt + anzt*gzt) - call print_output(0,jk,zero,one,xc0(:,k,j),-one,zero,(/zero,zero,zero/), & - ak0,zero,zero,(/zero,zero,zero/),zero,zero,zero,zero,zero,zero, & - 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 + ! st=0, index_rt=1, B=0, N=0, psin=-1, Xg=0, Yg=0, gradXg=0 + 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,index_rt,ddr,ddi,zero,zero,[zero,zero,zero]) end do end subroutine ic_gb @@ -954,7 +978,7 @@ contains real(wp_), dimension(6) :: yy,fk1,fk2,fk3,fk4 real(wp_) :: gr2 real(wp_), dimension(3) :: dgr2 - + gr2 = dot_product(dgr, dgr) dgr2 = 2 * matmul(ddgr, dgr) fk1 = yp @@ -1033,9 +1057,9 @@ contains if( abs(anpl) > anplth1) then if(abs(anpl) > anplth2) then error=ibset(error,pnpl+1) - else + else error=ibset(error,pnpl) - end if + end if end if end subroutine ywppla_upd @@ -1091,7 +1115,7 @@ contains 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 @@ -1182,7 +1206,7 @@ contains ggri(2,3,jk)=gyz ggri(3,3,jk)=gzz end do - + end subroutine gradi_upd @@ -1241,20 +1265,23 @@ contains - subroutine plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv, & - xg,yg,derxg,deryg,ajphi) - use const_and_precisions, only : zero,ccj=>mu0inv - use gray_params, only : iequil - use equilibrium, only : psia,equinum_fpol,equinum_psi,equian,sgnbphi - use coreprofiles, only : density + subroutine plas_deriv(xv, bres, xgcn, psinv, dens, btot, bv, derbv, & + xg, yg, derxg, deryg, ajphi) + use const_and_precisions, only : zero, ccj=>mu0inv + use gray_params, only : iequil + use equilibrium, only : psia, equinum_fpol, equinum_psi, equian, sgnbphi + use coreprofiles, only : density + implicit none -! arguments + + ! subroutine 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 + + ! local variables integer :: jv real(wp_) :: xx,yy,zz real(wp_) :: csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm @@ -1286,7 +1313,7 @@ contains yy = xv(2) zz = xv(3) -! cylindrical coordinates + ! cylindrical coordinates rr2 = xx**2 + yy**2 rr = sqrt(rr2) csphi = xx/rr @@ -1295,7 +1322,7 @@ contains bv(1) = -snphi*sgnbphi bv(2) = csphi*sgnbphi -! convert from cm to meters + ! convert from cm to meters zzm = 1.0e-2_wp_*zz rrm = 1.0e-2_wp_*rr @@ -1307,7 +1334,7 @@ contains call equinum_fpol(psinv,fpolv,dfpv) end if -! compute yg and derivative + ! compute yg and derivative if(psinv < zero) then bphi = fpolv/rrm btot = abs(bphi) @@ -1315,27 +1342,25 @@ contains return end if -! compute xg and derivative + ! 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 + + ! 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 + ! bvc(i) = B_i in cylindrical coordinates + bvc = [brr, bphi, bzz] -! bv(i) = B_i in cartesian coordinates + ! 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(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 @@ -1343,7 +1368,7 @@ contains dbvcdc(2,3) = dfpv*dpsidz/rrm dbvcdc(3,3) = ddpsidrz/rrm -! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) + ! 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) @@ -1359,22 +1384,22 @@ contains dphidx = -snphi/rrm dphidy = csphi/rrm -! dbv(iv,jv) = d Bcart(iv) / dxvcart(jv) + ! 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 + ! B magnitude and derivatives btot = norm2(bv) dbtot = matmul(bv, dbv)/btot yg = btot/Bres -! convert spatial derivatives from dummy/m -> dummy/cm -! to be used in rhs + ! 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 + ! 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 @@ -1385,21 +1410,24 @@ contains 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 + ! 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)) + !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,igrad) - use const_and_precisions, only : zero,one,half,two - use gray_params, only : idst + subroutine disp_deriv(anv, sox, xg, yg, derxg, deryg, bv, derbv, & + gr2, dgr2, dgr, ddgr, dery, anpl, anpr, & + ddr, ddi, dersdst, derdnm, igrad) + use const_and_precisions, only : zero, one, half, two + use gray_params, only : idst + implicit none -! arguments + + ! subroutine 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 @@ -1407,12 +1435,12 @@ contains 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 - 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 + + ! local variables + 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 + real(wp_), dimension(3) :: derdxv, danpldxv, derdnv, dbgr an2 = dot_product(anv, anv) anpl = dot_product(anv, bv) @@ -1435,17 +1463,28 @@ contains duh = one - xg - yg2 if(xg > zero) then + ! Derivatives of the cold plasma dispersion relation + ! + ! Λ = N²⊥ - N²s(X,Y,N∥) = 0 + ! + ! where N²s = 1 - X - XY²⋅(1 + N∥ ∓ √Δ)/[2(1 - X - Y²)] + ! Δ = (1 - N∥)² + 4N∥⋅(1 - X)/Y² + ! del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2) an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh + ! ∂(N²s)/∂X dan2sdxg = - half*yg2*(one - yg2)*(one + anpl2 + sox*del)/duh**2 & + sox*xg*anpl2/(del*duh) - one + ! ∂(N²s)/∂Y dan2sdyg = - xg*yg*(one - xg)*(one + anpl2 + sox*del)/duh**2 & + two*sox*xg*(one - xg)*anpl2/(yg*del*duh) + ! ∂(N²s)/∂N∥ dan2sdnpl = - xg*yg2*anpl/duh & - sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh) if(igrad > 0) then + ! Derivatives of the eikonal (beamtracing only) 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 @@ -1464,13 +1503,9 @@ contains end if end if - bdotgr = dot_product(bv, dgr) - 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 + bdotgr = dot_product(bv, dgr) + dbgr = matmul(dgr, derbv) + matmul(bv, ddgr) + danpldxv = matmul(anv, derbv) derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl + & igrad*dgr2) & @@ -1485,147 +1520,197 @@ contains + half*yg*dfdiadyg & + half*anpl*dfdiadnpl) - ! integration variable - if (idst == 0) then - ! optical path: s - denom = derdnm - else if (idst == 1) then - ! "time": c*t - denom = -derdom - else - ! real eikonal: S_R - denom = dot_product(anv, derdnv) - end if + ! integration variable σ + select case(idst) + case(0) + ! optical path: s + denom = derdnm + case(1) + ! "time": c⋅t + denom = -derdom + case(2) + ! real eikonal: S_R + denom = dot_product(anv, derdnv) + end select -! coefficient for integration in s -! ds/dst, where st is the integration variable + ! 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 + ! r.h.s. vector + dery(1:3) = derdnv(:)/denom ! ∂Λ/∂N̅ + dery(4:6) = -derdxv(:)/denom ! -∂Λ/∂x̅ -! 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) + ! dispersion relation + ddr = an2 - an2s - igrad*(gr2 - half*bdotgr**2*fdia) ! real part + ddi = dot_product(derdnv, dgr) ! imaginary part 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,error) - use const_and_precisions, only : 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 + subroutine alpha_effj(params, psinv, xg, yg, dens, tekev, ak0, bres, & + derdnm, anpl, anpr, sox, anprre, anprim, alpha, & + didp, nhmin, nhmax, iokhawa, error) + ! Computes the absorption coefficient and effective current + + use const_and_precisions, only : zero, pi, mc2=>mc2_ + use gray_params, only : ecrh_cd_parameters + 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) :: error -! local variables - real(wp_) :: rbavi,rrii,rhop - integer :: lrm,ithn,ierrcd - real(wp_) :: amu,ratiovgr,rbn,rbx - real(wp_) :: zeff,cst2,bmxi,bmni,fci + + ! subroutine arguments + + ! Inputs + + ! ECRH & CD parameters + type(ecrh_cd_parameters) :: params + ! poloidal flux ψ + real(wp_), intent(in) :: psinv + ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω + real(wp_), intent(in) :: xg, yg + ! density [10¹⁹ m⁻³], temperature [keV] + real(wp_), intent(in) :: dens, tekev + ! vacuum wavenumber k₀=ω/c, resonant B field + real(wp_), intent(in) :: ak0, bres + ! group velocity: |∂Λ/∂N̅| where Λ=c²/ω²⋅D_R + real(wp_), intent(in) :: derdnm + ! refractive index: N∥, N⊥ (cold dispersion) + real(wp_), intent(in) :: anpl, anpr + ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + real(wp_), intent(in) :: sox + + ! Outputs + + ! real/imaginary parts of N⊥ (warm dispersion) + real(wp_), intent(out) :: anprre, anprim + ! absorption coefficient, current density + real(wp_), intent(out) :: alpha, didp + ! lowest/highest resonant harmonic numbers + integer, intent(out) :: nhmin, nhmax + ! ECCD/overall error codes + integer, intent(out) :: iokhawa, error + + ! local variables + real(wp_) :: rbavi, rrii, rhop + integer :: lrm, ithn, ierrcd + real(wp_) :: amu, 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 + real(wp_) :: effjcd, effjcdav, btot + real(wp_) :: akim + complex(wp_) :: ex, ey, ez - alpha=zero - anprim=zero - anprre=zero - didp=zero - nhmin=0 - nhmax=0 - iokhawa=0 - error=0 + alpha = zero + anprim = zero + anprre = zero + didp = zero + nhmin = 0 + nhmax = 0 + iokhawa = 0 + error = 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,error,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 + ! Absorption computation - 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 - error=error+ierrcd - if(associated(eccdpar)) deallocate(eccdpar) - effjcdav=rbavi*effjcd - didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii) - end if - end if + ! Skip when the temperature is zero (no plasma) + if (tekev <= zero) return + + ! Skip when the lowest resonant harmonic number is zero + amu = mc2/tekev ! μ=(mc²/kT) + call harmnumber(yg, amu, anpl, nhmin, nhmax, params%iwarm) + if (nhmin <= 0) return + + ! Solve the dispersion relation + lrm = max(params%ilarm, nhmax) + call warmdisp(xg, yg, amu, anpl, anpr, sox, lrm, error, & + anprre, anprim, params%iwarm, params%imx, & + ex, ey, ez) + + ! The absoption coefficient is defined as + ! + ! α = 2 Im(k̅)⋅s̅ + ! + ! where s̅ = v̅_g/|v_g|, the direction of the energy flow. + ! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω we have: + ! + ! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅| + ! = -2N̅ / |∂Λ/∂N̅| (using the cold dispersion) + ! + ! Assuming Im(k∥)=0: + ! + ! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅| + ! + akim = ak0 * anprim ! imaginary part of k⊥ = k₀N⊥ + alpha = 4 * akim * anpr/derdnm + + if (alpha < zero) then + error = ibset(error, palph) + return end if + + ! Current drive computation + if (params%ieccd <= 0) return + + 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(params%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 + error = error + ierrcd + if (associated(eccdpar)) deallocate(eccdpar) + + ! current drive efficiency [A⋅m/W] + effjcdav = rbavi*effjcd + didp = sgnbphi * effjcdav / (2.0_wp_*pi*rrii) + end subroutine alpha_effj - subroutine set_pol(ywrk0,bres,sox,psipol0,chipol0,ext0,eyt0) - use const_and_precisions, only : 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 + subroutine set_pol(ywrk0, bres, sox, ipol, psipol0, chipol0, ext0, eyt0) + use const_and_precisions, only : degree, zero, one, half, im + use beamdata, only : nray,nrayth + use equilibrium, only : bfield + 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 + + ! subroutine arguments + real(wp_), dimension(6, nray), intent(in) :: ywrk0 + real(wp_), intent(in) :: bres, sox + integer, intent(in) :: ipol + 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 @@ -1646,7 +1731,7 @@ contains csphi=xmv(1)/rm snphi=xmv(2)/rm call bfield(rm,xmv(3),bphi,br,bz) -! bv(i) = B_i in cartesian coordinates + ! bv(i) = B_i in cartesian coordinates bv(1)=br*csphi-bphi*snphi bv(2)=br*snphi+bphi*csphi bv(3)=bz @@ -1697,7 +1782,7 @@ contains px = 0.5_wp_ - a = reshape(matr2dgrid,(/nr*nz/)) + a = reshape(matr2dgrid, [nr*nz]) rcon = 0.0_wp_ zcon = 0.0_wp_ @@ -1767,7 +1852,7 @@ contains end do if(ix<=0) return - + bb: do in=ix jx=lx(in) @@ -1842,7 +1927,7 @@ bb: do itm=itm+1 inext= i if(jfor/=0) exit aa - if(itm .gt. 3) then + if(itm > 3) then flag1=.true. exit aa end if @@ -1853,8 +1938,8 @@ bb: do if(.not.flag1) then lx(in)=0 - if(itm .eq. 1) exit - end if + if(itm == 1) exit + end if jfor=jx jx=lx(inext) @@ -1989,13 +2074,13 @@ bb: do real(wp_), dimension(nq,nq) :: btotal if (.not. unit_active(ubres)) return - + 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 + end do ! Btotal on psi grid btmx = -1.0e30_wp_ @@ -2058,7 +2143,7 @@ bb: do 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) @@ -2301,7 +2386,7 @@ bb: do 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 do write(upec, *) '' end subroutine print_pec diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 99b631f..fa27bd6 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -51,7 +51,7 @@ module gray_params type raytracing_parameters real(wp_) :: dst ! Integration step size real(wp_) :: rwmax ! Normalized maximum radius of beam power - integer :: nrayr, nrayth ! Radial/angular number of rays + integer :: nray, nrayr, nrayth ! Total/radial/angular number of rays integer :: igrad ! Complex eikonal switch integer :: nstep ! Max number of integration steps integer :: idst ! Choice of the integration variable diff --git a/src/main.f90 b/src/main.f90 index e87c26e..c1a9102 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -408,8 +408,7 @@ contains use const_and_precisions, only : zero, degree use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff use dispersion, only : expinit - use gray_params, only : gray_parameters, print_parameters, & - headw, headl + use gray_params, only : gray_parameters, print_parameters use beams, only : launchangles2n, xgygcoeff use magsurf_data, only : flux_average, dealloc_surfvec use beamdata, only : init_btr, dealloc_beam @@ -421,12 +420,12 @@ contains implicit none ! subroutine arguments - type(gray_parameters), intent(in) :: params - real(wp_), intent(in) :: pabs, icd - real(wp_), dimension(:), intent(in) :: jphi, jcd, dpdv, currins, pins + type(gray_parameters), intent(inout) :: params + 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 + ratjamx, ratjbmx ! local variables real(wp_) :: ak0, bres, xgcn