! This modules contains the core GRAY routines module gray_core use const_and_precisions, only : wp_ implicit none contains subroutine gray_main(params, equil, plasma, limiter, results, error) use const_and_precisions, only : zero, one, comp_tiny use polarization, only : ellipse_to_field use types, only : contour, table, wrap use gray_params, only : gray_parameters, gray_results, EQ_VACUUM, ABSORP_OFF use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma use gray_project, only : ray_projector use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, & store_beam_shape, find_flux_surfaces, & flux_surfaces, kinetic_profiles, flux_averages, & ec_resonance, inputs_maps use gray_errors, only : gray_error, is_critical, has_new_errors, & print_err_raytracing, print_err_ecrh_cd use beams, only : xgygcoeff, launchangles2n use beamdata, only : pweight use utils, only : vmaxmin use multipass, only : initbeam, initmultipass, turnoffray, & plasma_in, plasma_out, wall_out use logger, only : log_info, log_debug ! subroutine arguments type(gray_parameters), intent(inout) :: params class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma type(contour), intent(in) :: limiter type(gray_results), intent(out) :: results integer(kind=gray_error), intent(out) :: error ! local variables real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) character, dimension(2), parameter :: mode=['O','X'] integer :: sox real(wp_) :: 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 integer :: iox,nharm,nhf,iokhawa,index_rt, parent_index_rt integer :: ip,ib,iopmin,ipar,child_index_rt integer :: igrad_b,istop_pass,nbeam_pass logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff ! Error handing integer(kind=gray_error) :: prev_error ! from previous step integer(kind=gray_error) :: curr_error ! current step ! i: integration step, jk: global ray index integer :: i, jk ! buffer for formatting log messages character(256) :: msg type(ray_projector) :: projector real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg associate (nray => params%raytracing%nray, & nrayr => params%raytracing%nrayr, & nrayth => params%raytracing%nrayth, & nstep => params%raytracing%nstep, & npass => params%raytracing%ipass, & nbeam_tot => 2**(params%raytracing%ipass+1)-2, & nbeam_max => 2**(params%raytracing%ipass)) block ! ray variables real(wp_), dimension(6, nray) :: yw, ypw real(wp_), dimension(3, nray) :: gri real(wp_), dimension(3, 3, nray) :: ggri real(wp_), dimension(3, nrayth, nrayr) :: xc, du1 real(wp_), dimension(nray) :: ccci0, p0jk real(wp_), dimension(nray) :: tau0, alphaabs0 real(wp_), dimension(nray) :: dids0 integer, dimension(nray) :: iiv real(wp_), dimension(nray, nstep) :: psjki, ppabs, ccci complex(wp_), dimension(nray) :: ext, eyt ! multipass variables logical, dimension(nray) :: iwait integer, dimension(nray) :: iow, iop logical, dimension(nray, nbeam_tot) :: iroff real(wp_), dimension(6, nray, nbeam_max-2) :: yynext, yypnext real(wp_), dimension(6, nray) :: yw0, ypw0 real(wp_), dimension(nray, nbeam_tot) :: stnext real(wp_), dimension(nray) :: stv, p0ray real(wp_), dimension(nray, nbeam_tot) :: taus, cpls real(wp_), dimension(nray) :: tau1, etau1, cpl1, lgcpl1 real(wp_), dimension(0:nbeam_tot) :: psipv, chipv ! beam variables real(wp_), dimension(params%output%nrho) :: jphi_beam, pins_beam real(wp_), dimension(params%output%nrho) :: currins_beam, dpdv_beam, jcd_beam real(wp_), dimension(2) :: mode_ellipse ! ψ, χ of the current polarisation mode ! ======== set environment BEGIN ======== ! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1) call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn) ! Initialise the output tables call init_tables(params, results%tables) ! Compute the initial cartesian wavevector (anv0) ! from launch angles α,β and the position call launchangles2n(params%antenna, anv0) if (params%equilibrium%iequil /= EQ_VACUUM) then ! Initialise the output profiles call projector%init(params%output, equil) end if ! Allocate memory for the results... allocate(results%dpdv(params%output%nrho)) allocate(results%jcd(params%output%nrho)) ! ...and initialise them results%pabs = 0 results%icd = 0 results%dpdv = 0 results%jcd = 0 ! ========= set environment END ========= ! Pre-determinted tables results%tables%kinetic_profiles = kinetic_profiles(params, equil, plasma) results%tables%flux_averages = flux_averages(params, equil) results%tables%flux_surfaces = flux_surfaces(params, equil) results%tables%ec_resonance = ec_resonance(params, equil, bres) results%tables%inputs_maps = inputs_maps(params, equil, plasma, bres, xgcn) ! print ψ rational surfaces call find_flux_surfaces(qvals=[1.0_wp_, 1.5_wp_, 2.0_wp_], & equil=equil, tbl=results%tables%flux_surfaces) ! print initial position write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos call log_info(msg, mod='gray_core', proc='gray_main') write (msg, '("initial direction:",2(x,a,"=",g0.2))') & 'α', params%antenna%alpha, 'β', params%antenna%beta call log_info(msg, mod='gray_core', proc='gray_main') ! =========== main loop BEGIN =========== call initmultipass(params, iroff, yynext, yypnext, yw0, ypw0, & stnext, p0ray, taus, tau1, etau1, cpls, & cpl1, lgcpl1, psipv, chipv) sox=1 ! mode inverted for each beam iox=2 ! start with O: sox=-1, iox=1 call pweight(params%raytracing, params%antenna%power, p0jk) ! Set original polarisation psipv(0) = params%antenna%psi chipv(0) = params%antenna%chi call ellipse_to_field(psipv(0), chipv(0), ext, eyt) ! Error value for the whole simulation error = 0 nbeam_pass=1 ! max n of beam per pass 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 main_loop: do ip=1,params%raytracing%ipass write (msg, '("pass: ",g0)') ip call log_info(msg, mod='gray_core', proc='gray_main') pabs_pass = 0 icd_pass = 0 istop_pass = 0 ! stop flag for current pass nbeam_pass = 2*nbeam_pass ! max n of beams in current pass if(ip > 1) then du1 = 0 gri = 0 ggri = 0 if(ip == params%raytracing%ipass) cpl = [zero, zero] ! no successive passes end if ! =========== beam loop BEGIN =========== call log_debug('beam loop start', mod='gray_core', proc='gray_main') 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 child_index_rt = 2*index_rt + 1 ! * index_rt of O-mode child beam parent_index_rt = ceiling(index_rt / 2.0_wp_) - 1 ! * index_rt of parent beam call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, & pins_beam,currins_beam,dpdv_beam,jcd_beam) write(msg, '(" beam: ",g0," (",a1," mode)")') index_rt, mode(iox) call log_info(msg, mod='gray_core', proc='gray_main') if(iboff) then ! no propagation for current beam istop_pass = istop_pass +1 ! * +1 non propagating beam call log_info(" beam is off", mod='gray_core', proc='gray_main') cycle end if psjki = 0 ppabs = 0 ccci = 0 tau0 = 0 alphaabs0 = 0 dids0 = 0 ccci0 = 0 iiv = 1 if(ip == 1) then ! 1st pass igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass tau1 = 0 ! * tau from previous passes etau1 = 1 cpl1 = 1 ! * coupling from previous passes lgcpl1 = 0 p0ray = p0jk ! * initial beam power call compute_initial_conds(params%raytracing, params%antenna, & ! * initial conditions anv0, ak0, yw, ypw, stv, xc, du1, gri, ggri) do jk=1,params%raytracing%nray ! Save the step "zero" data call store_ray_data(params, equil, results%tables, & i=0, jk=jk, s=stv(jk), P0=p0jk(jk), pos=yw(1:3,jk), & psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, & N_pl=zero, N_pr=zero, N=yw(4:6,jk), N_pr_im=zero, & n_e=zero, T_e=zero, alpha=zero, tau=zero, dIds=zero, & nhm=0, nhf=0, iokhawa=0, index_rt=index_rt, & lambda_r=zero, lambda_i=zero, X=zero, Y=zero, & grad_X=[zero,zero,zero]) 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 (limiter%contains(rrm, zzm)) then ! * start propagation in/outside vessel? iow(jk) = 1 ! + inside else 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 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(params%raytracing%nray>1 .and. all(.not.iwait)) & ! iproj=0 ==> nfilp=8 call store_beam_shape(params%raytracing, results%tables, stv, yw) ! ======= propagation loop BEGIN ======= call log_debug(' propagation loop start', mod='gray_core', proc='gray_main') ! error value for the this step curr_error = 0 propagation_loop: do i=1,params%raytracing%nstep ! rotate current and previous step errors prev_error = curr_error curr_error = 0 ! 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) + params%raytracing%dst ! current ray step call rkstep(params, equil, plasma, 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(params%raytracing, yw, ak0, xc, & du1, gri, ggri, curr_error) if (has_new_errors(prev_error, curr_error)) then call print_err_raytracing(curr_error, i, jk) end if iopmin = 10 ! =========== 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 xv = yw(1:3,jk) anv = yw(4:6,jk) call ywppla_upd(params, equil, plasma, & xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, & xgcn,ypw(:,jk), psinv, dens, btot, bv, xg, yg, & derxg, anpl, anpr, ddr, ddi, dersdst, derdnm, & igrad_b) ! 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 continue current pass if(cpl(iox) < etaucr) then ! + IF low coupled power for current mode => de-activate derived rays call turnoffray(jk,ip+1,npass,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 == 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') write (msg,'(" polarisation: ψ=",g0.5,"°, χ=",g0.5,"°")') & psipol, chipol call log_debug(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) > 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 < 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) - 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,npass,2*ib-1,iroff) if(cpl(1).le.comp_tiny) cpl(1)=etaucr end if if(cpl(2) < etaucr) then ! . low coupled power for X-mode => de-activate derived rays call turnoffray(jk,ip+1,npass,2*ib,iroff) if(cpl(2).le.comp_tiny) cpl(2)=etaucr end if taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass cpls(jk,child_index_rt+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass if(jk == 1) then ! . polarization angles at plasma boundary for central ray psipv(child_index_rt:child_index_rt+1) = psipol chipv(child_index_rt:child_index_rt+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] end if end if else if(ext_pl) then ! ray exits plasma write (msg, '(" ray ", g0, " left plasma")') jk call log_debug(msg, mod='gray_core', proc='gray_main') call plasma_out(jk, equil, 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, equil, limiter, ins_pl, xv, anv, & params%raytracing%dst, 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(params, equil, plasma, & xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, & xgcn, ypw(:,jk), psinv, dens, btot, bv, xg, & yg, derxg, anpl, anpr, ddr, ddi, dersdst, & derdnm, igrad_b) ! * update derivatives after reflection 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 < 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, equil, xv, anv, bres, sox, cpl, & ! . ray re-enters plasma after reflection psipol, chipol, iop, ext, eyt, perfect=.false.) if(cpl(1) < etaucr) then ! . low coupled power for O-mode? => de-activate derived rays call turnoffray(jk,ip+1,npass,2*ib-1,iroff) if(cpl(1).le.comp_tiny) cpl(1)=etaucr end if if(cpl(2) < etaucr) then ! . low coupled power for X-mode? => de-activate derived rays call turnoffray(jk,ip+1,npass,2*ib,iroff) if(cpl(2).le.comp_tiny) cpl(2)=etaucr end if taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass cpls(jk,child_index_rt+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass if(jk == 1) then ! + polarization angles at plasma boundary for central ray psipv(child_index_rt:child_index_rt+1) = psipol chipv(child_index_rt:child_index_rt+1) = chipol end if end if end if end if iopmin = min(iopmin,iop(jk)) 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) ! Note: power check is τ + τ₀ + log(coupling O/X) < τ_critic if (params%ecrh_cd%iwarm /= ABSORP_OFF .and. ins_pl .and. & (tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check tekev = plasma%temp(psinv) block complex(wp_) :: Npr_warm call alpha_effj(params%ecrh_cd, equil, plasma, & psinv, xg, yg, dens, tekev, ak0, bres, & derdnm, anpl, anpr, sox, Npr_warm, alpha, & didp, nharm, nhf, iokhawa, curr_error) anprre = Npr_warm%re anprim = Npr_warm%im if (has_new_errors(prev_error, curr_error)) then call print_err_ecrh_cd(curr_error, i, jk, Npr_warm, alpha) end if end block else tekev=0 alpha=0 didp=0 anprim=0 anprre=anpr nharm=0 nhf=0 iokhawa=0 end if if(nharm>0) iiv(jk)=i psjki(jk,i) = psinv ! 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 driven: dI/ds = dI/dP⋅dP/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: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 store_ray_data(params, equil, results%tables, & 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 ! Accumulate errors from the latest step error = ior(error, curr_error) 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,npass,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,params%output%istpr) == 0) then if(params%raytracing%nray > 1 .and. all(.not.iwait)) & call store_beam_shape(params%raytracing, results%tables, stv, yw) end if ! test whether further trajectory integration is unnecessary call vmaxmin(tau1+tau0+lgcpl1, taumn, taumx) ! test on tau + coupling if (is_critical(curr_error)) then ! stop propagation for current beam istop_pass = istop_pass +1 ! * +1 non propagating beam if(ip < params%raytracing%ipass) call turnoffray(0,ip,npass,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 propagation_loop call log_debug(' propagation loop end', mod='gray_core', proc='gray_main') ! ======== propagation loop END ======== ! print all ray positions in local reference system if(params%raytracing%nray > 1 .and. all(.not.iwait)) & call store_beam_shape(params%raytracing, results%tables, & stv, yw, full=.true.) ! =========== 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, taumn, taumx) ! taumn,taumx for print if (params%equilibrium%iequil /= EQ_VACUUM) then ! compute power and current density profiles for all rays call projector%project(psjki, ppabs, ccci, iiv, dpdv_beam, & jphi_beam, jcd_beam, pins_beam, currins_beam) end if 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 < 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(:,child_index_rt)/cpl1, MASK=iop > 2) / & sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays) cpl_beam2 = 1 - cpl_beam1 ! * average X-mode coupling for next beam if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam cpl_cbeam1 = cpls(1,child_index_rt)/cpl1(1) cpl_cbeam2 = 1 - cpl_cbeam1 end if else ! last pass OR no ray re-entered plasma cpl_beam1 = 0 cpl_beam2 = 0 end if ! print final results for pass on screen call log_info(' partial results:', mod='gray_core', proc='gray_main') write(msg, '(" final step:", x,a,g0, x,a,g0.4)') 'n=', i, '(s, ct, Sr)=', stv(1) call log_info(msg, mod='gray_core', proc='gray_main') write(msg, '(3x,a,2(x,a,"=",g0.4))') 'optical depth:', 'τ_min', taumn, 'τ_max', taumx call log_info(msg, mod='gray_core', proc='gray_main') write(msg, '(3x,a,g0.3," MW")') 'absoption: P=', pabs_beam call log_info(msg, mod='gray_core', proc='gray_main') write(msg, '(3x,a,g0.3," kA")') 'current drive: I=', icd_beam*1.0e3_wp_ call log_info(msg, mod='gray_core', proc='gray_main') 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') if(iop(1) > 2) then write(msg, '(3x,a,(g0.4,", ",g0.4))') & 'coupling [ctr ray, O/X]:', cpl_cbeam1, cpl_cbeam2 ! central ray coupling for next O/X beams call log_info(msg, mod='gray_core', proc='gray_main') end if end if if (params%equilibrium%iequil /= EQ_VACUUM) then call store_ec_profiles( & results%tables%ec_profiles, projector%rho_p, projector%rho_t, & jphi_beam, jcd_beam, dpdv_beam, currins_beam, pins_beam, index_rt) call projector%statistics(equil, dpdv_beam, jphi_beam, & rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) ! *compute profiles width for current beam if (results%tables%summary%active) then call results%tables%summary%append([ & wrap(icd_beam), wrap(pabs_beam), wrap(jphip), wrap(dpdvp), & wrap(rhotj), wrap(rhotjava), wrap(rhotp), wrap(rhotpav), & wrap(drhotjava), wrap(drhotpav), wrap(ratjamx), wrap(ratjbmx), & wrap(stv(1)), wrap(mode_ellipse(1)), wrap(mode_ellipse(2)), & wrap(index_rt), wrap(jphimx), wrap(dpdvmx), wrap(drhotj), & wrap(drhotp), wrap(sum(p0ray)), wrap(cpl_beam1), wrap(cpl_beam2)]) ! *print 0D results for current beam end if end if ! ============ 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 ======= results%pabs = results%pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output] results%icd = results%icd + sum(icd_pass) ! print final results for pass on screen call log_info(' comulative results:', mod='gray_core', proc='gray_main') write(msg, '(" absoption [O,X mode] P=",g0.4,", ",g0.4," MW")') & pabs_pass(1), pabs_pass(2) call log_info(msg, mod='gray_core', proc='gray_main') 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 main_loop call log_debug('pass loop end', mod='gray_core', proc='gray_main') ! ============ main loop END ============ end block end associate end subroutine gray_main subroutine compute_initial_conds(rtx, beam, N_c, k0, y, yp, dist, & pos, grad_u, grad, hess) ! Computes the initial conditions for tracing a beam use const_and_precisions, only : zero, one, pi, im, degree use gray_params, only : raytracing_parameters, antenna_parameters use gray_params, only : STEP_ARCLEN, STEP_TIME, STEP_PHASE use utils, only : diag, rotate use beams, only : gaussian_eikonal use beamdata, only : fold_indices, tokamak2beam ! subroutine arguments ! Inputs type(raytracing_parameters), intent(in) :: rtx ! simulation parameters type(antenna_parameters), intent(in) :: beam ! beam parameters real(wp_), intent(in) :: N_c(3) ! vacuum refractive index real(wp_), intent(in) :: k0 ! vacuum wavenumber ! Outputs real(wp_), intent(out) :: y(6, rtx%nray) ! (x̅, N̅) for each ray real(wp_), intent(out) :: yp(6, rtx%nray) ! (dx̅/dσ, dN̅/dσ) for each ray real(wp_), intent(out) :: dist(rtx%nrayth * rtx%nrayr) ! distance travelled in σ for each ray real(wp_), intent(out) :: pos(3, rtx%nrayth, rtx%nrayr) ! ray positions (used in gradi_upd) real(wp_), intent(out) :: grad_u(3, rtx%nrayth, rtx%nrayr) ! variations Δu (used in gradi_upd) real(wp_), intent(out) :: grad(3, rtx%nray) ! ∇S_I, gradient of imaginary eikonal real(wp_), intent(out) :: hess(3, 3, rtx%nray) ! HS_I, Hessian of imaginary eikonal ! local variables real(wp_) :: x0_c(3), x0(3) ! beam centre, current ray position real(wp_) :: xi(2) ! position in the principal axes basis (ξ, η) integer :: j, k, jk ! ray indices, j: radial, k: angular real(wp_) :: dr, da ! ray grid steps, dρ: radial, dα: angular real(wp_) :: N(3) ! wavevector complex(wp_) :: Q(2,2) ! complex cuvature tensor real(wp_) :: M(3,3) ! local beam → tokamak change of basis matrix complex(wp_) :: S ! complex eikonal complex(wp_) :: grad_S(3) ! its gradient complex(wp_) :: hess_S(3,3) ! its Hessian matrix real(wp_) :: phi_k, phi_w ! ellipses rotation angles ! Compute the complex curvature tensor Q from the beam parameters block real(wp_) :: K(2,2), W(2,2) K = diag(beam%ri) ! curvature part W = diag(2/k0 * 1/beam%w**2) ! beam width part ! Switch from eigenbasis to local beam basis ! ! Notes on geometry: ! ! 1. The beam parameters are the eigenvalues of the K, W matrices. ! These are diagonal in the basis of the principal axes of the ! respective ellipse (curvature, beam width). ! ! 2. The φs are nothing more that the angles of the rotations that ! diagonalises the respective matrix. ! phi_k = beam%phi(2) * degree ! curvature rotation angle phi_w = beam%phi(1) * degree ! beam width rotation angle K = matmul(rotate(phi_k), matmul(K, rotate(-phi_k))) W = matmul(rotate(phi_w), matmul(W, rotate(-phi_w))) ! Combine into a single tensor Q = K - im*W end block ! Compute the beam → tokamak change-of-basis matrix M = tokamak2beam(N_c, inverse=.true.) ! Expanding the definitions (see `gaussian_eikonal`), the imaginary ! part of the eikonal becomes: ! ! S_I = Im(z + ½ r̅⋅Q(z)⋅r̅) ! = ½ r̅⋅ImQ(z)⋅r̅ ! = -½ r̅⋅W⋅r̅ ! = -½ 2/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅ ! = -1/k₀ ξ̅⋅diag(1/w_ξ², 1/w_η²)⋅ξ̅ ! ! Using polar coordinates ξ̅ = ρ [w_ξcosα, w_ηsinα] this ! simplifies to S_I = -ρ²/k₀. To distribue the rays uniformly ! on the curves of constant S_I (power flux) we define a grid ! such that the jk-th ray is traced starting from ! ! ξ̅_jk = ρ_j [w_ξcos(α_k), w_ηcos(α_k)] ! ! where ρ_j = (j-1) Δρ, ! α_k = (k-1) Δα, ! Δρ = ρ_max/(N_r - 1), ! Δα = 2π/N_θ. ! dr = merge(rtx%rwmax / (rtx%nrayr - 1), one, rtx%nrayr > 1) da = 2*pi / rtx%nrayth ! To estimate ∇S_I in `gradi_upd` we require the values of ! ∇u(ξ̅_jk), where u(x̅) ≡ ρ(x̅)/Δρ and ξ̅_jk is the ray ! position on the grid. Since S_I = -ρ²/k₀, we have ! ! ∇S_I(ξ̅_jk) = -(Δρ²/k₀) ∇u²(ξ̅_jk) ! = -(Δρ²/k₀) 2u(ξ̅_jk) ∇u(ξ̅_jk) ! = -(Δρ²/k₀) 2(ρ_j/Δρ) ∇u(ξ̅_jk) ! = -(Δρ²/k₀) 2(j-1) ∇u(ξ̅_jk) ! ! So, ∇u = k₀/[2Δρ²(1-j)] ∇S_I, for j>1. ! For the central ray (j=1) r=0, so ∇S_I=HS_I=0 jk = 1 grad(:,1) = 0 hess(:,:,1) = 0 x0_c = beam%pos y(:,1) = [x0_c, N_c] yp(:,1) = [N_c, 0*N_c] ! We consider the central ray as k degenerate rays ! with positions x̅=x₀_c, but ∇u computed at the points ! ξ̅_1k = Δρ [w_ξcos(α_k), w_ηcos(α_k)]. This choice ! simplifies the computation of ∇u in `gradi_upd`. do k = 1, rtx%nrayth xi = beam%w * dr * [cos((k-1)*da), sin((k-1)*da)] x0 = [matmul(rotate(phi_w), xi), zero] pos(:,k,1) = x0_c grad_u(:,k,1) = [-matmul(Q%im, x0(1:2) * k0/(2*dr**2)), zero] end do ! Compute x̅, ∇S_I, ∇u for all the other rays do concurrent (j=2:rtx%nrayr, k=1:rtx%nrayth) jk = fold_indices(rtx, j, k) ! Compute the ray position xi = beam%w * (j-1)*dr * [cos((k-1)*da), sin((k-1)*da)] x0 = [matmul(rotate(phi_w), xi), zero] ! Compute the Eikonal and its derivatives call gaussian_eikonal(Q, r=x0(1:2), z=x0(3), & S=S, grad_S=grad_S, hess_S=hess_S) ! Compute ∇u using the formula above grad_u(:,k,j) = matmul(M, k0/(2 * dr**2 * (1-j)) * grad_S%im) ! Switch from local beam to tokamak basis x0 = x0_c + matmul(M, x0) ! ray position grad(:,jk) = matmul(M, grad_S%im) ! gradient of S_I hess(:,:,jk) = matmul(M, matmul(hess_S%im, transpose(M))) ! Hessian of S_I ! Compute the refractive index N such that: ! ! 1. N² = 1 + ∇S_I² ! 2. N⋅∇S_I = 0 ! 3. N₁/N₂ = ∇S_R₁/∇S_R₂ ! ! Note: 1. and 2. are the quasioptical dispersion relation, ! 3. is necessary to obtain a unique solution. block real(wp_) :: den, r1, r2, gradS2 den = dot_product(grad_S(1:2)%re, grad_S(1:2)%im) gradS2 = norm2(grad_S%im)**2 if (den /= 0) then r1 = -grad_S(1)%re * grad_S(3)%im / den r2 = -grad_S(2)%re * grad_S(3)%im / den N = [r1, r2, one] * sqrt((1 + gradS2)/(1 + r1**2 + r2**2)) else if (grad_S(3)%im /= 0) then ! In this case the solution reduces to r1 = grad_S(1)%re r2 = grad_S(2)%re N = [r1, r2, zero] * sqrt((1 + gradS2)/(1 + r1**2 + r2**2)) else ! When even ∇S_I₃ = 0, the system is underdetermined: ! we pick the solution N₃ = 1 + ∇S_I², so that N₁=N₂=0. N = [zero, zero, sqrt(1 + gradS2)] end if end block ! Compute the actual initial conditions pos(:,k,j) = x0 y(:,jk) = [x0, matmul(M, N)] ! Compute the r.h.s. of the raytracing eqs. ! dx̅/dσ = + ∂Λ/∂N̅ / denom ! dN̅/dσ = - ∂Λ/∂x̅ / denom block real(wp_) :: denom ! Compute travelled distance in σ and r.h.s. normalisation select case(rtx%idst) case (STEP_ARCLEN) ! σ=s, arclength dist(jk) = 0 denom = norm2(N) case (STEP_TIME) ! σ=ct, time dist(jk) = 0 denom = 1 case (STEP_PHASE) ! σ=S_R, phase dist(jk) = S%re denom = norm2(N)**2 end select yp(1:3,jk) = matmul(M, N) / denom yp(4:6,jk) = matmul(M, matmul(hess_S%im, grad_S%im)) / denom end block end do end subroutine compute_initial_conds subroutine rkstep(params, equil, plasma, & sox, bres, xgcn, y, yp, dgr, ddgr, igrad) ! Runge-Kutta integrator use gray_params, only : gray_parameters use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma ! subroutine arguments type(gray_parameters), intent(in) :: params class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma real(wp_), intent(in) :: bres, xgcn real(wp_), intent(inout) :: y(6) real(wp_), intent(in) :: yp(6) real(wp_), intent(in) :: dgr(3) real(wp_), intent(in) :: ddgr(3,3) integer, intent(in) :: igrad,sox ! local variables real(wp_), dimension(6) :: k1, k2, k3, k4 associate (h => params%raytracing%dst) k1 = yp k2 = f(y + k1*h/2) k3 = f(y + k2*h/2) k4 = f(y + k3*h) y = y + h * (k1/6 + k2/3 + k3/3 + k4/6) end associate contains function f(y) real(wp_), intent(in) :: y(6) real(wp_) :: f(6) call rhs(params, equil, plasma, sox, bres, xgcn, y, dgr, ddgr, f, igrad) end function end subroutine rkstep subroutine rhs(params, equil, plasma, & sox, bres, xgcn, y, dgr, ddgr, dery, igrad) ! Compute right-hand side terms of the ray equations (dery) ! used in R-K integrator use gray_params, only : gray_parameters use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma ! subroutine arguments type(gray_parameters), intent(in) :: params class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma real(wp_), intent(in) :: y(6) real(wp_), intent(in) :: bres, xgcn real(wp_), intent(in) :: dgr(3) real(wp_), intent(in) :: ddgr(3,3) real(wp_), intent(out) :: dery(6) integer, intent(in) :: igrad, sox ! local variables real(wp_) :: dens,btot,xg,yg real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg real(wp_), dimension(3,3) :: derbv xv = y(1:3) call plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & bv, derbv, xg, yg, derxg, deryg) anv = y(4:6) call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, & bv, derbv, dgr, ddgr, igrad, dery) end subroutine rhs subroutine ywppla_upd(params, equil, plasma, & xv, anv, dgr, ddgr, sox, bres, xgcn, dery, & psinv, dens, btot, bv, xg, yg, derxg, anpl, & anpr, ddr, ddi, dersdst, derdnm, igrad) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update use gray_params, only : gray_parameters use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma ! subroutine rguments type(gray_parameters), intent(in) :: params class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma real(wp_), intent(in) :: xv(3), anv(3) real(wp_), intent(in) :: dgr(3) real(wp_), intent(in) :: ddgr(3,3) integer, intent(in) :: sox real(wp_), intent(in) :: bres,xgcn real(wp_), intent(out) :: dery(6) real(wp_), intent(out) :: psinv, dens, btot, xg, yg, anpl, anpr real(wp_), intent(out) :: ddr, ddi, dersdst, derdnm real(wp_), intent(out) :: bv(3) real(wp_), intent(out) :: derxg(3) integer, intent(in) :: igrad ! local variables real(wp_):: deryg(3) real(wp_):: derbv(3,3) call plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & bv, derbv, xg, yg, derxg, deryg, psinv) call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, dgr, & ddgr, igrad, dery, anpl, anpr, ddr, ddi, dersdst, derdnm) end subroutine ywppla_upd subroutine gradi_upd(params, ywrk, ak0, xc, du1, gri, ggri, error) use gray_params, only : raytracing_parameters use gray_errors, only : gray_error, unstable_beam, raise_error ! subroutine parameters type(raytracing_parameters), intent(in) :: params real(wp_), intent(in) :: ak0 real(wp_), dimension(6,params%nray), intent(in) :: ywrk real(wp_), dimension(3,params%nrayth,params%nrayr), intent(inout) :: xc, du1 real(wp_), dimension(3,params%nray), intent(out) :: gri real(wp_), dimension(3,3,params%nray), intent(out) :: ggri integer(kind=gray_error), intent(inout) :: error ! local variables real(wp_), dimension(3,params%nrayth,params%nrayr) :: xco, du1o integer :: jk, j, jm, jp, k, km, kp real(wp_) :: ux, uxx, uxy, uxz, uy, uyy, uyz, uz, uzz real(wp_) :: dr, dfuu, dffiu, gx, gxx, gxy, gxz, gy, gyy, gyz, gz, gzz real(wp_), dimension(3) :: dxv1, dxv2, dxv3, dgu real(wp_), dimension(3,3) :: dgg, dff ! update position and du1 vectors xco = xc du1o = du1 jk = 1 do j=1,params%nrayr do k=1,params%nrayth if(j>1) jk=jk+1 xc(1:3,k,j)=ywrk(1:3,jk) end do end do ! compute grad u1 for central ray j = 1 jp = 2 do k=1,params%nrayth if(k == 1) then km = params%nrayth else km = k-1 end if if(k == params%nrayth) then kp = 1 else kp = k+1 end if dxv1 = xc(:,k ,jp) - xc(:,k ,j) dxv2 = xc(:,kp,jp) - xc(:,km,jp) dxv3 = xc(:,k ,j) - xco(:,k ,j) call solg0(dxv1,dxv2,dxv3,dgu) du1(:,k,j) = dgu end do gri(:,1) = 0 ! compute grad u1 and grad(S_I) for all the other rays if (params%nrayr > 1) then dr = params%rwmax / (params%nrayr - 1) else dr = 1 end if dfuu=2*dr**2/ak0 jm=1 j=2 k=0 dffiu = dfuu do jk=2,params%nray k=k+1 if(k > params%nrayth) then jm = j j = j+1 k = 1 dffiu = dfuu*jm end if kp = k+1 km = k-1 if (k == 1) then km=params%nrayth else if (k == params%nrayth) then kp=1 end if dxv1 = xc(:,k ,j) - xc(:,k ,jm) dxv2 = xc(:,kp,j) - xc(:,km,j) dxv3 = xc(:,k ,j) - xco(:,k ,j) call solg0(dxv1, dxv2, dxv3, dgu) ! Check for numerical instability in the computation of ∇S_I. ! This can happen when the distances between the rays and the ! integration step size differ wildly; in such cases the linear ! system Δx⋅∇u = Δu becomes bad conditioned. ! ! Source: https://discourse.julialang.org/t/what-is-the-best-way-to-solve-a-ill-conditioned-linear-problem/95894/10 block use const_and_precisions, only : comp_eps use utils, only : singvals real(wp_) :: A(3,3), x(3), b(3), sigma(3) real(wp_) :: R, normA, back_error, forw_error A = reshape([dxv1, dxv2, dxv3], shape=[3,3], order=[2,1]) sigma = singvals(A) b = [1, 0, 0] x = dgu R = maxval(sigma) / minval(sigma) ! condition number normA = maxval(sigma) ! L₂ norm of A back_error = norm2(matmul(A, x) - b) / (normA * norm2(b)) forw_error = back_error * R * norm2(x) ! TODO: figure out how to tie this threshold ! to the derivatives of the dispersion relation if (forw_error / comp_eps > 5000) then error = raise_error(error, unstable_beam) end if end block du1(:,k,j) = dgu gri(:,jk) = dgu(:)*dffiu end do ! compute derivatives of grad u and grad(S_I) for rays jk>1 ggri(:,:,1) = 0 jm=1 j=2 k=0 dffiu = dfuu do jk=2,params%nray k=k+1 if(k > params%nrayth) then jm=j j=j+1 k=1 dffiu = dfuu*jm end if kp=k+1 km=k-1 if (k == 1) then km=params%nrayth else if (k == params%nrayth) then kp=1 end if dxv1 = xc(:,k ,j) - xc(:,k ,jm) dxv2 = xc(:,kp,j) - xc(:,km,j) dxv3 = xc(:,k ,j) - xco(:,k ,j) dff(:,1) = du1(:,k ,j) - du1(:,k ,jm) dff(:,2) = du1(:,kp,j) - du1(:,km,j) dff(:,3) = du1(:,k ,j) - du1o(:,k ,j) call solg3(dxv1,dxv2,dxv3,dff,dgg) ! derivatives of u ux = du1(1,k,j) uy = du1(2,k,j) uz = du1(3,k,j) uxx = dgg(1,1) uyy = dgg(2,2) uzz = dgg(3,3) uxy = (dgg(1,2) + dgg(2,1))/2 uxz = (dgg(1,3) + dgg(3,1))/2 uyz = (dgg(2,3) + dgg(3,2))/2 ! derivatives of S_I and Grad(S_I) gx = ux*dffiu gy = uy*dffiu gz = uz*dffiu gxx = dfuu*ux*ux + dffiu*uxx gyy = dfuu*uy*uy + dffiu*uyy gzz = dfuu*uz*uz + dffiu*uzz gxy = dfuu*ux*uy + dffiu*uxy gxz = dfuu*ux*uz + dffiu*uxz gyz = dfuu*uy*uz + dffiu*uyz ggri(1,1,jk)=gxx ggri(2,1,jk)=gxy ggri(3,1,jk)=gxz ggri(1,2,jk)=gxy ggri(2,2,jk)=gyy ggri(3,2,jk)=gyz ggri(1,3,jk)=gxz ggri(2,3,jk)=gyz ggri(3,3,jk)=gzz end do end subroutine gradi_upd subroutine solg0(dxv1,dxv2,dxv3,dgg) ! solution of the linear system of 3 eqs : dgg . dxv = dff ! input vectors : dxv1, dxv2, dxv3, dff ! output vector : dgg ! dff=(1,0,0) ! arguments real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3 real(wp_), dimension(3), intent(out) :: dgg ! local variables real(wp_) :: denom,aa1,aa2,aa3 aa1 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3)) aa2 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3)) aa3 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3)) denom = dxv1(1)*aa1 - dxv2(1)*aa2 + dxv3(1)*aa3 dgg(1) = aa1/denom dgg(2) = -(dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3))/denom dgg(3) = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2))/denom end subroutine solg0 subroutine solg3(dxv1,dxv2,dxv3,dff,dgg) ! rhs "matrix" dff, result in dgg ! arguments real(wp_), dimension(3), intent(in) :: dxv1,dxv2,dxv3 real(wp_), dimension(3,3), intent(in) :: dff real(wp_), dimension(3,3), intent(out) :: dgg ! local variables real(wp_) denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 a11 = (dxv2(2)*dxv3(3) - dxv3(2)*dxv2(3)) a21 = (dxv1(2)*dxv3(3) - dxv3(2)*dxv1(3)) a31 = (dxv1(2)*dxv2(3) - dxv2(2)*dxv1(3)) a12 = (dxv2(1)*dxv3(3) - dxv3(1)*dxv2(3)) a22 = (dxv1(1)*dxv3(3) - dxv3(1)*dxv1(3)) a32 = (dxv1(1)*dxv2(3) - dxv2(1)*dxv1(3)) a13 = (dxv2(1)*dxv3(2) - dxv3(1)*dxv2(2)) a23 = (dxv1(1)*dxv3(2) - dxv3(1)*dxv1(2)) a33 = (dxv1(1)*dxv2(2) - dxv2(1)*dxv1(2)) denom = dxv1(1)*a11 - dxv2(1)*a21 + dxv3(1)*a31 dgg(:,1) = ( dff(:,1)*a11 - dff(:,2)*a21 + dff(:,3)*a31)/denom dgg(:,2) = (-dff(:,1)*a12 + dff(:,2)*a22 - dff(:,3)*a32)/denom dgg(:,3) = ( dff(:,1)*a13 - dff(:,2)*a23 + dff(:,3)*a33)/denom end subroutine solg3 subroutine plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & bv, derbv, xg, yg, derxg, deryg, psinv_opt) use const_and_precisions, only : cm use gray_equil, only : abstract_equil, vacuum use gray_plasma, only : abstract_plasma ! subroutine arguments class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma real(wp_), dimension(3), intent(in) :: xv real(wp_), intent(in) :: xgcn, bres real(wp_), intent(out) :: dens, btot, xg, yg real(wp_), dimension(3), intent(out) :: bv, derxg, deryg real(wp_), dimension(3,3), intent(out) :: derbv real(wp_), optional, intent(out) :: psinv_opt ! local variables integer :: jv real(wp_) :: xx,yy,zz real(wp_) :: psinv real(wp_) :: csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm real(wp_), dimension(3) :: dbtot,bvc real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv real(wp_) :: brr,bphi,bzz,dxgdpsi real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddensdpsi xg = 0 yg = 99._wp_ psinv = -1._wp_ dens = 0 btot = 0 derxg = 0 deryg = 0 bv = 0 derbv = 0 select type (equil) type is (vacuum) ! copy optional output if (present(psinv_opt)) psinv_opt = psinv return end select dbtot = 0 dbv = 0 dbvcdc = 0 dbvcdc = 0 dbvdc = 0 xx = xv(1) yy = xv(2) zz = xv(3) ! cylindrical coordinates rr2 = xx**2 + yy**2 rr = sqrt(rr2) csphi = xx/rr snphi = yy/rr bv(1) = -snphi*equil%sgn_bphi bv(2) = csphi*equil%sgn_bphi ! convert from cm to meters zzm = 1.0e-2_wp_*zz rrm = 1.0e-2_wp_*rr call equil%pol_flux(rrm, zzm, psinv, dpsidr, dpsidz, & ddpsidrr, ddpsidzz, ddpsidrz) call equil%pol_curr(psinv, fpolv, dfpv) ! copy optional output if (present(psinv_opt)) psinv_opt = psinv ! compute yg and derivative if(psinv < 0) then bphi = fpolv/rrm btot = abs(bphi) yg = btot/bres return end if ! compute xg and derivative call plasma%density(psinv, dens, ddensdpsi) xg = xgcn*dens dxgdpsi = xgcn*ddensdpsi ! B = f(psi)/R e_phi+ grad psi x e_phi/R bphi = fpolv/rrm brr = -dpsidz*equil%psi_a/rrm bzz = +dpsidr*equil%psi_a/rrm ! bvc(i) = B_i in cylindrical coordinates bvc = [brr, bphi, bzz] ! bv(i) = B_i in cartesian coordinates bv(1)=bvc(1)*csphi - bvc(2)*snphi bv(2)=bvc(1)*snphi + bvc(2)*csphi bv(3)=bvc(3) ! dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) dbvcdc(1,1) = -ddpsidrz*equil%psi_a/rrm - brr/rrm dbvcdc(1,3) = -ddpsidzz*equil%psi_a/rrm dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm dbvcdc(2,3) = dfpv*dpsidz/rrm dbvcdc(3,1) = +ddpsidrr*equil%psi_a/rrm - bzz/rrm dbvcdc(3,3) = +ddpsidrz*equil%psi_a/rrm ! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi dbvdc(2,1) = dbvcdc(1,1)*snphi + dbvcdc(2,1)*csphi dbvdc(3,1) = dbvcdc(3,1) dbvdc(1,2) = -bv(2) dbvdc(2,2) = bv(1) dbvdc(3,2) = dbvcdc(3,2) ! = 0 dbvdc(1,3) = dbvcdc(1,3)*csphi - dbvcdc(2,3)*snphi dbvdc(2,3) = dbvcdc(1,3)*snphi + dbvcdc(2,3)*csphi dbvdc(3,3) = dbvcdc(3,3) drrdx = csphi drrdy = snphi dphidx = -snphi/rrm dphidy = csphi/rrm ! dbv(iv,jv) = d Bcart(iv) / dxvcart(jv) dbv(:,1) = drrdx*dbvdc(:,1) + dphidx*dbvdc(:,2) dbv(:,2) = drrdy*dbvdc(:,1) + dphidy*dbvdc(:,2) dbv(:,3) = dbvdc(:,3) ! B magnitude and derivatives btot = norm2(bv) ! dbtot(i) = d |B| / dxvcart(i) dbtot = matmul(bv, dbv)/btot yg = btot/Bres ! convert spatial derivatives from dummy/m -> dummy/cm ! to be used in rhs ! bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z deryg = 1.0e-2_wp_*dbtot/Bres bv = bv/btot do jv=1,3 derbv(:,jv) = cm * (dbv(:,jv) - bv(:)*dbtot(jv))/btot end do derxg(1) = cm * drrdx*dpsidr*dxgdpsi derxg(2) = cm * drrdy*dpsidr*dxgdpsi derxg(3) = cm * dpsidz *dxgdpsi end subroutine plas_deriv subroutine disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, & dgr, ddgr, igrad, dery, anpl_, anpr, ddr, ddi, & dersdst, derdnm_) ! Computes the dispersion relation, derivatives and other ! related quantities ! ! The mandatory outputs are used for both integrating the ray ! trajectory (`rhs` subroutine) and updating the ray state ! (`ywppla_upd` suborutine); while the optional ones are used for ! computing the absoprtion and current drive. use const_and_precisions, only : zero, half use gray_params, only : gray_parameters, STEP_ARCLEN, & STEP_TIME, STEP_PHASE ! subroutine arguments ! Inputs type(gray_parameters), intent(in) :: params ! refractive index N̅ vector, b̅ = B̅/B magnetic field direction real(wp_), dimension(3), intent(in) :: anv, bv ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X integer, intent(in) :: sox ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω real(wp_), intent(in) :: xg, yg ! gradients of X, Y real(wp_), dimension(3), intent(in) :: derxg, deryg ! gradients of the complex eikonal real(wp_), dimension(3), intent(in) :: dgr ! ∇S_I real(wp_), dimension(3, 3), intent(in) :: ddgr ! ∇∇S_I ! gradient of the magnetic field direction real(wp_), dimension(3, 3), intent(in) :: derbv ! ∇b̅ ! raytracing/beamtracing switch integer, intent(in) :: igrad ! Outputs ! the actual derivatives: (∂Λ/∂x̅, -∂Λ/∂N̅) real(wp_), dimension(6), intent(out) :: dery ! additional quantities: ! refractive index real(wp_), optional, intent(out) :: anpl_, anpr ! N∥, N⊥ ! real and imaginary part of the dispersion real(wp_), optional, intent(out) :: ddr, ddi ! Λ, ∂Λ/∂N̅⋅∇S_I ! Jacobian ds/dσ, where s arclength, σ integration variable real(wp_), optional, intent(out) :: dersdst ! |∂Λ/∂N̅| real(wp_), optional, intent(out) :: derdnm_ ! Note: assign values to missing optional arguments results in a segfault. ! Since some are always needed anyway, we store them here and copy ! them later, if needed: real(wp_) :: anpl, derdnm ! local variables real(wp_) :: gr2, yg2, anpl2, del, dnl, duh, dan2sdnpl, an2, an2s real(wp_) :: dan2sdxg, dan2sdyg, denom real(wp_) :: derdom, dfdiadnpl, dfdiadxg, dfdiadyg, fdia, bdotgr real(wp_), dimension(3) :: derdxv, danpldxv, derdnv, dbgr an2 = dot_product(anv, anv) ! N² anpl = dot_product(anv, bv) ! N∥ = N̅⋅B̅ ! Shorthands used in the expressions below yg2 = yg**2 ! Y² anpl2 = anpl**2 ! N∥² dnl = 1 - anpl2 ! 1 - N∥² duh = 1 - xg - yg2 ! UH denom (duh=0 on the upper-hybrid resonance) ! Compute/copy optional outputs if (present(anpr)) anpr = sqrt(max(an2 - anpl2, zero)) ! N⊥ if (present(anpl_)) anpl_ = anpl ! N∥ an2s = 1 dan2sdxg = 0 dan2sdyg = 0 dan2sdnpl = 0 del = 0 fdia = 0 dfdiadnpl = 0 dfdiadxg = 0 dfdiadyg = 0 if(xg > 0) then ! Derivatives of the cold plasma refractive index ! ! N²s = 1 - X - XY²⋅(1 + N∥² ± √Δ)/[2(1 - X - Y²)] ! ! where Δ = (1 - N∥²)² + 4N∥²⋅(1 - X)/Y² ! + for the X mode, - for the O mode ! √Δ del = sqrt(dnl**2 + 4.0_wp_*anpl2*(1 - xg)/yg2) ! ∂(N²s)/∂X ! Note: this term is nonzero for X=0, but it multiplies terms ! proportional to X or ∂X/∂ψ which are zero outside the plasma. dan2sdxg = - half*yg2*(1 - yg2)*(1 + anpl2 + sox*del)/duh**2 & + sox*xg*anpl2/(del*duh) - 1 ! ∂(N²s)/∂Y dan2sdyg = - xg*yg*(1 - xg)*(1 + anpl2 + sox*del)/duh**2 & + 2*sox*xg*(1 - xg)*anpl2/(yg*del*duh) ! ∂(N²s)/∂N∥ dan2sdnpl = - xg*yg2*anpl/duh & - sox*xg*anpl*(2*(1 - xg) - yg2*dnl)/(del*duh) if(igrad > 0) then ! Derivatives used in the complex eikonal terms (beamtracing only) block real(wp_) :: ddelnpl2, ddelnpl2x, ddelnpl2y, derdel ! ∂²Δ/∂N∥² ddelnpl2 = 2*(2*(1 - xg)*(1 + 3.0_wp_*anpl2**2) & - yg2*dnl**3)/yg2/del**3 ! ∂²(N²s)/∂N∥² fdia = - xg*yg2*(1 + half*sox*ddelnpl2)/duh ! Intermediates results used right below derdel = 2*(1 - xg)*anpl2*(1 + 3*anpl2**2) & - dnl**2*(1 + 3*anpl2)*yg2 derdel = 4*derdel/(yg*del)**5 ddelnpl2y = 2*(1 - xg)*derdel ddelnpl2x = yg*derdel ! ∂³(N²s)/∂N∥³ dfdiadnpl = 24*sox*xg*(1 - xg)*anpl*(1 - anpl2**2)/(yg2*del**5) ! ∂³(N²s)/∂N∥²∂X dfdiadxg = - yg2*(1 - yg2)/duh**2 - sox*yg2*((1 - yg2) & *ddelnpl2 + xg*duh*ddelnpl2x)/(2*duh**2) ! ∂³(N²s)/∂N∥²∂Y dfdiadyg = - 2*yg*xg*(1 - xg)/duh**2 & - sox*xg*yg*(2*(1 - xg)*ddelnpl2 & + yg*duh*ddelnpl2y)/(2*duh**2) end block end if end if ! ∇(N∥) = ∇(N̅⋅b̅) = N̅⋅∇b̅ danpldxv = matmul(anv, derbv) ! ∂Λ/∂x̅ = - ∇(N²s) = -∂Λ/∂X ∇X - ∂Λ/∂Y ∇Y - ∂Λ/∂N∥ ∇(N∥) derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl) ! ∂Λ/∂N̅ = 2N̅ - ∂(N²s)/∂N̅ = 2N̅ - ∂(N²s)/∂N∥ b̅ ! Note: we used the identity ∇f(v̅⋅b̅) = f' ∇(v̅⋅b̅) = f'b̅. derdnv = 2*anv - dan2sdnpl*bv ! ∂Λ/∂ω = ∂N²/∂ω - ∂N²s/∂X⋅∂X/∂ω - ∂N²s/∂Y⋅∂Y/∂ω - ∂N²s/∂N∥⋅∂N∥/∂ω ! Notes: 1. N depends on ω: N²=c²k²/ω² ⇒ ∂N²/∂ω = -2N²/ω ! N∥=ck∥/ω ⇒ ∂N∥/∂ω = -N∥/ω ! 2. derdom is actually ω⋅∂Λ/∂ω, see below for the reason. ! 3. N gains a dependency on ω because Λ(∇S, ω) is computed ! on the constrains Λ=0. derdom = -2*an2 + 2*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl if (igrad > 0) then ! Complex eikonal terms added to the above expressions block real(wp_) :: dgr2(3) bdotgr = dot_product(bv, dgr) ! b̅⋅∇S_I gr2 = dot_product(dgr, dgr) ! |∇S_I|² dgr2 = 2 * matmul(ddgr, dgr) ! ∇|∇S_I|² = 2 ∇∇S_I⋅∇S_I ! ∇(b̅⋅∇S_I) = ∇b̅ᵗ ∇S_I + b̅ᵗ ∇∇S_I ! Notes: 1. We are using the convention ∇b̅_ij = ∂b_i/∂x_j. ! 2. Fortran doesn't distinguish between column/row vectors, ! so matmul(A, v̅) ≅ Av̅ and matmul(v̅, A) ≅ v̅ᵗA ≅ Aᵗv̅. dbgr = matmul(dgr, derbv) + matmul(bv, ddgr) ! ∂Λ/∂x̅ += - ∇(|∇S_I|²) + ½ ∇[(b̅⋅∇S_I)² ∂²N²s/∂N∥²] derdxv = derdxv - dgr2 + fdia*bdotgr*dbgr + half*bdotgr**2 & * (derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl) ! ∂Λ/∂N̅ += ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅ ! Note: we used again the identity ∇f(v̅⋅b̅) = f'b̅. derdnv = derdnv + half*bdotgr**2*dfdiadnpl*bv ! ∂Λ/∂ω += ∂|∇S_I|²/∂ω + ½∂(b⋅∇S_I)²/∂ω + ½(b̅⋅∇S_I)² ∂/∂ω (∂²N²s/∂N∥²) ! Note: as above ∇S_I gains a dependency on ω derdom = derdom + 2*gr2 - bdotgr**2 & * (fdia + xg*dfdiadxg + half*yg*dfdiadyg & + half*anpl*dfdiadnpl) end block end if derdnm = norm2(derdnv) ! Denominator of the r.h.s, depending on the ! choice of the integration variable σ: ! ! dx̅/dσ = + ∂Λ/∂N̅ / denom ! dN̅/dσ = - ∂Λ/∂x̅ / denom ! select case (params%raytracing%idst) case (STEP_ARCLEN) ! σ=s, arclength ! denom = |∂Λ/∂N̅| ! Proof: Normalising ∂Λ/∂N̅ (∥ to the group velocity) ! simply makes σ the arclength parameter s. denom = derdnm case (STEP_TIME) ! σ=ct, time ! denom = -ω⋅∂Λ/∂ω ! Proof: The Hamilton equations are ! ! dx̅/dt = + ∂H/∂k̅ (H=ℏΩ, p=ℏk̅) ! dp̅/dt = - ∂H/∂x̅ ! ! where Ω(x̅, k̅) is implicitely defined by the dispersion ! D(k̅, ω, x̅) = 0. By differentiating the latter we get: ! ! ∂Ω/∂x̅ = - ∂D/∂x̅ / ∂D/∂ω = - ∂Λ/∂x̅ / ∂Λ/∂ω ! ∂Ω/∂k̅ = - ∂D/∂k̅ / ∂D/∂ω = - ∂Λ/∂k̅ / ∂Λ/∂ω ! ! with Λ=D/k₀, k₀=ω/c. Finally, substituting k̅=k₀N̅: ! ! dx̅/d(ct) = + ∂Λ/∂N̅ / (-ω⋅∂Λ/∂ω) ! dN̅/d(ct) = - ∂Λ/∂x̅ / (-ω⋅∂Λ/∂ω) ! denom = -derdom case (STEP_PHASE) ! s=S_R, phase ! denom = N̅⋅∂Λ/∂N̅ ! Note: This parametrises the curve by the value ! of the (real) phase, with the wave frozen in time. ! Proof: By definition N̅ = k̅/k₀ = ∇S. Using the gradient ! theorem on the ray curve parametrised by the arclength ! ! S(s) = ∫₀^x̅(s) N̅⋅dl̅ = ∫₀^s N̅(σ)⋅dx̅/ds ds ! ! Differentiating gives dS = N̅⋅dx̅/ds ds, so ! ! dS = N̅⋅∂Λ/∂N̅ / |∂Λ/∂N̅| ds ⇒ ! ! dx̅/dS = + ∂Λ/∂N̅ / N̅⋅∂Λ/∂N̅ ! dN̅/dS = - ∂Λ/∂x̅ / N̅⋅∂Λ/∂N̅ ! denom = dot_product(anv, derdnv) end select ! Jacobian ds/dσ, where s is the arclength, ! for computing integrals in dσ, like: ! ! τ = ∫α(s)ds = ∫α(σ)(ds/dσ)dσ ! if (present(dersdst)) dersdst = derdnm/denom ! |∂Λ/∂N̅|: used in the α computation (see alpha_effj) if (present(derdnm_)) derdnm_ = derdnm ! r.h.s. vector dery(1:3) = derdnv(:)/denom ! +∂Λ/∂N̅ dery(4:6) = -derdxv(:)/denom ! -∂Λ/∂x̅ if (present(ddr) .or. present(ddi)) then ! Dispersion relation (geometric optics) ! ddr ← Λ = N² - N²s(X,Y,N∥) = 0 an2s = 1 - xg - half*xg*yg2*(1 + anpl2 + sox*del)/duh ddr = an2 - an2s end if if (present(ddr) .and. igrad > 0) then ! Dispersion relation (complex eikonal) ! ddr ← Λ = N² - N²s - |∇S_I|² + ½ (b̅⋅∇S_I)² ∂²N²s/∂N∥² ddr = ddr - gr2 - half*bdotgr**2*fdia ! real part end if ! Note: we have to return ddi even for igrad=0 because ! it's printed to udisp unconditionally if (present(ddi)) then ! Dispersion relation (complex eikonal) ! ddi ← ∂Λ/∂N̅⋅∇S_I ddi = merge(dot_product(derdnv, dgr), zero, igrad > 0) ! imaginary part end if end subroutine disp_deriv subroutine alpha_effj(params, equil, plasma, psinv, X, Y, density, & temperature, k0, Bres, derdnm, Npl, Npr_cold, & sox, Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error) ! Computes the absorption coefficient and effective current use const_and_precisions, only : pi, mc2=>mc2_ use gray_params, only : ecrh_cd_parameters use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma use dispersion, only : harmnumber, warmdisp use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl use gray_errors, only : gray_error, negative_absorption, raise_error ! subroutine arguments ! Inputs ! ECRH & CD parameters type(ecrh_cd_parameters), intent(in) :: params ! MHD equilibrium, plasma object class(abstract_equil), intent(in) :: equil class(abstract_plasma), intent(in) :: plasma ! poloidal flux ψ real(wp_), intent(in) :: psinv ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω real(wp_), intent(in) :: X, Y ! densityity [10¹⁹ m⁻³], temperature [keV] real(wp_), intent(in) :: density, temperature ! vacuum wavenumber k₀=ω/c, resonant B field real(wp_), intent(in) :: k0, Bres ! group velocity: |∂Λ/∂N̅| where Λ=c²/ω²⋅D_R real(wp_), intent(in) :: derdnm ! refractive index: N∥, N⊥ (cold dispersion) real(wp_), intent(in) :: Npl, Npr_cold ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X integer, intent(in) :: sox ! Outputs ! orthogonal refractive index N⊥ (solution of the warm dispersion) complex(wp_), intent(out) :: Npr ! 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(kind=gray_error), intent(inout) :: iokhawa, error ! local variables real(wp_) :: rbavi, rrii, rhop integer :: nlarmor, ithn, ierrcd real(wp_) :: mu, rbn, rbx real(wp_) :: zeff, cst2, bmxi, bmni, fci real(wp_), dimension(:), allocatable :: eccdpar real(wp_) :: effjcd, effjcdav, Btot complex(wp_) :: e(3) alpha = 0 Npr = 0 dIdp = 0 nhmin = 0 nhmax = 0 iokhawa = 0 ! Absorption computation ! Skip when the temperature is zero (no plasma) if (temperature <= 0) return ! Skip when the lowest resonant harmonic number is zero mu = mc2/temperature ! μ=(mc²/kT) call harmnumber(Y, mu, Npl**2, params%iwarm == 1, nhmin, nhmax) if (nhmin <= 0) return ! Solve the full dispersion only when needed if (params%iwarm /= 4 .or. params%ieccd /= 0) then nlarmor = max(params%ilarm, nhmax) if (params%ieccd /= 0) then ! Compute the polarisation vector only when current drive is on call warmdisp(X, Y, mu, Npl, Npr_cold, sox, error, Npr, e, & model=params%iwarm, nlarmor=nlarmor, & max_iters=abs(params%imx), fallback=params%imx < 0) else call warmdisp(X, Y, mu, Npl, Npr_cold, sox, error, Npr, & model=params%iwarm, nlarmor=nlarmor, & max_iters=abs(params%imx), fallback=params%imx < 0) end if end if ! Compute α from the solution of the dispersion relation ! 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̅ / ∂Λ/∂ω, using the cold dispersion ! relation, we have that ! ! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅| ! = [2N̅ - ∂(N²s)/∂N∥ b̅ ! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅| ! ! Assuming Im(k∥)=0: ! ! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅| ! block real(wp_) :: k_im k_im = k0 * Npr%im ! imaginary part of k⊥ alpha = 4 * k_im*Npr_cold / derdnm end block if (alpha < 0) then error = raise_error(error, negative_absorption) return end if ! Current drive computation if (params%ieccd <= 0) return zeff = plasma%zeff(psinv) ithn = 1 if (nlarmor > nhmin) ithn = 2 rhop = sqrt(psinv) call equil%flux_average(rhop, R_J=rrii, B_avg=rbavi, & B_min=bmni, B_max=bmxi, f_c=fci) rbavi = rbavi / bmni Btot = Y*Bres rbn = Btot/bmni rbx = Btot/bmxi select case(params%ieccd) case(1) ! Cohen model call setcdcoeff(zeff, rbn, rbx, cst2, eccdpar) call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & ithn, cst2, fjch, eccdpar, effjcd, iokhawa, ierrcd) case(2) ! No trapping call setcdcoeff(zeff, cst2, eccdpar) call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & ithn, cst2, fjch0, eccdpar, effjcd, iokhawa, ierrcd) case default ! Neoclassical model call setcdcoeff(zeff, rbx, fci, mu, rhop, equil%spline_cd_eff, cst2, eccdpar) call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & ithn, cst2, fjncl, eccdpar, effjcd, iokhawa, ierrcd) end select error = error + ierrcd ! current drive efficiency R* = / [A⋅m/W] effjcdav = rbavi*effjcd dIdp = equil%sgn_bphi * effjcdav / (2*pi*rrii) end subroutine alpha_effj end module gray_core