diff --git a/Makefile b/Makefile index 4eb8947..e49dd70 100644 --- a/Makefile +++ b/Makefile @@ -197,9 +197,9 @@ $(OBJDIR)/%.d: $(SRCDIR)/%.f90 | $(OBJDIR) $(OBJDIR)/%.o: $(SRCDIR)/%.f90 | $(OBJDIR) $(INCDIR) $(FC) -cpp $(CPPFLAGS) $(FFLAGS) -c '$<' -o '$@' -# Generate Fortran code from shell script -$(INCDIR)/%.inc: $(SRCDIR)/%.sh | $(INCDIR) - sh '$<' > '$@' +# Generate Fortran code from awk script +$(INCDIR)/%.inc: $(SRCDIR)/%.awk | $(INCDIR) + awk -f '$<' > '$@' # Create directories $(DIRS): diff --git a/input/gray.ini b/input/gray.ini index 60a0aee..8acec10 100644 --- a/input/gray.ini +++ b/input/gray.ini @@ -3,7 +3,7 @@ nrayr = 1 nrayth = 16 ; Normalized maximum radius of the beam power -rwmax = 1 +rwmax = 1.0 ; Switch between simple raytracing (0) and beamtracing (1) igrad = 0 @@ -13,29 +13,28 @@ igrad = 0 ; G-EQDSK file; when negative on a simple limiter at R=`rwall`, see below. ipass = 1 -; Whether to compute the wave polarisation -ipol = 0 +; Whether to compute the wave polarisation (from chi, psi) +ipol = .false. ; Step size (cm) for the numerical integration dst = 0.1 ; Max number of integration steps nstep = 12000 ; Choice of the integration variable -; 0: path length (s) -; 1: "time" (actually c⋅t) -; 2: real part of the eikonal (S_r) -idst = 0 +; STEP_ARCLEN: arc length (s) +; STEP_TIME: time (actually c⋅t) +; STEP_PHASE: phase (actually real part of eikonal S_r=k₀⋅φ) +idst = STEP_ARCLEN [ecrh_cd] ; Choice of the power absorption model -; 0: no absorption at all -; 1: weakly relativistic -; 2: fully relativistic (faster variant) -; 3: fully relativistic (slower variant) -; 4: tenuous plasma (very fast) -; Note: iwarm>0 is required for current driver -iwarm = 2 +; ABSORP_OFF: no absorption at all +; ABSORP_WEAK: weakly relativistic +; ABSORP_FULL: fully relativistic (faster variant) +; ABSORP_FULL_ALT: fully relativistic (slower variant) +; Note: iwarm /= ABSORP_OFF is required for current drive +iwarm = ABSORP_FULL ; Order of the electron Larmor radius expansion ; (used by some absorption models) @@ -48,11 +47,11 @@ ilarm = 5 imx = -20 ; Current drive model -; 0: no current drive at all -; 1: Cohen -; 2: no trapping -; 3: Neoclassical -ieccd = 3 +; CD_OFF: no current drive at all +; CD_COHEN: Cohen +; CD_NO_TRAP: no trapping +; CD_NEOCLASSIC: Neoclassical +ieccd = CD_NEOCLASSIC [antenna] @@ -61,22 +60,22 @@ alpha = 45 ; Poloidal angle (positive → up) beta = 0 ; Toroidal angle (positive → right) ; Injected power (MW) -power = 1 +power = 1.0 ; Polarisation mode -; 1: ordinary (O) -; 2: extraordinary (X) -iox = 1 +; MODE_O: ordinary (O) +; MODE_X: extraordinary (X) +iox = MODE_X ; Alternatively, parameters of the polarisation ellipse ; χ: angle between the principal axes and the (x,y) axes ; ψ: atan(ε), where ε is the ellipticity chi = 0 psi = 0 -; Beam description kind -; 0: simple beam shape -; 1: 1D table -; 2: 2D table +; Beam parameters format +; BEAM_0D: fixed beam parameters +; BEAM_1D: 1D steering angle table +; BEAM_2D: 2D steering angles table ibeam = 0 ; Filepath of the beam data (relative to this file) @@ -85,40 +84,40 @@ filenm = "beamdata.txt" [equilibrium] ; MHD equilibrium kind -; 0: vacuum (i.e. no plasma at all) -; 1: analytical -; 2: G-EQDSK format - data valid on the whole domain -; 3: G-EQDSK format - data valid only inside the LCFS -iequil = 3 +; EQ_VACUUM: vacuum (i.e. no plasma at all) +; EQ_ANALYTICAL: analytical model +; EQ_EQDSK_FULL: G-EQDSK format - data valid on the whole domain +; EQ_EQDSK_PARTIAL: G-EQDSK format - data valid only inside the LCFS +iequil = EQ_EQDSK_FULL ; Filepath of the equilibrium data (relative to this file) filenm = "magneticdata.eqdsk" ; COCOS index icocos = 0 -; Normalisation of the poloidal function -; 0: G-EQDSK, ψ → |ψ - ψ(edge)|/|ψ(axis) - ψ(edge)| -; 1: no normalisation -ipsinorm = 0 +; Whether the poloidal function is normalised (G-EQDSK) +; false: is not normalised, ψ → |ψ - ψ(edge)|/|ψ(axis) - ψ(edge)| +; true: is already normalised +ipsinorm = false ; G-EQDSK format parameters ; Whether header starts with a description, a.k.a identification string -idesc = 1 +idesc = true ; Whether the records have variable length ; Note: some non-compliant programs output numbers formatted with variable length ; instead of using the single (5e16.9) specifier. -ifreefmt = 0 +ifreefmt = false ; Position of the X point -; -1: bottom -; 0: no X point -; +1: top -ixp = 0 +; X_IS_MISSING: No X point +; X_AT_TOP: At the top of the plasma +; X_AT_BOTTOM: At the bottom of the plasma +ixp = X_IS_MISSING ; Tension of splines ; Note: 0 means perfect interpolation ssplps = 0.005 ; for ψ(R,Z), normalised poloidal flux -ssplf = 0.01 ; for I(ψ)=R⋅B_T, poloidal current function +ssplf = 0.01 ; for F(ψ)=R⋅B_T, poloidal current function ; Sign of toroidal field/current (used when COCOS index is 0,10) ; When viewing from above: +1 → counter-clockwise, -1 → clockwise @@ -131,16 +130,16 @@ factb = 1 [profiles] ; (input) plasma profiles parameters -; Profiles kind -; 0: analytical -; 1: numerical -iprof = 1 +; Plasma profiles kind +; PROF_ANALYTIC: analytical model +; PROF_NUMERIC: numerical data (ρ, n_e, T_e, table) +iprof = PROF_NUMERIC -; Profile radial coordinate -; 0: ρ_t = √Φ (where Φ is the normalised toroidal flux) -; 1: ρ_p = √ψ (where ψ is the normalised poloidal flux) -; 2: normalised poloidal flux ψ -irho = 0 +; Plasma profiles radial coordinate +; RHO_TOR: ρ_t = √Φ (where Φ is the normalised toroidal flux) +; RHO_POL: ρ_p = √ψ (where ψ is the normalised poloidal flux) +; RHO_PSI: normalised poloidal flux ψ +irho = RHO_TOR ; Filepath of the equilibrium (relative to this file) filenm = "profiles.txt" @@ -153,11 +152,12 @@ sspld = 0.1 factte = 1 factne = 1 -; Choice of model for rescaling the temperature/density with the -; magnetic field -; 1: costant Greenwald density (ab=1) -; 2: no rescaling (ab=0) -iscal = 2 +; Choice of model for rescaling the temperature/density +; with the magnetic field (if factb ≠ 0) +; SCALE_OFF: don't rescale at all +; SCALE_COLLISION: scale while preserving collisionality +; SCALE_GREENWALD: scale while preserving the Greenwald fraction +iscal = SCALE_OFF [output] @@ -165,8 +165,8 @@ iscal = 2 ; ECRH&CD profiles grid: ; Radial coordinate -; 1: ρ_p = √ψ (where ψ is the normalised poloidal flux) -; 2: ρ_t = √Φ (where Φ is the normalised toroidal flux) +; RHO_TOR: ρ_t = √Φ (where Φ is the normalised toroidal flux) +; RHO_POL: ρ_p = √ψ (where ψ is the normalised poloidal flux) ipec = 1 ; Number of points nrho = 501 diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index e45eebd..ee51ae8 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -114,7 +114,7 @@ contains end if ! get size of main arrays and allocate them - if (params%idesc == 1) then + if (params%idesc) then read (u,'(a48,3i4)') string,idum,nr,nz else read (u,*) nr, nz @@ -132,7 +132,7 @@ contains data%qpsi(nr)) ! Store 0D data and main arrays - if (params%ifreefmt==1) then + if (params%ifreefmt) then read (u, *) dr, dz, data%rvac, rleft, zmid read (u, *) data%rax, data%zax, psiaxis, psiedge, xdum read (u, *) xdum, xdum, xdum, xdum, xdum @@ -166,7 +166,7 @@ contains ! Load plasma boundary data if(nbnd > 0) then allocate(data%rbnd(nbnd), data%zbnd(nbnd)) - if (params%ifreefmt == 1) then + if (params%ifreefmt) then read(u, *) (data%rbnd(i), data%zbnd(i), i=1,nbnd) else read(u, '(5e16.9)') (data%rbnd(i), data%zbnd(i), i=1,nbnd) @@ -176,7 +176,7 @@ contains ! Load limiter data if(nlim > 0) then allocate(data%rlim(nlim), data%zlim(nlim)) - if (params%ifreefmt == 1) then + if (params%ifreefmt) then read(u, *) (data%rlim(i), data%zlim(i), i=1,nlim) else read(u, '(5e16.9)') (data%rlim(i), data%zlim(i), i=1,nlim) @@ -201,7 +201,7 @@ contains ! Normalize psin data%psia = psiedge - psiaxis - if(params%ipsinorm == 0) data%psin = (data%psin - psiaxis)/data%psia + if(.not. params%ipsinorm) data%psin = (data%psin - psiaxis)/data%psia end subroutine read_eqdsk diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 54ebb26..5e100b0 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -150,7 +150,7 @@ contains iroff,yynext,yypnext,yw0,ypw0, & stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) - if(params%raytracing%ipol == 0) then + if(.not. params%raytracing%ipol) then if(params%antenna%iox == 2) then ! only X mode on 1st pass cpl0 = [zero, one] else ! only O mode on 1st pass @@ -305,7 +305,7 @@ contains 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 + if(.not. params%raytracing%ipol) then ! + IF single mode propagation cpl = cpl0 p0ray(jk) = p0ray(jk)*cpl(iox) else if(cpl(iox) < etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays @@ -1842,7 +1842,8 @@ contains ! subroutine arguments real(wp_), dimension(6, nray), intent(in) :: ywrk0 real(wp_), intent(in) :: bres - integer, intent(in) :: sox, ipol + integer, intent(in) :: sox + logical, intent(in) :: ipol real(wp_), intent(inout) :: psipol0, chipol0 complex(wp_), dimension(nray), intent(out) :: ext0, eyt0 @@ -1860,7 +1861,7 @@ contains k=1 end if - if(ipol == 0) then + if(.not. ipol) then xmv=ywrk0(1:3,jk)*0.01_wp_ ! convert from cm to m anv=ywrk0(4:6,jk) rm=sqrt(xmv(1)**2+xmv(2)**2) diff --git a/src/gray_params.awk b/src/gray_params.awk new file mode 100644 index 0000000..026de74 --- /dev/null +++ b/src/gray_params.awk @@ -0,0 +1,69 @@ +#/usr/bin/env awk -f + +# This script generates the very repetitive code for parsing the +# GRAY parameters from a string. The output is embedded +# into gray_params.f90 using the #include cpp directive. + +BEGIN { + # set file to open + ARGV[1] = "src/gray_params.f90" + ARGC = 2 + + # start enum conversion switch + print "select case (value)" +} + +# match an enum constant +/enumerator :: [A-Z0-9_]+/ { + # generate conversion case + print " case ('"$3"')" + print " temp = '"$5"'" +} + +# match start of a parameters set +/type [[:alnum:]_]+_parameters$/ { + if (!done_enum) { + # close enum switch + print " case default" + print " temp = value" + print "end select\n" + # start params switch + print "select case (name)" + done_enum = 1 + } + # set the current set + split($2, words, "_parameters") + if (words[1] != "gray") set = words[1] + # change separator for the next step + FS = "::" +} + +# match a parameter +/ :: [[:alnum:]_]+/ { + if (!set) next + # strip comments + sub(/!.*$/, "") + for (i = 1; i <= split($2, ids, ","); i++) { + # extract identifier + match(ids[i], /^ *[[:alnum:]_]+/) + id = substr(ids[i], RSTART + 1, RLENGTH - 1) + # generate parameter case + print " case ('"set"."id"')" + print " read (temp, *, iostat=err) params%"set"%"id + print " if (err > 0) err = ERR_VALUE" + } +} + +# match end of a parameters set +/end type$/ { + set = "" + FS = " " +} + +END { + # close the parameters switch + print " case default" + print " ! unknown parameter" + print " err = ERR_UNKNOWN" + print "end select" +} diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 537fc44..0cdf4c3 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -7,15 +7,96 @@ module gray_params integer, parameter :: lenfnm = 256 integer, parameter :: headw = 132, headl = 21 + ! Polarisation mode + enum, bind(c) + enumerator :: polar_enum = -1 + enumerator :: MODE_O = 1 ! ordinary mode (O) + enumerator :: MODE_X = 2 ! extraordinary mode (X) + end enum + + ! Beam parameters format + enum, bind(c) + enumerator :: beam_enum = -1 + enumerator :: BEAM_0D = 0 ! fixed beam parameters + enumerator :: BEAM_1D = 1 ! 1D steering angle table + enumerator :: BEAM_2D = 2 ! 2D steering angles table + end enum + + ! Position of the X point + enum, bind(c) + enumerator :: x_point_enum = -1 + enumerator :: X_IS_MISSING = 0 ! no X point + enumerator :: X_AT_TOP = +1 ! at the top of the plasma + enumerator :: X_AT_BOTTOM = -1 ! at the bottom of the plasma + end enum + + ! MHD equilibrium kind + enum, bind(c) + enumerator :: equil_enum = -1 + enumerator :: EQ_VACUUM = 0 ! vacuum (i.e. no plasma at all) + enumerator :: EQ_ANALYTICAL = 1 ! analytical model + enumerator :: EQ_EQDSK_FULL = 2 ! G-EQDSK format - data valid on the whole domain + enumerator :: EQ_EQDSK_PARTIAL = 3 ! G-EQDSK format - data valid only inside the LCFS + end enum + + ! Temperature/density scaling model + enum, bind(c) + enumerator :: scaling_enum = -1 + enumerator :: SCALE_COLLISION = 0 ! preserve collisionality + enumerator :: SCALE_GREENWALD = 1 ! preserve Greenwald fraction + enumerator :: SCALE_OFF = 2 ! don't rescale at all + end enum + + ! Radial profile coordinate + enum, bind(c) + enumerator :: rho_enum = -1 + enumerator :: RHO_TOR = 0 ! ρ_t = √Φ (where Φ is the normalised toroidal flux) + enumerator :: RHO_POL = 1 ! ρ_p = √ψ (where ψ is the normalised poloidal flux) + enumerator :: RHO_PSI = 2 ! normalised poloidal flux ψ + end enum + + ! Plasma profiles kind + enum, bind(c) + enumerator :: profiles_enum = -1 + enumerator :: PROF_ANALYTIC = 0 ! analytical model + enumerator :: PROF_NUMERIC = 1 ! numerical data (ρ, n_e, T_e, table) + end enum + + ! Choice of the integration variable + enum, bind(c) + enumerator :: step_enum = -1 + enumerator :: STEP_ARCLEN = 0 ! arc length (s) + enumerator :: STEP_TIME = 1 ! time (actually c⋅t) + enumerator :: STEP_PHASE = 2 ! phase (actually real part of eikonal S_r=k₀⋅φ) + end enum + + ! Absorption model + enum, bind(c) + enumerator :: absorption_enum = -1 + enumerator :: ABSORP_OFF = 0 ! no absorption at all + enumerator :: ABSORP_WEAK = 1 ! weakly relativistic + enumerator :: ABSORP_FULL = 2 ! fully relativistic (faster variant) + enumerator :: ABSORP_FULL_ALT = 3 ! fully relativistic (slower variant) + end enum + + ! Current drive model + enum, bind(c) + enumerator :: current_drive_enum = -1 + enumerator :: CD_OFF = 0 ! no current drive at all + enumerator :: CD_COHEN = 1 ! Cohen + enumerator :: CD_NO_TRAP = 2 ! no trapping + enumerator :: CD_NEOCLASSIC = 3 ! Neoclassical + end enum + ! Antenna/wave launcher parameters type antenna_parameters ! From gray_params.data: - real(wp_) :: alpha, beta ! Launching angles - real(wp_) :: power ! Initial power - real(wp_) :: psi, chi ! Initial polarisation ellipse parameters - integer :: iox ! Initial wave mode - integer :: ibeam ! Beam kind - character(len=lenfnm) :: filenm ! beamdata.txt filename + real(wp_) :: alpha, beta ! Launching angles + real(wp_) :: power = 1.0_wp_ ! Initial power + real(wp_) :: psi = 0, chi = 0 ! Initial polarisation ellipse parameters + integer(kind(polar_enum)) :: iox ! Initial wave mode + integer(kind(beam_enum)) :: ibeam ! Beam parameters format + character(len=lenfnm) :: filenm ! beamdata.txt filename ! From beamdata.txt: real(wp_) :: fghz ! EC frequency @@ -27,54 +108,57 @@ module gray_params ! MHD equilibrium parameters type equilibrium_parameters - real(wp_) :: ssplps, ssplf ! Spline tensions for ψ(R,Z), I(ψ) - real(wp_) :: factb ! Rescaling factor for B - integer :: sgnb, sgni ! Sign of B, I - integer :: ixp ! Position of X point - integer :: iequil ! Equilibrium kind - integer :: icocos ! COCOS index - integer :: ipsinorm ! Normalisation of ψ - integer :: idesc, ifreefmt ! G-EQDSK quirks - character(len=lenfnm) :: filenm ! Equilibrium data filepath + real(wp_) :: ssplps = 0.005_wp_ ! Spline tension for ψ(R,Z) + real(wp_) :: ssplf = 0.01_wp_ ! Spline tension for F(ψ) + real(wp_) :: factb = 1.0_wp_ ! Rescaling factor for B + integer :: sgnb = 0, sgni = 0 ! Sign of B, I + integer(kind(x_point_enum)) :: ixp = X_IS_MISSING ! Position of X point + integer(kind(equil_enum)) :: iequil = EQ_EQDSK_FULL ! Equilibrium kind + integer :: icocos = 3 ! COCOS index + logical :: ipsinorm = .false. ! Whether ψ is normalised + logical :: idesc = .true. ! G-EQDSK quirks + logical :: ifreefmt = .false. ! + character(len=lenfnm) :: filenm ! Equilibrium data filepath end type ! Kinetic plasma profiles type profiles_parameters - real(wp_) :: psnbnd ! Plasma boundary (ψ: ne(ψ)=0) - real(wp_) :: sspld ! Density spline tension - real(wp_) :: factne, factte ! Rescale factors for n_e, T_e - integer :: iscal ! Rescaling model - integer :: irho ! Radial coordinate - integer :: iprof ! Profile kind - character(len=lenfnm) :: filenm ! Profiles data filepath + real(wp_) :: psnbnd = 1.01_wp_ ! Plasma boundary (ψ: ne(ψ)=0) + real(wp_) :: sspld = 0.1_wp_ ! Density spline tension + real(wp_) :: factne = 1.0_wp_ ! Rescale factor for n_e + real(wp_) :: factte = 1.0_wp_ ! Rescale factor for T_e + integer(kind(scaling_enum)) :: iscal = SCALE_OFF ! Rescaling model for n_e,T_e + integer(kind(rho_enum)) :: irho ! Radial coordinate + integer(kind(profiles_enum)) :: iprof ! Profile kind + character(len=lenfnm) :: filenm ! Profiles data filepath end type ! Raytracing parameters type raytracing_parameters - real(wp_) :: dst ! Integration step size - real(wp_) :: rwmax ! Normalized maximum radius of beam power - 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 - integer :: ipass ! Number of plasma passes - integer :: ipol ! Whether to compute wave polarisation + real(wp_) :: dst ! Integration step size + real(wp_) :: rwmax = 1 ! Normalized maximum radius of beam power + integer :: nray, nrayr, nrayth ! Total/radial/angular number of rays + integer :: igrad = 0 ! Complex eikonal switch + integer :: nstep = 12000 ! Max number of integration steps + integer(kind(step_enum)) :: idst = STEP_TIME ! Choice of the integration variable + integer :: ipass = 1 ! Number of plasma passes + logical :: ipol = .false. ! Whether to compute wave polarisation (from chi, psi) end type ! EC resonant heating & Current Drive parameters type ecrh_cd_parameters - integer :: iwarm ! Choice of power absorption model - integer :: ieccd ! Choice of current drive model - integer :: ilarm ! Larmor-radius expansion order - integer :: imx ! Max iterations for solving the dispersion relation + integer(kind(absorption_enum)) :: iwarm = ABSORP_FULL ! Choice of power absorption model + integer(kind(current_drive_enum)) :: ieccd = CD_NEOCLASSIC ! Choice of current drive model + integer :: ilarm = 5 ! Larmor-radius expansion order + integer :: imx = -20 ! Max iterations for solving the dispersion relation end type ! Output data parameters type output_parameters - integer :: ipec ! Grid spacing for profiles (profili EC) - integer :: nrho ! Number of grid steps for EC profiles + 1 - integer :: istpr ! Subsampling factor for beam cross-section (units 8, 12) - integer :: istpl ! Subsampling factor for outer rays (unit 33) + integer :: ipec = RHO_TOR ! Radial coordinate of the EC profiles + integer :: nrho = 501 ! Number of grid steps for EC profiles + 1 + integer :: istpr = 5 ! Subsampling factor for beam cross-section (units 8, 12) + integer :: istpl = 5 ! Subsampling factor for outer rays (unit 33) end type ! Other parameters @@ -134,7 +218,7 @@ module gray_params end type ! Global variables exposed for gray_core - integer, save :: iequil, iprof, ipol + integer, save :: iequil, iprof integer, save :: iwarm, ilarm, imx, ieccd integer, save :: igrad, idst, ipass integer, save :: istpr0, istpl0 @@ -228,7 +312,7 @@ contains write(strout(17), '("# RFL rwall:",1x,e12.5)') params%misc%rwall ! code parameters - write(strout(18), '("# COD igrad idst ipass ipol:",4(1x,i4))') & + write(strout(18), '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') & params%raytracing%igrad, params%raytracing%idst, & params%raytracing%ipass, params%raytracing%ipol write(strout(19), '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') & @@ -259,15 +343,12 @@ contains character(*), intent(in) :: name, value integer(kind(ini_error)) :: err + character(len=:), allocatable :: temp + ! hope for the best... err = ERR_SUCCESS - select case (name) #include "gray_params.inc" - case default - ! unknown parameter - err = ERR_UNKNOWN - end select end function update_parameter @@ -311,7 +392,7 @@ contains integer, intent(out) :: err ! local variables - integer :: u + integer :: u, temp(3) u = get_free_unit() @@ -324,14 +405,20 @@ contains read(u, *) params%raytracing%nrayr, params%raytracing%nrayth, & params%raytracing%rwmax - read(u, *) params%raytracing%igrad, params%raytracing%ipass, & - params%raytracing%ipol + read(u, *) params%raytracing%igrad, params%raytracing%ipass, temp(1) + + ! convert integers to logicals + params%raytracing%ipol = temp(1) > 0 + read(u, *) params%raytracing%dst, params%raytracing%nstep, & params%raytracing%idst read(u, *) params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, params%ecrh_cd%imx read(u, *) params%ecrh_cd%ieccd + ! restrict to 0-3 + params%ecrh_cd%ieccd = min(params%ecrh_cd%ieccd, 3) + read(u, *) params%antenna%alpha, params%antenna%beta read(u, *) params%antenna%power read(u, *) params%antenna%iox, params%antenna%psi, params%antenna%chi @@ -340,8 +427,13 @@ contains read(u, *) params%equilibrium%iequil read(u, *) params%equilibrium%filenm - read(u, *) params%equilibrium%icocos, params%equilibrium%ipsinorm, & - params%equilibrium%idesc, params%equilibrium%ifreefmt + read(u, *) params%equilibrium%icocos, temp(1:3) + + ! convert integers to logicals + params%equilibrium%ipsinorm = temp(1) > 0 + params%equilibrium%idesc = temp(2) > 0 + params%equilibrium%ifreefmt = temp(3) > 0 + read(u, *) params%equilibrium%ixp, params%equilibrium%ssplps, & params%equilibrium%ssplf read(u, *) params%equilibrium%sgnb, params%equilibrium%sgni, & @@ -358,6 +450,9 @@ contains read(u, *) params%output%ipec, params%output%nrho read(u, *) params%output%istpr, params%output%istpl + ! remap to 0,1 (as in irho) + params%output%ipec = mod(params%output%ipec, 2) + close(u) end subroutine read_gray_params @@ -384,7 +479,6 @@ contains mod="gray_params", proc="set_globals") end if - ipol = params%raytracing%ipol igrad = params%raytracing%igrad idst = params%raytracing%idst ipass = params%raytracing%ipass diff --git a/src/gray_params.sh b/src/gray_params.sh deleted file mode 100644 index ae21d72..0000000 --- a/src/gray_params.sh +++ /dev/null @@ -1,30 +0,0 @@ -#!/bin/sh - -# This script generates the very repetitive code for parsing the -# GRAY parameters from a string. The output is embedded -# into gray_params.f90 using a CPP include directive. - -sets='antenna equilibrium profiles raytracing ecrh_cd output misc' - -# shellcheck disable=SC2034 -{ -antenna='alpha beta power psi chi iox ibeam filenm fghz pos w ri phi' -equilibrium='ssplps ssplf factb sgnb sgni ixp iequil icocos ipsinorm idesc ifreefmt filenm' -profiles='psnbnd sspld factne factte iscal irho iprof filenm' -raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol' -ecrh_cd='iwarm ilarm imx ieccd' -output='ipec nrho istpr istpl' -misc='rwall' -} - -deref() { eval "echo \$$1"; } - -for set in $sets; do - for param in $(deref "$set"); do - cat < 0) err = ERR_VALUE -EOF - done -done diff --git a/src/multipass.f90 b/src/multipass.f90 index 46078f9..8a8e056 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -165,7 +165,8 @@ contains subroutine initmultipass(i,iox,iroff,yynext,yypnext,yw0,ypw0,stnext,p0ray, & taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) ! initialization before pass loop ! arguments - integer, intent(in) :: i, iox ! ipol, mode active on 1st pass + logical, intent(in) :: i ! ipol + integer, intent(in) :: iox ! mode active on 1st pass logical, dimension(:,:), intent(out), pointer :: iroff ! global ray status (F = active, T = inactive) real(wp_), dimension(:), intent(out), pointer :: p0ray,tau1,etau1,cpl1,lgcpl1, & psipv,chipv @@ -173,7 +174,7 @@ contains real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext ! iroff = .false. ! global ray status (F = active, T = inactive) - if(i.eq.0) call turnoffray(0,1,3-iox,iroff) ! ipol = 0 => stop other mode (iox=1/2 -> stop ib=2/1 at first pass) + if(.not. i) call turnoffray(0,1,3-iox,iroff) ! !ipol => stop other mode (iox=1/2 -> stop ib=2/1 at first pass) yynext = zero ! starting beam coordinates (1) yypnext = zero ! starting beam coordinates (2) yw0 = zero ! temporary beam coordinates (1) diff --git a/src/pec.f90 b/src/pec.f90 index cbaa5ac..91c5e1d 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -26,8 +26,8 @@ contains ! rt_in present: read input grid ! else: build equidistant grid dimension nnd +! ipec=0 rho_tor grid ! ipec=1 rho_pol grid -! ipec=2 rho_tor grid call dealloc_pec allocate(rhop_tab(nnd),rhot_tab(nnd),rtabpsi1(0:nnd),dvol(nnd),darea(nnd), & ratjav(nnd),ratjbv(nnd),ratjplv(nnd)) diff --git a/tests/06-ITER-startup/__init__.py b/tests/06-ITER-startup/__init__.py index 0cb130d..75a812d 100644 --- a/tests/06-ITER-startup/__init__.py +++ b/tests/06-ITER-startup/__init__.py @@ -21,4 +21,4 @@ class Mixed(GrayTest, TestCase): Horizontal polarisation at launch (both X, O mode) ''' reference = get_basedir(__name__) / 'outputs' / 'mixed' - gray_params = {"raytracing.ipol": 1} + gray_params = {"raytracing.ipol": True}