diff --git a/docs/examples/a10_ideal_example/slayer.in b/docs/examples/a10_ideal_example/slayer.in index 8d2d5ba0..bceddf8d 100644 --- a/docs/examples/a10_ideal_example/slayer.in +++ b/docs/examples/a10_ideal_example/slayer.in @@ -2,16 +2,16 @@ !!! FOR READING IN ASCII TABLE, SET INPUT_FLAG input_flag=f ! reads profile quantities (n_e, t_e, etc.) from an ascii table - infile='/fusion/projects/codes/gpec/GPEC-1.5/docs/examples/a10_ideal_example/a10_prof1.txt' ! Path to ascii table of profile quantities read when using input_flag - + infile='' ! Path to ascii table of profile quantities read when using input_flag + ncfile='' ! Path to netCDF file of STRIDE outputs !!! FOR MANUALLY ENTERING KINETIC QUANTITIES @ RATIONAL SURFACE, SET PARAMS_FLAG - params_flag=t ! calculates normalized parameters from profile quantities + params_flag=f ! calculates normalized parameters from profile quantities mm=2 ! poloidal mode number nn=1 ! toroidal mode number n_e=7.77E+17 ! electron density [m^-3] t_e=25.8 ! electron temperature [eV] t_i=25.8 ! ion temperature [eV] - sval=729.0 ! magnetic shear at the layer + sval=2.0 ! magnetic shear at the layer bt=1.0 ! toroidal field [T] rs=0.17 ! minor radius of resonant surface [m] R0=2.0 ! major radius of magnetic axis [m] @@ -21,7 +21,7 @@ qval=2.0 ! q of resonant surface mu_i=2.0 ! ion mass ratio to proton ?? zeff=2.0 ! plasma Z_effective - inpr=0.0001 ! Prantdl number + inpr=0.1 ! Prantdl number inpe=0.0 ! Electron viscosity !!! IF INPUT_FLAG & PARAMS_FLAG ARE FALSE, MANUALLY ENTER DIMENSIONLESS QUANTITIES @ RATIONAL SURFACE @@ -56,7 +56,8 @@ bin_flag=t ! writes results to binary files netcdf_flag=f ! writes results to netcdf files stability_flag=f ! calculate delta dependence on complex Q - bal_flag=t ! calculate the resonant field penetration threshold from torque balance + est_gamma_flag=t ! Calculate estimated growth rates on each rational surface + !bal_flag=f ! calculate the resonant field penetration threshold from torque balance / &SLAYER_DIAGNOSE riccati_out=f ! writes LSDOE Riccati integration to an ascii file diff --git a/input/slayer.in b/input/slayer.in index 65d4d452..24077b95 100644 --- a/input/slayer.in +++ b/input/slayer.in @@ -1,45 +1,60 @@ &SLAYER_INPUT - - !!! FOR READING IN ASCII TABLE, SET INPUT_FLAG = TRUE input_flag=f ! reads profile quantities (n_e, t_e, etc.) from an ascii table - infile='' ! Path to ascii table of profile quantities read when using input_flag - ncfile='' ! Path to netCDF file of STRIDE outputs - - !!! FOR MANUALLY ENTERING KINETIC QUANTITIES @ RATIONAL SURFACE, SET PARAMS_FLAG = TRUE + infile='' ! Path to ascii table of profile quantities read when using input_flag + ncfile='' ! Path to netCDF file of STRIDE outputs params_flag=f ! calculates normalized parameters from profile quantities + ! >>> + ! >>> Input rational surface quantities + ! >>> mm=2 ! poloidal mode number nn=1 ! toroidal mode number - n_e=7.77E+17 ! electron density [m^-3] - t_e=25.8 ! electron temperature [eV] - t_i=25.8 ! ion temperature [eV] - sval=2.0 ! magnetic shear at the layer - bt=1.0 ! toroidal field [T] - rs=0.17 ! minor radius of resonant surface [m] - R0=2.0 ! major radius of magnetic axis [m] - omega=5.45E+04 ! ExB frequency [1/s] + n_e=5.50E+19 ! electron density [m^-3] + t_e=1000.0 ! electron temperature [eV] + t_i=1000.0 ! ion temperature [eV] + sval=1.5 ! magnetic shear at the layer + bt=2.0 ! toroidal field [T] + rs=0.5 ! minor radius of resonant surface [m] + R0=1.67 ! major radius of magnetic axis [m] + omega=5.00E+03 ! ExB frequency [Hz] l_t=0.11 ! temperature gradient scale length l_n=0.25 ! density gradient scale length qval=2.0 ! q of resonant surface - mu_i=2.0 ! ion mass ratio to proton ?? + mu_i=2.0 ! ion mass ratio to proton mass, 2 for deuterium plasma zeff=2.0 ! plasma Z_effective - inpr=10.0 ! Prantdl number - inpr_prof=10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0 ! Prandtl number profile - inpe=0.0 ! Electron viscosity (if > 0.0, Pe_flag=t will be overridden and this value used instead) - - !!! IF INPUT_FLAG & PARAMS_FLAG ARE FALSE, MANUALLY ENTER DIMENSIONLESS QUANTITIES @ RATIONAL SURFACE HERE - inQ=23.0 ! normalized ExB frequency - inQ_e=2.0 ! normalized electron diamagnetic frequency - inQ_i=-2.0 ! normalized ion diamagnetic frequency - inc_beta=0.7 ! dimensionless measure of the plasma pressure - inds=6.0 ! normalized ion sound radius - intau=1.0 ! ion temperature divided by electron temperature - Q0=4.0 ! unnecessary placeholder variable for inQ? - delta_n_p=0 ! delta offset used for jxb calculations. Default is (1e-2,1e-2). + dr_val=-0.1 ! GGJ resistive interchange criterion D_R + dgeo_val=10.0 ! J.W. Connor geometric prefactor for Delta_crit + chi_p_prof=0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2 ! Rational surface anomalous perpendicular energy diffusivity profile + chi_t_prof=0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2 ! Rational surface anomalous perpendicular ion momentum diffusivity profile + kappa_prof=0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 ! Rational surface thermal conductivity profile (NOTE: set > 0 to override chi_p_prof) + ! >>> + ! >>> Input normalized layer parameters. If read_eq=f, each will individually override internal calcuation using above inputs + ! >>> + inpr=0.0 ! Prantdl number. Set /= 0.0 to override internal calculation + inpe=0.0 ! electron viscosity + inQ=1.0 ! normalized ExB frequency. TO DO: update for mode rotation frequency calculation! + inQ_e=0.0 ! normalized electron diamagnetic frequency. Set /= 0.0 to override internal calculation + inQ_i=0.0 ! normalized ion diamagnetic frequency. Set /= 0.0 to override internal calculation + inc_beta=0.0 ! dimensionless measure of the plasma pressure. Set /= 0.0 to override internal calculation + inds=0.0 ! normalized ion sound radius, or normalized ion skin depth. Set /= 0.0 to override internal calculation + intau=1.0 ! ratio of ion to electron temperature + Q0=1.0 ! unnecessary placeholder variable for inQ? + delta_prime=(0.0,0.0) ! input Delta' or (Delta' - Delta_crit) for growth rate calculation, set /= 0.0 to override internal calculation + delta_n_p=(0.01,0.01) ! delta offset used for jxb calculations. Default for latter is (1e-2,1e-2). + ingamma=(0.0,0.0) ! input (Re(Q),Im(Q)) to calculate inner layer Delta. Set /= 0.0 to override internal root finding / + &SLAYER_CONTROL inum=400 ! resolution to find error field thresholds. jnum=100 ! resolution for 2d scan along with Q,omega. - knum=100 ! resolution for 2d scan alont with the other. + knum=100 ! resolution for 2d scan along with the other. + Q_num=100 ! resolution for stability scan along Re(Q),Im(Q) axes + scan_width=1.5 ! stability scan width for Re(Q),Im(Q) scan + AMR_passes=4 + msing_max=2 ! number of surfaces to include in coupled tearing growth rate calculation + dc_type="toroidal" ! Delta_crit type, options are "toroidal", "lar", and "rfitzp" + read_eq=f ! read in equilibrium quantities from STRIDE and kinetic file. Set all normalized params to 0.0 + Pperp_Ptor_flag=t ! Use R. Fitzpatrick formalism (inc. anomalous diffusion) for growth rate calculation + coupling_flag=t ! Solve determinant problem and output 2D scan QPscan_flag=f ! scan (Q,P) space for delta and torque. Qscan_flag=f ! scan Q space QPescan_flag=f ! scan (Q,Pe) space for delta and torque. @@ -50,18 +65,27 @@ nbtscan_flag=f ! scan (n,bt) space for error fields. parflow_flag=f ! set parallel flow on PeOhmOnly_flag=t ! only include Pe from Ohm's law - Pe_flag=f ! Not operational yet. If true, it will calculate classical inpe value from inpr and include in SLAYER run + Pe_flag=f ! if true, calculate classical inpe value from inpr and include in SLAYER run layfac=0.02 ! layfac*EXP(ifac*ATAN2(AIMAG(Q-Q_e),REAL(Q-Q_e)) is added to Q_e if ABS(Q-Q_e) $(BINDIR)/rmatch) # SLAYER +../slayer/version.inc: + @ver="`git describe --tags`"; \ + if [ ! -f $@ ] || ! grep -F "'$$ver'" $@ >/dev/null 2>&1; then \ + echo ">>> Updating version file"; \ + echo " CHARACTER(len=*), PARAMETER :: version ='$$ver'" | tee $@; \ + else \ + echo ">>> Git version unchanged!"; \ + fi + $(call print_success,version file updated: $@) + SLAYER_LIBS = \ + -lpentrc \ -llsode \ -lequil \ + -lharvest \ $(MATHLIBS) \ $(NETCDFLIBS) \ $(NETCDF_EXTRA_LIBS) @@ -317,7 +329,7 @@ SLAYER_LIBDIRS = \ -L$(LIBDIR) \ -L$(MATHDIR) \ -L$(NETCDFDIR) -slayer: v neededdeps $(SLAYER_OBJS) lsode equil +slayer: v neededdeps $(SLAYER_OBJS) lsode equil pentrc harvest mkdir -p $(BINDIR) rm -f ../slayer/slayer $(FCSLAYER) -o ../slayer/slayer $(SLAYER_OBJS) $(SLAYER_LIBDIRS) $(SLAYER_LIBS) $(LDFLAGS) diff --git a/pentrc/inputs.f90 b/pentrc/inputs.f90 index bff8aeec..e50d2e9d 100644 --- a/pentrc/inputs.f90 +++ b/pentrc/inputs.f90 @@ -217,6 +217,7 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log) do i=1,5 tmp%fs(0:,i) = table(1:,i+1) enddo + call spline_fit(tmp,"extrap") if(write_log) print *,"Formed temporary spline" diff --git a/rmatch/match.f b/rmatch/match.f index be4fa63a..8b989064 100644 --- a/rmatch/match.f +++ b/rmatch/match.f @@ -92,12 +92,16 @@ MODULE match_mod $ deltac_flag=.FALSE.,deltaj_flag=.FALSE., $ match_flag=.FALSE. LOGICAL :: bin_rpecsol=.FALSE.,out_rpecsol=.FALSE. + LOGICAL :: detgrid_flag=.FALSE. CHARACTER(10) :: model="deltac" INTEGER :: msing,totmsing,nstep=32,qscan_ising=1 INTEGER :: scan_nstep, scan_estep INTEGER :: nroot=1,iroot,totnsol,ising_output=1,itermax=500 + INTEGER :: detgrid_nre=64, detgrid_nim=64 REAL(r8) :: eta(20),dlim=1000,massden(20),rotation(20)=0,ntor=1 REAL(r8) :: scan_x0,scan_x1,relax_fac,scan_e0,scan_e1 + REAL(r8) :: detgrid_re_min=-1.0, detgrid_re_max=1.0, + $ detgrid_im_min=-1.0, detgrid_im_max=1.0 REAL(r8), DIMENSION(:), ALLOCATABLE :: taur_save REAL(r8), DIMENSION(:), ALLOCATABLE :: zo_out,zi_in COMPLEX(r8) :: initguess @@ -137,7 +141,10 @@ SUBROUTINE match_run $ deflate,nroot,match_flag,ising_output, $ match_sol,matrix_diagnose,fulldomain, $ coil,itermax,relax_fac,init_scan_flag, - $ scan_e0,scan_e1,eqscan_flag,scan_estep + $ scan_e0,scan_e1,eqscan_flag,scan_estep, + $ detgrid_flag,detgrid_nre,detgrid_nim, + $ detgrid_re_min,detgrid_re_max, + $ detgrid_im_min,detgrid_im_max NAMELIST/rmatch_output/ bin_rpecsol,out_rpecsol NAMELIST/nyquist_input/nyquist 10 FORMAT(1x,"Eigenvalue=",1p,2e11.3) @@ -246,7 +253,12 @@ SUBROUTINE match_run c----------------------------------------------------------------------- c scan eigen value (Q) for different inner models. c----------------------------------------------------------------------- - IF(qscan_flag) CALL match_qscan + IF(qscan_flag) CALL match_qscan +c----------------------------------------------------------------------- +c 2D complex-Q grid scan of the full match_delta determinant +c (patch added for Julia↔Fortran apples-to-apples comparison). +c----------------------------------------------------------------------- + IF(detgrid_flag) CALL match_detgrid c----------------------------------------------------------------------- c nyquist plot. c----------------------------------------------------------------------- @@ -475,6 +487,17 @@ FUNCTION match_delta(guess,mat) RESULT(det) zi_in(ising)=0 q_in(ising)=q_deltac sol=0 +c -- PATCH (Julia/Fortran benchmark) -- +c Dump per-surface Δ for each match_delta call so we can compare +c Julia's Galerkin output to Fortran's at the same Q. Writes to +c delta_per_surface.out (appended each call; nuke before a run). + OPEN(UNIT=42, FILE='delta_per_surface.out', + $ POSITION='APPEND') + WRITE(42,'(2e22.14, i5, 4e22.14)') + $ REAL(guess_modify), AIMAG(guess_modify), ising, + $ REAL(deltar(ising,1)), AIMAG(deltar(ising,1)), + $ REAL(deltar(ising,2)), AIMAG(deltar(ising,2)) + CLOSE(UNIT=42) END SELECT c----------------------------------------------------------------------- c construct the matching matrix. @@ -963,7 +986,52 @@ SUBROUTINE match_qscan c----------------------------------------------------------------------- CALL program_stop("Normal termination for q scan.") END SUBROUTINE match_qscan - + +c----------------------------------------------------------------------- +c subprogram 7b. match_detgrid. +c scan match_delta on a 2D Q_re × Q_im grid and write an ascii dump +c (detgrid.out, format "qre qim re_det im_det") plus a binary +c (detgrid.bin). Patch added for Julia↔Fortran apples-to-apples +c comparison of the coupled GGJ dispersion relation. +c----------------------------------------------------------------------- + SUBROUTINE match_detgrid + INTEGER :: ire, iim + REAL(r8) :: qre, qim, dre, dim + COMPLEX(r8) :: guess, det + COMPLEX(r8), DIMENSION(4*msing,4*msing) :: mat + 10 FORMAT(1p,4e20.11) + WRITE(*,*) "DETGRID: 2D scan of match_delta(Q, mat)" + WRITE(*,'(2x,a,i4,a,i4)') "nre=",detgrid_nre," nim=",detgrid_nim + WRITE(*,'(2x,a,2e12.3)') "Re(Q) range:", + $ detgrid_re_min, detgrid_re_max + WRITE(*,'(2x,a,2e12.3)') "Im(Q) range:", + $ detgrid_im_min, detgrid_im_max + OPEN(UNIT=bin_unit,FILE="detgrid.bin",STATUS="REPLACE", + $ FORM="UNFORMATTED") + CALL ascii_open(out_unit,"detgrid.out","REPLACE") + WRITE(out_unit,'(a)') "# Q_re Q_im Re(det) Im(det)" + WRITE(bin_unit) detgrid_nre, detgrid_nim, msing + dre = (detgrid_re_max - detgrid_re_min) / MAX(1, detgrid_nre - 1) + dim = (detgrid_im_max - detgrid_im_min) / MAX(1, detgrid_nim - 1) + DO iim = 1, detgrid_nim + qim = detgrid_im_min + (iim-1) * dim + DO ire = 1, detgrid_nre + qre = detgrid_re_min + (ire-1) * dre + guess = CMPLX(qre, qim, r8) + det = match_delta(guess, mat) + WRITE(out_unit,10) qre, qim, REAL(det,r8), AIMAG(det) + WRITE(bin_unit) REAL(qre,4), REAL(qim,4), + $ REAL(REAL(det,r8),4), REAL(AIMAG(det),4) + ENDDO + IF (MOD(iim, MAX(1,detgrid_nim/10)) == 0) + $ WRITE(*,'(4x,a,i4,a,i4)') "row ",iim," /",detgrid_nim + ENDDO + CALL ascii_close(out_unit) + WRITE(bin_unit) + CLOSE(UNIT=bin_unit) + WRITE(*,*) "DETGRID: wrote detgrid.out and detgrid.bin" + CALL program_stop("Normal termination for detgrid scan.") + END SUBROUTINE match_detgrid c----------------------------------------------------------------------- c subprogram 8. match_delta_jardin. c finite differential method of GGJ. diff --git a/slayer/delta.f b/slayer/delta.f index 44f3d7c4..8ef9bff2 100644 --- a/slayer/delta.f +++ b/slayer/delta.f @@ -1,32 +1,81 @@ MODULE delta_mod +c----------------------------------------------------------------------- +c delta_mod: Riccati-based tearing-mode layer Delta solvers. +c +c Contains three Riccati formulations: +c riccati - standard formulation (non-stiff, lsode mf=10) +c riccati_del_s - del_s formulation (stiff, lsode mf=21) +c riccati_f - Fitzpatrick P_perp/P_tor formulation (mf=21) +c +c Each solver integrates a Riccati ODE for W(x) from a large-|x| +c asymptotic boundary condition inward to x ~ 0, then extracts +c Delta = pi / W'(0). +c +c Associated ODE subroutines (w_der, w_der_del_s, w_der_f) and +c Jacobian subroutines (jac_del_s, jac_f) follow below. +c +c Module-level flags: +c riccati_out - write W(x) profile to binary + text file +c parflow_flag - include parallel electron flow terms in w_der +c PeOhmOnly_flag - retain only Pe-Ohm coupling in w_der +c----------------------------------------------------------------------- USE sglobal_mod IMPLICIT NONE - LOGICAL :: riccati_out,parflow_flag,PeOhmOnly_flag +c --- module-level control flags + LOGICAL :: riccati_out ! write W(x) profile to binary + text + LOGICAL :: parflow_flag ! enable parallel-flow terms in w_der + LOGICAL :: PeOhmOnly_flag ! Pe-Ohm-only coupling in w_der CONTAINS c----------------------------------------------------------------------- -c calculate delta based on riccati w_der formulation. +c riccati: compute Delta via Riccati integration of w_der. +c +c Integrates W(x) from x_start inward to x_min using lsode +c (non-stiff, mf=10). Optionally applies the layfac singularity +c guard when Q is near Q_e. Returns Delta = pi / W'(x_min). +c +c If riccati_out = .TRUE., writes the W(x) profile to +c slayer_riccati_profile_n.{bin,out}. +c c----------------------------------------------------------------------- FUNCTION riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe, $ iinQ,inx,iny) - REAL(r8),INTENT(IN) :: inQ,inQ_e,inQ_i,inpr,inpe,inc_beta,inds - REAL(r8),INTENT(IN) :: intau - REAL(r8),INTENT(IN),OPTIONAL :: iinQ,inx - COMPLEX(r8), INTENT(IN), OPTIONAL :: iny +c --- input arguments: physical parameters for this surface + REAL(r8),INTENT(IN) :: inQ ! real part of Q + REAL(r8),INTENT(IN) :: inQ_e ! electron diamagnetic freq + REAL(r8),INTENT(IN) :: inQ_i ! ion diamagnetic freq + REAL(r8),INTENT(IN) :: inpr ! Prandtl number + REAL(r8),INTENT(IN) :: inpe ! electron Prandtl number + REAL(r8),INTENT(IN) :: inc_beta ! c_beta parameter + REAL(r8),INTENT(IN) :: inds ! magnetic shear ds + REAL(r8),INTENT(IN) :: intau ! ion-to-electron temp ratio +c --- optional arguments + REAL(r8),INTENT(IN),OPTIONAL :: iinQ ! imaginary part of Q + REAL(r8),INTENT(IN),OPTIONAL :: inx ! override starting x + COMPLEX(r8),INTENT(IN),OPTIONAL :: iny ! override starting W +c --- function result COMPLEX(r8) :: riccati - INTEGER :: istep,neq,itol,itask,istate,liw,lrw,iopt,mf - - REAL(r8) :: xintv,x,xout,rtol,jac,xmin - COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: y,dy +c --- lsode solver control + INTEGER :: istep ! integration step counter + INTEGER :: neq ! number of equations (=2) + INTEGER :: itol,itask ! lsode tolerance/task flags + INTEGER :: istate,iopt,mf ! lsode state/option/method flags + INTEGER :: liw,lrw ! lsode work array sizes + REAL(r8) :: x,xout ! current and target x + REAL(r8) :: xmin ! inner integration bound + REAL(r8) :: rtol ! relative tolerance + REAL(r8) :: jac ! dummy Jacobian (unused for mf=10) +c --- work arrays + COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: y,dy ! W and dW + INTEGER, DIMENSION(:), ALLOCATABLE :: iwork ! lsode int work + REAL(r8), DIMENSION(:), ALLOCATABLE :: atol,rwork ! lsode real work - INTEGER, DIMENSION(:), ALLOCATABLE :: iwork - REAL(r8), DIMENSION(:), ALLOCATABLE :: xfac,atol,rwork - +c --- copy input arguments to module-level globals for w_der Q=inQ IF(present(iinQ)) Q=inQ+ifac*iinQ Q_e=inQ_e @@ -36,76 +85,415 @@ FUNCTION riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe, c_beta=inc_beta ds=inds tau=intau - + +c --- singularity guard: displace Q away from Q_e when too close IF ((layfac>0).AND.(ABS(Q-Q_e)xout) istep=istep+1 CALL lsode(w_der,neq,y,x,xout,itol,rtol,atol, $ itask,istate,iopt,rwork,lrw,iwork,liw,jac,mf) - WRITE(bin_unit)REAL(x,4),REAL(REAL(y),4),REAL(AIMAG(y),4) + WRITE(bin_unit)REAL(x,4),REAL(REAL(y),4),REAL(AIMAG(y),4) WRITE(out2_unit,'(1x,3(es17.8e3))') x,REAL(y),AIMAG(y) - ENDDO + ENDDO CLOSE(bin_unit) CLOSE(out2_unit) ELSE +c single-shot integration to x_min istep = 1 itask = 1 CALL lsode(w_der,neq,y,x,xout,itol,rtol,atol, $ itask,istate,iopt,rwork,lrw,iwork,liw,jac,mf) - ENDIF - ! w=0 when Q=Q_e. Why? - +c --- extract Delta from final W derivative at x_min +c NOTE: W -> 0 when Q -> Q_e (see layfac guard above). CALL w_der(neq,x,y,dy) riccati=pi/dy(1) - DEALLOCATE(atol,y,dy,iwork,rwork) + DEALLOCATE(atol,y,dy,iwork,rwork) END FUNCTION riccati c----------------------------------------------------------------------- -c riccati integration. +c riccati_del_s: compute Delta via the del_s Riccati formulation. +c +c Uses a stiff solver (lsode mf=21) with user-supplied Jacobian +c (jac_del_s). Integrates W(q) from large q inward to q_min. +c Returns Delta = -(pi / sqrt(1+1/tau)) * W'(q_min). +c +c Module-level variables tau, D_norm, and P_tor must be valid +c before entry; callers are responsible for setting them. +c Q_e, Q_i, and P_perp are set from arguments inside this +c function. Q, c_beta, and d_beta are not referenced here. +c Note: slayer.f passes Q_e_arr (not Q_arr) as the first +c argument -- this is intentional. +c----------------------------------------------------------------------- + FUNCTION riccati_del_s(inQ_e,inQ_i,inpr,inx,iny) + +c --- input arguments + REAL(r8),INTENT(IN) :: inQ_e ! electron diamagnetic freq + REAL(r8),INTENT(IN) :: inQ_i ! ion diamagnetic freq + REAL(r8),INTENT(IN) :: inpr ! mapped to P_perp (see below) +c --- optional arguments + REAL(r8),INTENT(IN) :: inx ! starting q for integration + COMPLEX(r8),INTENT(IN),OPTIONAL :: iny ! override starting W +c --- function result + COMPLEX(r8) :: riccati_del_s + +c --- lsode solver control + INTEGER :: istep ! integration step counter + INTEGER :: neq ! number of equations (=2) + INTEGER :: itol,itask ! lsode tolerance/task flags + INTEGER :: istate,iopt,mf ! lsode state/option/method flags + INTEGER :: liw,lrw ! lsode work array sizes + REAL(r8) :: x ! secondary x variable (set from inx) + REAL(r8) :: xout ! target integration endpoint + REAL(r8) :: xmin ! inner integration bound + REAL(r8) :: rtol ! relative tolerance + REAL(r8) :: jac ! dummy (overridden by jac_del_s) + REAL(r8) :: my_q ! integration variable (large -> small) + REAL(r8) :: P_hat ! normalized P_perp + REAL(r8) :: alpha ! boundary condition coefficient +c --- work arrays + COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: W,dW_dq ! W and dW/dq + INTEGER, DIMENSION(:), ALLOCATABLE :: iwork ! lsode int work + REAL(r8), DIMENSION(:), ALLOCATABLE :: atol,rwork ! lsode real work + +c --- configure lsode: stiff BDF method with user Jacobian (mf=21) + neq = 2 + itol = 2 + rtol = 1e-10 + ALLOCATE(atol(neq),W(1),dW_dq(1)) + atol(:) = 1e-10 + itask = 2 + istate = 1 + iopt = 1 ! enable optional inputs (iwork(6)) + mf = 21 ! stiff, user-supplied Jacobian (jac_del_s) + liw = 20*2 + lrw = 22+9*neq+neq**2 ! stiff work array size + ALLOCATE(iwork(liw+neq),rwork(lrw)) + +c --- set maximum internal steps + iwork=0 + iwork(6)=50000 ! MXSTEP: max internal steps per call + rwork=0 + +c --- set starting integration point + my_q=inx ! start backwards integration at large q + xmin=1e-5 + x=inx + xout=xmin + +c --- copy input arguments to module-level globals for w_der_del_s + Q_e = inQ_e + Q_i = inQ_i + P_perp = inpr + + P_hat = P_perp / D_norm**6.0 + +c --- asymptotic boundary condition at large q + alpha = (P_hat/(1+1/tau))**0.5 + W(1) = -alpha*my_q**2 - 0.5 + IF(present(iny)) W(1)=iny + +c --- integrate W(q) from q_start inward to q_min via lsode + IF (riccati_out) THEN +c profile output: step-by-step integration with file writes + istep = 1 + itask = 2 + OPEN(UNIT=bin_unit,FILE='slayer_riccati_profile_n'// + $ TRIM(sn_str)//'.bin',STATUS='UNKNOWN', + $ POSITION='REWIND',FORM='UNFORMATTED') + + OPEN(UNIT=out2_unit,FILE='slayer_riccati_profile_n'// + $ TRIM(sn_str)//'.out',STATUS='UNKNOWN') + WRITE(out2_unit,'(1x,3(a17))'),"x","RE(y)","IM(y)" + DO WHILE (my_q>xout) + istep=istep+1 + CALL lsode(w_der_del_s,neq,W,my_q,xout,itol,rtol,atol, + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_del_s,mf) + WRITE(bin_unit)REAL(my_q,4),REAL(REAL(W),4),REAL(AIMAG(W),4) + WRITE(out2_unit,'(1x,3(es17.8e3))')my_q,REAL(W),AIMAG(W) + ENDDO + CLOSE(bin_unit) + CLOSE(out2_unit) + ELSE +c single-shot integration to q_min + istep = 1 + itask = 1 + CALL lsode(w_der_del_s,neq,W,my_q,xout,itol,rtol,atol, + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_del_s,mf) + ENDIF + +c --- extract Delta from final W derivative at q_min + CALL w_der_del_s(neq,my_q,W,dW_dq) + riccati_del_s=-( pi/((1+1/tau)**0.5) )*dW_dq(1) + DEALLOCATE(atol,W,dW_dq,iwork,rwork) + + END FUNCTION riccati_del_s +c----------------------------------------------------------------------- +c jacobian for riccati_del_s(): pd = dF/dW for stiff lsode. +c----------------------------------------------------------------------- + SUBROUTINE jac_del_s(neq, my_q, W, ml, mu, pd, nrpd) + INTEGER, INTENT(IN) :: neq, ml, mu, nrpd + REAL(r8), INTENT(IN) :: my_q + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: W + COMPLEX(r8), DIMENSION(nrpd,neq), INTENT(INOUT) :: pd + pd(1,1) = 1.0/my_q - 2.0d0*W(1)/my_q + END SUBROUTINE jac_del_s +c----------------------------------------------------------------------- +c w_der_del_s: ODE right-hand side dW/dq for riccati_del_s. +c Implements the del_s dispersion relation using normalised +c quantities Q_hat, P_perp_hat, P_tor_hat. +c----------------------------------------------------------------------- + SUBROUTINE w_der_del_s(neq,my_q,W,dW_dq) + + INTEGER, INTENT(IN) :: neq + REAL(r8), INTENT(IN) :: my_q + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: W + COMPLEX(r8), DIMENSION(neq), INTENT(OUT) :: dW_dq + REAL(r8) :: Q_hat, P_tor_hat, P_perp_hat + COMPLEX(r8) :: E,F + +c --- normalise physical quantities + Q_hat = (Q_e*(1+tau)/tau) / D_norm**4.0 + P_perp_hat = P_perp / D_norm**6.0 + P_tor_hat = P_tor / D_norm**6.0 +c --- build the E and F dispersion coefficients + E = (-(Q_hat**2)/(1+1/tau)) - ifac*Q_hat*(P_perp_hat+ + $ P_tor_hat)*(my_q**2) + P_perp_hat*P_tor_hat*(my_q**4) + F = P_perp_hat - ifac*Q_hat + (1+1/tau)*P_tor_hat*my_q**2 + +c --- Riccati ODE: dW/dq = W/q - W^2/q + q*E/F + dW_dq(1)=W(1)/my_q - (W(1)**2)/my_q + (my_q*E)/F + RETURN + END SUBROUTINE w_der_del_s +c----------------------------------------------------------------------- +c riccati_f: compute Delta via Fitzpatrick P_perp / P_tor +c Riccati formulation. +c +c Uses a stiff solver (lsode mf=21) with user-supplied Jacobian +c (jac_f). Boundary conditions are set analytically from the +c large-p asymptotic behaviour; the branch depends on whether +c D_norm^2 exceeds iota_e * P_perp / P_tor^(2/3). +c +c----------------------------------------------------------------------- + FUNCTION riccati_f() + +c --- function result + COMPLEX(r8) :: riccati_f + +c --- lsode solver control + INTEGER :: istep ! integration step counter + INTEGER, PARAMETER :: neq = 2 ! number of equations + INTEGER :: itol,itask ! lsode tolerance/task flags + INTEGER :: istate,iopt,mf ! lsode state/option/method flags + INTEGER, PARAMETER :: liw = 42 ! 20*2 + neq + INTEGER, PARAMETER :: lrw = 44 ! 22 + 9*neq + neq**2 + REAL(r8) :: xout ! target integration endpoint + REAL(r8) :: xmin ! inner integration bound + REAL(r8) :: rtol ! relative tolerance + REAL(r8) :: jac ! dummy (overridden by jac_f) + REAL(r8) :: my_p ! integration variable p (large -> small) + REAL(r8) :: bk ! asymptotic coefficient b_k +c --- boundary-condition intermediates + COMPLEX(r8) :: ak ! asymptotic coefficient a_k + COMPLEX(r8) :: ck ! asymptotic coefficient c_k + COMPLEX(r8) :: xk ! asymptotic coefficient x_k + COMPLEX(r8) :: W_bound ! boundary value for W(p_start) +c --- work arrays (fixed-size, stack-allocated to avoid heap overhead) + COMPLEX(r8) :: W(1), dWdp(1) ! W and dW/dp + REAL(r8) :: atol(neq) ! absolute tolerance + INTEGER :: iwork(liw) ! lsode int work + REAL(r8) :: rwork(lrw) ! lsode real work + +c --- configure lsode: stiff BDF method with user Jacobian (mf=21) + itol = 2 + rtol = 1.0d-10 + atol(:) = 1.0d-10 + itask = 2 + istate = 1 + iopt = 1 ! enable optional inputs (iwork(6)) + mf = 21 ! stiff, user-supplied Jacobian (jac_f) + +c --- set maximum internal steps (only iwork(5:6) matter for lsode) + iwork(:) = 0 + iwork(6) = 50000 ! MXSTEP: max internal steps per call + rwork(:) = 0.0d0 + + xmin=1e-6 + xout=xmin + +c --- compute starting p and W boundary condition +c Branch on asymptotic regime: D_norm^2 vs iota_e*P_perp/P_tor^(2/3) +c --- branch 1: D_norm^2 > iota_e * P_perp / P_tor^(2/3) +c large-D_norm regime: p scales with (P_tor*D_norm^2/(iota_e*...)) + IF (D_norm**2 > iota_e*P_perp / P_tor**(2.0d0/3.0d0)) THEN + my_p = ( (P_tor*D_norm**2)/(iota_e*P_tor*P_perp) )**0.25d0 + my_p = MAX(my_p, 6.0d0) + + ak = -(g_tmp + ifac*Q_e) + bk = (iota_e*P_perp*P_tor)/(P_tor*D_norm**2) + + ck = bk*(1 + (g_tmp+ifac*Q_i) + $ *((P_tor+P_perp)/(P_tor*P_perp)) + $ - (P_perp + (g_tmp+ifac*Q_i)*D_norm**2) + $ *(iota_e/(P_tor*D_norm**2))) + + xk = (ck - SQRT(bk)*(1 - SQRT(bk)*ak)) / (2.0d0*SQRT(bk)) + + W_bound = xk - SQRT(bk)*my_p + ELSE +c --- branch 2: D_norm^2 <= iota_e * P_perp / P_tor^(2/3) +c small-D_norm regime: p scales with 1/P_tor^(1/6) + my_p = 1.0d0 / P_tor**(1.0d0/6.0d0) + my_p = MAX(my_p, 6.0d0) + + ak = -(g_tmp + ifac*Q_e) + bk = P_tor + ck = -ifac*(Q_e - Q_i)*(P_tor/P_perp) + (g_tmp + ifac*Q_i) + xk = (ak*bk - ck)/(2.0d0*SQRT(bk)) + + W_bound = -1.0d0 + xk*my_p - SQRT(bk)*my_p**3 + END IF + + W(1) = W_bound + +c --- integrate W(p) from p_start inward to p_min via lsode + IF (riccati_out) THEN +c profile output: step-by-step integration with file writes + istep = 1 + itask = 2 + OPEN(UNIT=bin_unit,FILE='slayer_riccati_profile_n'// + $ TRIM(sn_str)//'.bin',STATUS='UNKNOWN', + $ POSITION='REWIND',FORM='UNFORMATTED') + + OPEN(UNIT=out2_unit,FILE='slayer_riccati_profile_n'// + $ TRIM(sn_str)//'.out',STATUS='UNKNOWN') + WRITE(out2_unit,'(1x,3(a17))'),"x","RE(y)","IM(y)" + DO WHILE (my_p>xout) + istep=istep+1 + CALL lsode(w_der_f,neq,W,my_p,xout,itol,rtol,atol, + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_f,mf) + WRITE(bin_unit)REAL(my_p,4),REAL(REAL(W),4), + $ REAL(AIMAG(W),4) + WRITE(out2_unit,'(1x,3(es17.8e3))')my_p,REAL(W),AIMAG(W) + ENDDO + CLOSE(bin_unit) + CLOSE(out2_unit) + ELSE +c single-shot integration to p_min + istep = 1 + itask = 1 + CALL lsode(w_der_f,neq,W,my_p,xout,itol,rtol,atol, + $ itask,istate,iopt,rwork,lrw,iwork,liw,jac_f,mf) + ENDIF + +c --- extract Delta from final W derivative at p_min + CALL w_der_f(neq,my_p,W,dWdp) + riccati_f = pi / dWdp(1) + + END FUNCTION riccati_f +c----------------------------------------------------------------------- +c jacobian for riccati_f(): pd = dF/dW for stiff lsode. +c----------------------------------------------------------------------- + SUBROUTINE jac_f(neq, my_p, W, ml, mu, pd, nrpd) + INTEGER, INTENT(IN) :: neq, ml, mu, nrpd + REAL(r8), INTENT(IN) :: my_p + COMPLEX(r8) :: fA_p, denom + REAL(r8) :: p2 + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: W + COMPLEX(r8), DIMENSION(nrpd,neq), INTENT(INOUT) :: pd + + p2 = my_p * my_p + denom = g_tmp + ifac*Q_e + p2 + fA_p = (denom - 2.0d0*p2) / denom + + pd(1,1) = (-fA_p/my_p) - (2.0d0*W(1))/my_p + END SUBROUTINE jac_f +c----------------------------------------------------------------------- +c w_der_f: ODE right-hand side dW/dp for riccati_f. +c Implements the Fitzpatrick P_perp / P_tor dispersion relation. +c Coefficients fA, fB, fC are evaluated at the current p. +c----------------------------------------------------------------------- + SUBROUTINE w_der_f(neq,my_p,W,dWdp) + + INTEGER, INTENT(IN) :: neq + REAL(r8), INTENT(IN) :: my_p + COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: W + COMPLEX(r8), DIMENSION(neq), INTENT(OUT) :: dWdp + COMPLEX(r8) :: fA, fA_prime, fB, fC + COMPLEX(r8) :: denom ! cached g_tmp + i*Q_e + p^2 + REAL(r8) :: p2, p4 ! cached p^2, p^4 + REAL(r8) :: D2 ! cached D_norm^2 + +c cache powers and shared denominator + p2 = my_p * my_p + p4 = p2 * p2 + D2 = D_norm * D_norm + denom = g_tmp + ifac*Q_e + p2 + +c evaluate coefficients at the current p + fA = p2 / denom + fA_prime = (denom - 2.0d0*p2) / denom ! (g+iQe - p^2)/(g+iQe + p^2) + + fB = g_tmp*(g_tmp + ifac*Q_i) + $ + (g_tmp + ifac*Q_i)*(P_perp + P_tor)*p2 + $ + (P_perp*P_tor)*p4 + + fC = g_tmp + ifac*Q_e + $ + (P_perp + (g_tmp + ifac*Q_i)*D2)*p2 + $ + (1.0d0/iota_e)*P_tor*D2*p4 + + dWdp(1) = -(fA_prime/my_p)*W(1) - W(1)*W(1)/my_p + $ + (fB/(fA*fC))*(p2*my_p) + + RETURN + END SUBROUTINE w_der_f +c----------------------------------------------------------------------- +c W derivative for riccati() c----------------------------------------------------------------------- SUBROUTINE w_der(neq,x,y,dy) @@ -122,42 +510,42 @@ SUBROUTINE w_der(neq,x,y,dy) COMPLEX(r8) :: C2p COMPLEX(r8) :: A1 COMPLEX(r8) :: A2 - COMPLEX(r8), PARAMETER :: ifac=(0,1) + !COMPLEX(r8), PARAMETER :: ifac=(0,1) IF (parflow_flag) THEN C1=((1 + tau)*x**2*pe* $ (-(((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))* - $ (4*ds**2*pr*(1 + tau)*x**3 + + $ (4*ds**2*pr*(1 + tau)*x**3 + $ 2*x*(c_beta**2 + ifac*ds**2*(Q - Q_i)))) - $ /(ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ /(ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) $ **2) + ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (-((4*ds**2*pr*(1 + tau)*x**3 + - $ 2*x*(c_beta**2 + + $ (-((4*ds**2*pr*(1 + tau)*x**3 + + $ 2*x*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ - $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)) + + $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)) + $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i)))) $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**2))/ - $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + + $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)))/ - $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))))/ - $ (ifac*Q + pr*x**2 + x**2*pe - + $ (ifac*Q + pr*x**2 + x**2*pe - $ (ds**2*(1 + tau)*x**6*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -165,106 +553,106 @@ SUBROUTINE w_der(neq,x,y,dy) $ + (ifac*(1 + tau)*x**2*pe*Q_e)/ $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e))) - + C1p=((1 + tau)*x**2*pe* $ ((2*(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))* - $ (4*ds**2*pr*(1 + tau)*x**3 + + $ (4*ds**2*pr*(1 + tau)*x**3 + $ 2*x*(c_beta**2 + ifac*ds**2*(Q - Q_i))) $ **2)/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))** $ 3 - ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))* - $ (12*ds**2*pr*(1 + tau)*x**2 + + $ (12*ds**2*pr*(1 + tau)*x**2 + $ 2*(c_beta**2 + ifac*ds**2*(Q - Q_i))))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))** $ 2 - (2*(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (-((4*ds**2*pr*(1 + tau)*x**3 + + $ (-((4*ds**2*pr*(1 + tau)*x**3 + $ 2*x* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) - $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))+ + $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))+ $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**2)* - $ (4*ds**2*pr*(1 + tau)*x**3 + + $ (4*ds**2*pr*(1 + tau)*x**3 + $ 2*x*(c_beta**2 + ifac*ds**2*(Q - Q_i))))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))** $ 2 - (2*(2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))* - $ (4*ds**2*pr*(1 + tau)*x**3 + + $ (4*ds**2*pr*(1 + tau)*x**3 + $ 2*x*(c_beta**2 + ifac*ds**2*(Q - Q_i))))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))** $ 2 + ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (-((12*ds**2*pr*(1 + tau)*x**2 + + $ (-((12*ds**2*pr*(1 + tau)*x**2 + $ 2*(c_beta**2 + ifac*ds**2*(Q - Q_i)) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)) $ + (2*(2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (4*ds**2*pr*(1 + tau)*x**3 + + $ (4*ds**2*pr*(1 + tau)*x**3 + $ 2*x* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**2 $ - (2*(2*c_beta**2*x + 4*ds**2*pr*tau*x**3)**2* - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**3 $ + ((2*c_beta**2 + 12*ds**2*pr*tau*x**2)* - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**2)) - $ /(ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + + $ /(ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + $ (2*(2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (-((4*ds**2*pr*(1 + tau)*x**3 + + $ (-((4*ds**2*pr*(1 + tau)*x**3 + $ 2*x* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) - $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))+ + $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))+ $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**2)) - $ /(ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + + $ /(ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + $ ((2*c_beta**2 + 12*ds**2*pr*tau*x**2)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))))/ - $ (ifac*Q + pr*x**2 + x**2*pe - + $ (ifac*Q + pr*x**2 + x**2*pe - $ (ds**2*(1 + tau)*x**6*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -273,14 +661,14 @@ SUBROUTINE w_der(neq,x,y,dy) $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e))) $ - ((1 + tau)*x**2*pe* - $ (2*pr*x + 2*x*pe + + $ (2*pr*x + 2*x*pe + $ (ds**2*(1 + tau)*x**6*pe**2* $ (2*x + (4*ds**2*(1 + tau)*x**3*pe)/ $ c_beta**2))/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ - $ c_beta**2 + - $ ifac*(Q - Q_e))**2) - + $ c_beta**2 + + $ ifac*(Q - Q_e))**2) - $ (6*ds**2*(1 + tau)*x**5*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -295,41 +683,41 @@ SUBROUTINE w_der(neq,x,y,dy) $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e))) $ *(-(((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))* - $ (4*ds**2*pr*(1 + tau)*x**3 + + $ (4*ds**2*pr*(1 + tau)*x**3 + $ 2*x*(c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) $ **2) + ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (-((4*ds**2*pr*(1 + tau)*x**3 + + $ (-((4*ds**2*pr*(1 + tau)*x**3 + $ 2*x* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) - $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)) + + $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)) + $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**2)) - $ /(ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + + $ /(ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))))/ - $ (ifac*Q + pr*x**2 + x**2*pe - + $ (ifac*Q + pr*x**2 + x**2*pe - $ (ds**2*(1 + tau)*x**6*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -340,41 +728,41 @@ SUBROUTINE w_der(neq,x,y,dy) $ c_beta**2 + ifac*(Q - Q_e))) $ **2 + (2*(1 + tau)*x*pe* $ (-(((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))* - $ (4*ds**2*pr*(1 + tau)*x**3 + + $ (4*ds**2*pr*(1 + tau)*x**3 + $ 2*x*(c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) $ **2) + ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (-((4*ds**2*pr*(1 + tau)*x**3 + + $ (-((4*ds**2*pr*(1 + tau)*x**3 + $ 2*x* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) - $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))+ + $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))+ $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**2)) - $ /(ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + + $ /(ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))))/ - $ (ifac*Q + pr*x**2 + x**2*pe - + $ (ifac*Q + pr*x**2 + x**2*pe - $ (ds**2*(1 + tau)*x**6*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -382,7 +770,7 @@ SUBROUTINE w_der(neq,x,y,dy) $ + (ifac*(1 + tau)*x**2*pe*Q_e)/ $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e))) - + C2=((1 + tau)*x**2*pe* $ ((ds**2*x**4*pe)/ $ (c_beta**2* @@ -392,13 +780,13 @@ SUBROUTINE w_der(neq,x,y,dy) $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e)) $ + ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)))/ - $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))))/ - $ (ifac*Q + pr*x**2 + x**2*pe - + $ (ifac*Q + pr*x**2 + x**2*pe - $ (ds**2*(1 + tau)*x**6*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -406,15 +794,15 @@ SUBROUTINE w_der(neq,x,y,dy) $ + (ifac*(1 + tau)*x**2*pe*Q_e)/ $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e))) - + C2p=((1 + tau)*x**2*pe* $ (-((ds**2*x**4*pe* $ (2*x + (4*ds**2*(1 + tau)*x**3*pe)/ $ c_beta**2))/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ - $ c_beta**2 + - $ ifac*(Q - Q_e))**2)) + + $ c_beta**2 + + $ ifac*(Q - Q_e))**2)) + $ (4*ds**2*x**3*pe)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -425,40 +813,40 @@ SUBROUTINE w_der(neq,x,y,dy) $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e)) $ **2 - ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4))* - $ (4*ds**2*pr*(1 + tau)*x**3 + + $ (4*ds**2*pr*(1 + tau)*x**3 + $ 2*x*(c_beta**2 + ifac*ds**2*(Q - Q_i))))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))** $ 2 + ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (-((4*ds**2*pr*(1 + tau)*x**3 + + $ (-((4*ds**2*pr*(1 + tau)*x**3 + $ 2*x* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) - $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)) + + $ /(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)) + $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2* $ (c_beta**2 + ifac*ds**2*(Q - Q_i))) $ )/(ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)**2)) - $ /(ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + + $ /(ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))) + $ ((2*c_beta**2*x + 4*ds**2*pr*tau*x**3)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))))/ - $ (ifac*Q + pr*x**2 + x**2*pe - + $ (ifac*Q + pr*x**2 + x**2*pe - $ (ds**2*(1 + tau)*x**6*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -467,14 +855,14 @@ SUBROUTINE w_der(neq,x,y,dy) $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e))) $ - ((1 + tau)*x**2*pe* - $ (2*pr*x + 2*x*pe + + $ (2*pr*x + 2*x*pe + $ (ds**2*(1 + tau)*x**6*pe**2* $ (2*x + (4*ds**2*(1 + tau)*x**3*pe)/ $ c_beta**2))/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ - $ c_beta**2 + - $ ifac*(Q - Q_e))**2) - + $ c_beta**2 + + $ ifac*(Q - Q_e))**2) - $ (6*ds**2*(1 + tau)*x**5*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -496,15 +884,15 @@ SUBROUTINE w_der(neq,x,y,dy) $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e)) $ + ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))))/ - $ (ifac*Q + pr*x**2 + x**2*pe - + $ (ifac*Q + pr*x**2 + x**2*pe - $ (ds**2*(1 + tau)*x**6*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -522,15 +910,15 @@ SUBROUTINE w_der(neq,x,y,dy) $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e)) $ + ((ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)* - $ (1 - (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + - $ x**2*(c_beta**2 + + $ (1 - (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + + $ x**2*(c_beta**2 + $ ifac*ds**2*(Q - Q_i)))/ $ (ifac*Q + c_beta**2*x**2 + ds**2*pr*tau*x**4)))/ - $ (ds**2*pr*(1 + tau)*x**4 + - $ ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + + $ ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i)))))/ - $ (ifac*Q + pr*x**2 + x**2*pe - + $ (ifac*Q + pr*x**2 + x**2*pe - $ (ds**2*(1 + tau)*x**6*pe**2)/ $ (c_beta**2* $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -543,24 +931,24 @@ SUBROUTINE w_der(neq,x,y,dy) C1p=0 C2=0 C2p=0 - ENDIF - + ENDIF + IF (PeOhmOnly_flag) THEN - G=((c_beta**2*pr*x**4 - Q*(Q - Q_i) + + G=((c_beta**2*pr*x**4 - Q*(Q - Q_i) + $ ifac*(c_beta**2 + pr)*x**2*(Q - Q_i))/ - $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))))*(x**2.0) ELSE - G=(x**2*pe + - $ (c_beta**2*pr*x**4 - Q*(Q - Q_i) + + G=(x**2*pe + + $ (c_beta**2*pr*x**4 - Q*(Q - Q_i) + $ ifac*(c_beta**2 + pr)*x**2*(Q - Q_i))/ - $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + + $ (ds**2*pr*(1 + tau)*x**4 + ifac*(Q - Q_e) + $ x**2*(c_beta**2 + ifac*ds**2*(Q - Q_i))))*(x**2.0) ENDIF - + C3=x**2/(x**2 + (ds**2*(1 + tau)*x**4*pe)/ $ c_beta**2 + ifac*(Q - Q_e)) - + C3p=-((x**2*(2*x + (4*ds**2*(1 + tau)*x**3*pe)/ $ c_beta**2))/ $ (x**2 + (ds**2*(1 + tau)*x**4*pe)/ @@ -570,11 +958,11 @@ SUBROUTINE w_der(neq,x,y,dy) $ c_beta**2 + (0,1)*(Q - Q_e)) A1=(C1 + (C3p/C3)*(C2 + 1) + C2p)/(C2 + 1) - + A2=(C1p + C1*(C3p/C3) - G/C3)/(C2 + 1) - + dy(1)=(-A1 + 1/x)*y(1) - y(1)*y(1)/x - A2*x - + RETURN END SUBROUTINE w_der c----------------------------------------------------------------------- @@ -586,7 +974,7 @@ SUBROUTINE w_der_temp(neq,x,y,dy) REAL(r8), INTENT(IN) :: x COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: y COMPLEX(r8), DIMENSION(neq), INTENT(OUT) :: dy - COMPLEX(r8), PARAMETER :: ifac=(0,1) + !COMPLEX(r8), PARAMETER :: ifac=(0,1) dy(1)=(2.0*x/(ifac*(Q-Q_e)+x**2.0)-1.0/x)*y(1)-y(1)*y(1)/x $ +x*(ifac*(Q-Q_e)+x**2.0) @@ -619,7 +1007,7 @@ FUNCTION directred(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) INTEGER, DIMENSION(:), ALLOCATABLE :: iwork REAL(r8), DIMENSION(:), ALLOCATABLE :: xfac,atol,rwork - + Q=inQ Q_e=inQ_e Q_i=inQ_i @@ -638,10 +1026,10 @@ FUNCTION directred(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) istate = 1 iopt = 0 mf = 10 - liw = 20 + liw = 20 lrw = 22+16*neq ALLOCATE(iwork(liw),rwork(lrw)) - + xintv = 0.1 x=0.1 istep=1 @@ -668,12 +1056,12 @@ END FUNCTION directred c direct integration (obsolete). c----------------------------------------------------------------------- SUBROUTINE phi_der(neq,x,y,dy) - + INTEGER, INTENT(IN) :: neq REAL(r8), INTENT(IN) :: x COMPLEX(r8), DIMENSION(neq), INTENT(IN) :: y COMPLEX(r8), DIMENSION(neq), INTENT(OUT) :: dy - COMPLEX(r8), PARAMETER :: ifac=(0,1) + !COMPLEX(r8), PARAMETER :: ifac=(0,1) dy(1)=(1+ifac*(Q-Q_e)*x**2.0)/x**2.0*y(2) dy(2)=(-Q*(Q-Q_i)*x**4.0+ifac*(Q-Q_i)*(pr+c_beta**2.0)*x**2.0 @@ -681,6 +1069,6 @@ SUBROUTINE phi_der(neq,x,y,dy) $ +(c_beta**2.0+ifac*(Q-Q_i)*ds**2.0)*x**6.0 $ +(1+tau)*pr*ds**2.0*x**4.0)*y(1) RETURN - END SUBROUTINE phi_der - - END MODULE delta_mod + END SUBROUTINE phi_der + + END MODULE delta_mod \ No newline at end of file diff --git a/slayer/growthrates.f b/slayer/growthrates.f new file mode 100644 index 00000000..502fffae --- /dev/null +++ b/slayer/growthrates.f @@ -0,0 +1,716 @@ + MODULE growthrates_mod +c----------------------------------------------------------------------- +c growthrates_mod: Growth-rate scanning and AMR dispersion solvers. +c +c Split from gslayer_mod, this module contains every subroutine +c after gpec_slayer that was formerly inside gslayer.f. It +c provides I/O helpers, array utilities, the dispersion +c determinant, and both AMR scanner variants used by the main +c SLAYER driver (slayer.f). +c +c Subprograms: +c 1. output_gamma - write results to netCDF +c 2. allocate_inputs - allocate slayer_inputs_type +c 3. allocate_outputs - allocate slayer_outputs_type +c 4. shrink_array - trim over-allocated scan arrays +c 5. grow_array - expand scan arrays dynamically +c 6. calc_determinant - 2x2 / 3x3 complex determinant +c 7. dispersion_det - coupled dispersion determinant +c 8. get_or_compute_v2 - hash-cached dispersion eval +c 9. dispersion_AMR_v2 - AMR scan (cell-based storage) +c 10. check_cell_crossing_sub - zero-crossing test +c 11. subdivide_cell_sub - cell refinement +c----------------------------------------------------------------------- + + USE omp_lib + + USE sglobal_mod, ONLY: out_unit,r8,mu0,m_p,chag,lnLamb, + $ Q_e,Q_i,pr,pe,c_beta,ds,tau, + $ eta,visc,rho_s,lu,omega_e,omega_i, + $ delta_n,Q, + $ ifac,g_tmp,pi, ! used by AMR + $ tauk,iota_e,D_norm,P_perp,P_tor,delta_eff, ! used by det + $ amr_cell_type,amr_cells,n_amr_cells, ! v2 types + $ Q_store,D_store,n_pts, ! output arrays + $ MAX_PTS,HASH_SZ,HASH_SCALE, ! v1 constants + $ hash_head,hash_next, ! v1 hash + $ MAX_CELLS, ! v2 limit + $ slayer_inputs_type,slayer_outputs_type, + $ deltas_outputs_type, + $ tau_r,dc_tmp,dc_type, + $ sn_str,sm_str + USE delta_mod, ONLY: riccati,riccati_f,riccati_out, + $ parflow_flag,PeOhmOnly_flag + USE params_mod + USE layerinputs_mod + USE slayer_netcdf_mod + + IMPLICIT NONE + +c --- reconnection regulariser used in psi0 / JxB expressions; +c expose via namelist to make user-configurable. + REAL(r8), PARAMETER :: DELTA_N_PERT = 1.0e-2_r8 + + CONTAINS + +c----------------------------------------------------------------------- +c subprogram 1. output_gamma. +c Write growth-rate results to netCDF via slayer_netcdf_out. +c Passes the input/output structured types and AMR results +c through to the netCDF writer (slayer_netcdf_mod). +c----------------------------------------------------------------------- + SUBROUTINE output_gamma(est_gamma_flag,m_AMR,sl_in,sl_out, + $ all_deltas_out) + + LOGICAL, INTENT(IN) :: est_gamma_flag ! single-surface mode? + INTEGER, INTENT(IN) :: m_AMR ! AMR pass count + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + TYPE(slayer_outputs_type), INTENT(IN) :: sl_out + TYPE(deltas_outputs_type), INTENT(IN) :: + $ all_deltas_out(SIZE(sl_in%qval_arr)) + + CALL slayer_netcdf_out(SIZE(sl_in%qval_arr),m_AMR,est_gamma_flag, + $ sl_in,sl_out,all_deltas_out) + + END SUBROUTINE output_gamma +c----------------------------------------------------------------------- +c subprogram 2. allocate_inputs. +c Allocate all per-surface arrays inside slayer_inputs_type. +c----------------------------------------------------------------------- + SUBROUTINE allocate_inputs(n_k,sl_in) + + INTEGER, INTENT(IN) :: n_k ! number of surfaces + TYPE(slayer_inputs_type), INTENT(INOUT) :: sl_in + + ALLOCATE(sl_in%qval_arr(n_k),sl_in%omegas_arr(n_k), + $ sl_in%Q_e_arr(n_k),sl_in%Q_i_arr(n_k),sl_in%psi_n_arr(n_k), + $ sl_in%Re_dp_arr(n_k),sl_in%Im_dp_arr(n_k), + $ sl_in%d_crit_arr(n_k),sl_in%P_tor_arr(n_k), + $ sl_in%P_perp_arr(n_k),sl_in%tau_arr(n_k), + $ sl_in%D_norm_arr(n_k), + $ sl_in%d_beta_arr(n_k),sl_in%gammafac_arr(n_k), + $ sl_in%c_beta_arr(n_k),sl_in%lu_arr(n_k),sl_in%Qconv_arr(n_k)) + RETURN + END SUBROUTINE allocate_inputs +c----------------------------------------------------------------------- +c subprogram 3. allocate_outputs. +c Allocate per-surface arrays inside slayer_outputs_type. +c----------------------------------------------------------------------- + SUBROUTINE allocate_outputs(n_k,sl_out) + + INTEGER, INTENT(IN) :: n_k ! number of surfaces + TYPE(slayer_outputs_type), INTENT(INOUT) :: sl_out + + ALLOCATE(sl_out%dels_db_arr(n_k),sl_out%gamma_sol_arr(n_k), + $ sl_out%gamma_est_arr(n_k),sl_out%br_th_arr(n_k) ) + RETURN + END SUBROUTINE allocate_outputs +c----------------------------------------------------------------------- +c subprogram 4. shrink_array. +c Trim an over-allocated REAL(r8) array down to new_size using +c MOVE_ALLOC (no copy of trailing elements). +c----------------------------------------------------------------------- + SUBROUTINE shrink_array(arr, new_size) + + REAL(r8), ALLOCATABLE, INTENT(INOUT) :: arr(:) + INTEGER, INTENT(IN) :: new_size ! target size + REAL(r8), ALLOCATABLE :: temp(:) ! temporary buffer + + ALLOCATE(temp(new_size)) + temp(1:new_size) = arr(1:new_size) + CALL move_alloc(temp, arr) + END SUBROUTINE shrink_array +c----------------------------------------------------------------------- +c subprogram 5. grow_array. +c Expand a REAL(r8) array from old_size to new_size, preserving +c existing data via MOVE_ALLOC. +c----------------------------------------------------------------------- + SUBROUTINE grow_array(arr, old_size, new_size) + + REAL(r8), ALLOCATABLE, INTENT(INOUT) :: arr(:) + INTEGER, INTENT(IN) :: old_size ! current valid element count + INTEGER, INTENT(IN) :: new_size ! target allocation size + REAL(r8), ALLOCATABLE :: temp(:) ! temporary buffer + + ALLOCATE(temp(new_size)) + temp(1:old_size) = arr(1:old_size) + CALL move_alloc(temp, arr) + END SUBROUTINE grow_array +c----------------------------------------------------------------------- +c subprogram 6. calc_determinant. +c Compute the determinant of a 2x2 or 3x3 complex matrix. +c Returns (0,0) and sets status=-1 for unsupported sizes. +c status=0 on success, -1 when nk is neither 2 nor 3. +c----------------------------------------------------------------------- + SUBROUTINE calc_determinant(matk, nk, detk, status) + + IMPLICIT NONE + +c --- arguments + INTEGER, INTENT(IN) :: nk ! matrix rank (2 or 3) + COMPLEX(r8), DIMENSION(nk,nk), INTENT(IN) :: matk ! input matrix + COMPLEX(r8), INTENT(OUT) :: detk ! determinant result + INTEGER, INTENT(OUT) :: status ! 0=success, -1=unsupported rank + + status = 0 ! Initialize status as success + + SELECT CASE (nk) + CASE (2) + ! 2x2 determinant: ad - bc + detk = matk(1,1) * matk(2,2) - matk(1,2) * matk(2,1) + + CASE (3) + ! 3x3 determinant using cofactor expansion along first row + detk = matk(1,1)*(matk(2,2)*matk(3,3)-matk(2,3) + $ *matk(3,2))-matk(1,2)*(matk(2,1)*matk(3,3) + $ -matk(2,3)*matk(3,1))+matk(1,3)*(matk(2,1) + $ *matk(3,2)-matk(2,2)*matk(3,1)) + + CASE default + ! Unsupported matrix size + detk = CMPLX(0.0_r8, 0.0_r8, KIND=r8) + status = -1 + + END SELECT + RETURN + END SUBROUTINE calc_determinant +c----------------------------------------------------------------------- +c subprogram 7. dispersion_det. +c Compute the coupled dispersion determinant for n_k surfaces. +c +c For n_k = 1 (single surface): +c Evaluate riccati_f() (uses module g_tmp), de-normalise by lu^(1/3), and +c return Deltaprime - delta(Q). +c +c For n_k = 2 or 3 (coupled surfaces): +c Build the diagonal delta(Q) matrix, subtract from dp_matrix, +c and return det(dp_matrix - delta_Q). +c +c----------------------------------------------------------------------- + FUNCTION dispersion_det(g_in,n_k,sl_in,msing_max) + +c --- arguments + COMPLEX(r8), INTENT(IN) :: g_in ! complex growth rate + INTEGER, INTENT(IN) :: n_k ! number of surfaces + INTEGER, INTENT(IN) :: msing_max ! max surfaces to include + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in +c --- function result and locals + COMPLEX(r8) :: dispersion_det ! returned determinant + COMPLEX(r8) :: det_val ! intermediate determinant + COMPLEX(r8) :: tmp_delta ! single-surface delta + COMPLEX(r8), ALLOCATABLE :: delta_Q(:,:) ! diagonal delta matrix + COMPLEX(r8), ALLOCATABLE :: result_matrix(:,:) ! dp - delta_Q + INTEGER :: k ! surface loop index + INTEGER :: det_status ! calc_determinant error status + +c --- single-surface branch + IF (msing_max < 2) THEN +c set module-level variables for riccati_f + Q_e = sl_in%Q_e_arr(1) + Q_i = sl_in%Q_i_arr(1) + P_perp = sl_in%P_perp_arr(1) + P_tor = sl_in%P_tor_arr(1) + tau = sl_in%tau_arr(1) + D_norm = sl_in%D_norm_arr(1) + c_beta = sl_in%c_beta_arr(1) + tauk = sl_in%Qconv_arr(1) + iota_e = Q_e / (Q_e - Q_i) +c BENCHMARK PATCH: ensure pr is set so riccati_f's +c atol = 1e-7*pr**0.4 does not inherit a stale/zero module value. + pr = P_perp + + g_tmp = g_in + tmp_delta=riccati_f() +c de-normalise delta by lu^(1/3) + det_val=tmp_delta* + $ (sl_in%lu_arr(1)**(1.0_r8/3.0_r8)) + +c return Deltaprime - delta(Q) + dispersion_det = sl_in%Re_dp_arr(1) - det_val + +c --- coupled-surface branch (2 or 3 surfaces) + ELSEIF ((msing_max == 2) .OR. (msing_max == 3)) THEN + ALLOCATE(delta_Q(msing_max,msing_max)) + delta_Q=CMPLX(0.0_r8, 0.0_r8, KIND=r8) + DO k=1,msing_max +c set module-level variables for this surface + Q_e = sl_in%Q_e_arr(k) + Q_i = sl_in%Q_i_arr(k) + P_perp = sl_in%P_perp_arr(k) + P_tor = sl_in%P_tor_arr(k) + tau = sl_in%tau_arr(k) + D_norm = sl_in%D_norm_arr(k) + c_beta = sl_in%c_beta_arr(k) + tauk = sl_in%Qconv_arr(k) + iota_e = Q_e / (Q_e - Q_i) +c BENCHMARK PATCH: ensure pr is set so riccati_f's +c atol = 1e-7*pr**0.4 does not inherit a stale/zero module value. + pr = P_perp + +c evaluate riccati_f at rescaled growth rate, de-normalise +c rescale g_in to this surface's normalisation + g_tmp = (g_in*sl_in%Qconv_arr(1))/tauk + delta_Q(k,k)=riccati_f() + delta_Q(k,k)=delta_Q(k,k)* + $ sl_in%lu_arr(k)**(1.0_r8/3.0_r8) + END DO + +c Coupled dispersion: det(Delta' - diag(Delta_c) - delta_Q). +c Delta_crit is a *local* layer-intrinsic offset at each +c rational surface, so it modifies only the diagonal of +c Delta'; outer-region off-diagonal couplings are untouched. +c This reduces to the single-surface matching +c (Delta'_kk - Delta_crit_k) = S^(1/3) * Delta_s +c on the diagonal. When dc_type='none', d_crit_arr is zero +c and the subtraction is a no-op, so this is backward +c compatible with the original coupled behaviour. +c +c BENCHMARK PATCH: take only the (1:msing_max, 1:msing_max) +c sub-block of dp_matrix so callers can compute the coupled +c dispersion on a subset of the rational surfaces that STRIDE +c actually found. This is needed when STRIDE finds 4+ surfaces +c (e.g. DIII-D 147131 with sas_flag=t, qhigh=1e3) but the +c coupled-dispersion code only supports msing_max ≤ 3. + result_matrix = sl_in%dp_matrix(1:msing_max,1:msing_max) + $ - delta_Q + DO k = 1, msing_max + result_matrix(k,k) = result_matrix(k,k) + $ - CMPLX(sl_in%d_crit_arr(k), 0.0_r8, KIND=r8) + END DO + CALL calc_determinant(result_matrix, msing_max, det_val, + $ det_status) + IF (det_status /= 0) THEN + WRITE(*,*) 'ERROR: calc_determinant unsupported rank=', + $ msing_max + STOP 'dispersion_det: calc_determinant failed' + END IF + DEALLOCATE(delta_Q, result_matrix) + dispersion_det = det_val + ELSE + WRITE(*,*) "Error: no support for msing > 3" + STOP "dispersion_det: unsupported n_k" + END IF + END FUNCTION dispersion_det + +c----------------------------------------------------------------------- +c get_or_compute_v2: hash-cached dispersion evaluation for AMR v2. +c Applies ifac Wick rotation: g_tmp = q_in * ifac. +c----------------------------------------------------------------------- + SUBROUTINE get_or_compute_v2(q_in, idx_out, n_k, + $ sl_in, msing_max, + $ coupling_flag, full) + + IMPLICIT NONE + +c --- arguments + COMPLEX(r8), INTENT(IN) :: q_in + INTEGER, INTENT(OUT) :: idx_out + INTEGER, INTENT(IN) :: n_k, msing_max + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + LOGICAL, INTENT(IN) :: coupling_flag + LOGICAL, INTENT(OUT) :: full + +c --- locals + INTEGER :: h, curr + COMPLEX(r8) :: delta_val + INTEGER(8) :: ix8, iy8, h8 + + full = .FALSE. +c --- 1. compute hash bucket + ix8 = NINT(REAL(q_in) * HASH_SCALE, KIND=8) + iy8 = NINT(AIMAG(q_in) * HASH_SCALE, KIND=8) + h8 = MOD(ABS(ix8 * 73856093_8 + iy8 * 19349663_8), + $ INT(HASH_SZ, 8)) + 1_8 + h = INT(h8) + +c --- 2. search hash chain for existing point + curr = hash_head(h) + DO WHILE (curr /= 0) + IF (ABS(Q_store(curr) - q_in) < 1.0d-8) THEN + idx_out = curr + RETURN + END IF + curr = hash_next(curr) + END DO + +c --- 3. not found: evaluate with ifac rotation and store + n_pts = n_pts + 1 + IF (n_pts > MAX_PTS) THEN + n_pts = n_pts - 1 + full = .TRUE. + idx_out = -1 + RETURN + END IF + + idx_out = n_pts + Q_store(idx_out) = q_in + + IF (coupling_flag) THEN +c dispersion_det sets g_tmp per-surface internally; +c pass q_in*ifac directly as g_in argument. + delta_val = dispersion_det(q_in * ifac, n_k, sl_in, + $ msing_max) + ELSE + g_tmp = q_in * ifac + delta_val = riccati_f() + delta_val = delta_val - delta_eff + END IF + D_store(idx_out) = delta_val + +c --- 4. insert into hash chain (prepend) + hash_next(idx_out) = hash_head(h) + hash_head(h) = idx_out + + END SUBROUTINE get_or_compute_v2 + +c----------------------------------------------------------------------- +c dispersion_AMR_v2: cell-based adaptive mesh refinement scanner. +c Unlike v1 (hash-based point deduplication), v2 stores complete +c cells (TYPE amr_cell_type), each carrying 4 corner Q- and D- +c values. Refinement subdivides cells that contain a zero in +c Re(D) or Im(D) and re-evaluates the dispersion relation at the +c 5 new midpoints. +c +c All dispersion evaluations go through get_or_compute_v2, which +c caches results in Q_store / D_store via a hash table. This +c eliminates redundant evaluations for shared corners (initial grid) +c and shared edge-midpoints (refinement). At completion, Q_store +c and D_store are trimmed to n_pts unique output points. +c +c----------------------------------------------------------------------- + SUBROUTINE dispersion_AMR_v2(n_k, sl_in, msing_max, + $ scan_width, Q_num, AMR_passes, + $ coupling_flag) + + IMPLICIT NONE + +c --- arguments + INTEGER, INTENT(IN) :: n_k ! number of rational surfaces + INTEGER, INTENT(IN) :: msing_max ! max surfaces for coupling + INTEGER, INTENT(IN) :: Q_num ! grid points per axis + INTEGER, INTENT(IN) :: AMR_passes ! refinement passes + REAL(r8), INTENT(IN) :: scan_width ! half-width of Re/Im scan window + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + LOGICAL, INTENT(IN) :: coupling_flag ! coupled dispersion_det? +c --- locals + TYPE(amr_cell_type), ALLOCATABLE :: new_cells(:) + TYPE(amr_cell_type), ALLOCATABLE :: swap_tmp(:) ! for pointer swap + INTEGER :: i, j, c, corner, pass ! loop counters + REAL(r8) :: step ! grid spacing + REAL(r8) :: x, y ! real / imag grid coords + LOGICAL :: cross_real, cross_imag ! zero-crossing flags + LOGICAL :: pts_full ! MAX_PTS reached flag + INTEGER :: n_new_cells ! count during refinement + INTEGER :: cells_to_refine ! cells flagged per pass + INTEGER :: cells_kept ! cells kept per pass + INTEGER :: idx_tmp ! hash-cache index + COMPLEX(r8), ALLOCATABLE :: temp_Q(:) ! for trimming output + COMPLEX(r8), ALLOCATABLE :: temp_D(:) ! for trimming output + +c --- 1. initialise cell storage and hash cache + + IF (ALLOCATED(amr_cells)) DEALLOCATE(amr_cells) + IF (ALLOCATED(Q_store)) DEALLOCATE(Q_store) + IF (ALLOCATED(D_store)) DEALLOCATE(D_store) + IF (ALLOCATED(hash_head)) DEALLOCATE(hash_head) + IF (ALLOCATED(hash_next)) DEALLOCATE(hash_next) + + ALLOCATE(amr_cells(MAX_CELLS)) + ALLOCATE(new_cells(MAX_CELLS)) + ALLOCATE(Q_store(MAX_PTS)) + ALLOCATE(D_store(MAX_PTS)) + ALLOCATE(hash_head(HASH_SZ)) + ALLOCATE(hash_next(MAX_PTS)) + + hash_head = 0 + hash_next = 0 + n_pts = 0 + n_amr_cells = 0 + step = (2.0d0 * scan_width) / DBLE(Q_num - 1) + +c --- 2. build initial coarse grid of (Q_num-1)^2 cells + + DO i = 1, Q_num - 1 + DO j = 1, Q_num - 1 + x = -scan_width + DBLE(i-1) * step + y = -scan_width + DBLE(j-1) * step + + n_amr_cells = n_amr_cells + 1 + + IF (n_amr_cells > MAX_CELLS) THEN + WRITE(*,*) 'ERROR: Exceeded MAX_CELLS in init' + STOP 'dispersion_AMR_v2: MAX_CELLS in init' + END IF + +c corner order: BL=1, BR=2, TL=3, TR=4 + amr_cells(n_amr_cells)%Q(1) = CMPLX(x, y, KIND=r8) + amr_cells(n_amr_cells)%Q(2) = CMPLX(x+step, y, KIND=r8) + amr_cells(n_amr_cells)%Q(3) = CMPLX(x, y+step, KIND=r8) + amr_cells(n_amr_cells)%Q(4) = CMPLX(x+step, y+step, + $ KIND=r8) + +c evaluate dispersion at each corner (hash-cached) + DO corner = 1, 4 + CALL get_or_compute_v2( + $ amr_cells(n_amr_cells)%Q(corner), + $ idx_tmp, n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) GOTO 800 + amr_cells(n_amr_cells)%D(corner) = + $ D_store(idx_tmp) + END DO + + amr_cells(n_amr_cells)%needs_refine = .FALSE. + END DO + END DO + +c --- 3. refinement passes: subdivide cells with zero crossings + DO pass = 1, AMR_passes + WRITE(*,'(A,I2,A,I7,A)') ' Pass ', pass, + $ ': Processing ', n_amr_cells, ' cells' + +c flag cells that span a zero in Re(D) or Im(D) + cells_to_refine = 0 + DO c = 1, n_amr_cells + CALL check_cell_crossing_sub(amr_cells(c), + $ cross_real, cross_imag) + amr_cells(c)%needs_refine = (cross_real .OR. cross_imag) + IF (amr_cells(c)%needs_refine) THEN + cells_to_refine = cells_to_refine + 1 + END IF + END DO + +c build new cell list: subdivide flagged, keep the rest + n_new_cells = 0 + cells_kept = 0 + + DO c = 1, n_amr_cells + IF (amr_cells(c)%needs_refine) THEN + CALL subdivide_cell_sub( + $ amr_cells(c), + $ new_cells, n_new_cells, + $ MAX_CELLS, n_k, sl_in, + $ msing_max, coupling_flag, + $ pts_full) + IF (pts_full) GOTO 800 + ELSE + n_new_cells = n_new_cells + 1 + IF (n_new_cells > MAX_CELLS) THEN + WRITE(*,*) 'ERROR: Exceeded MAX_CELLS in refine' + STOP 'dispersion_AMR_v2: MAX_CELLS in refine' + END IF + new_cells(n_new_cells) = amr_cells(c) + cells_kept = cells_kept + 1 + END IF + END DO + +c swap arrays for next pass (pointer swap, no element copy) + CALL MOVE_ALLOC(new_cells, swap_tmp) + CALL MOVE_ALLOC(amr_cells, new_cells) ! old amr_cells becomes new_cells + CALL MOVE_ALLOC(swap_tmp, amr_cells) ! filled array becomes amr_cells + n_amr_cells = n_new_cells + + END DO + GOTO 810 + + 800 CONTINUE + WRITE(*,'(A,I8,A)') + $ ' WARNING: MAX_PTS (', MAX_PTS, + $ ') reached during AMR v2.' + WRITE(*,'(A)') + $ ' Saving existing results.' + + 810 CONTINUE +c --- 4. output: Q_store/D_store already populated by hash cache. +c Trim to exact size n_pts and deallocate hash infrastructure. + + ALLOCATE(temp_Q(n_pts)) + ALLOCATE(temp_D(n_pts)) + temp_Q(1:n_pts) = Q_store(1:n_pts) + temp_D(1:n_pts) = D_store(1:n_pts) + CALL MOVE_ALLOC(temp_Q, Q_store) + CALL MOVE_ALLOC(temp_D, D_store) + + IF (ALLOCATED(hash_head)) DEALLOCATE(hash_head) + IF (ALLOCATED(hash_next)) DEALLOCATE(hash_next) + IF (ALLOCATED(new_cells)) DEALLOCATE(new_cells) +c keep amr_cells allocated for potential post-run inspection + + WRITE(*,*) 'AMR v2 Complete. Unique output points:', n_pts + WRITE(*,'(A,2ES14.6)') ' D_store checksum (Re,Im):', + $ SUM(REAL(D_store(1:n_pts))), + $ SUM(AIMAG(D_store(1:n_pts))) + WRITE(*,'(A,2ES14.6)') ' D_store(1) sample:', + $ REAL(D_store(1)), AIMAG(D_store(1)) + + RETURN + END SUBROUTINE dispersion_AMR_v2 + +c----------------------------------------------------------------------- +c check_cell_crossing_sub: test whether a cell’s 4 corner D-values +c span a zero crossing in Re(D) and/or Im(D). Used by +c dispersion_AMR_v2 to decide which cells to refine. +c----------------------------------------------------------------------- + SUBROUTINE check_cell_crossing_sub(cell, cross_real, cross_imag) + + IMPLICIT NONE + + TYPE(amr_cell_type), INTENT(IN) :: cell + LOGICAL, INTENT(OUT) :: cross_real, cross_imag + + REAL(r8) :: r_vals(4), i_vals(4) ! corner Re/Im values + REAL(r8) :: r_min, r_max, i_min, i_max + INTEGER :: k + + DO k = 1, 4 + r_vals(k) = REAL(cell%D(k), KIND=r8) + i_vals(k) = AIMAG(cell%D(k)) + END DO + + r_min = MINVAL(r_vals) + r_max = MAXVAL(r_vals) + cross_real = (r_min * r_max <= 0.0d0) + + i_min = MINVAL(i_vals) + i_max = MAXVAL(i_vals) + cross_imag = (i_min * i_max <= 0.0d0) + + RETURN + END SUBROUTINE check_cell_crossing_sub + + +c----------------------------------------------------------------------- +c subdivide_cell_sub: split a parent cell into 4 child cells by +c computing 5 midpoints (bottom-mid, top-mid, left-mid, right-mid, +c centre) and evaluating the dispersion relation at each. The 4 +c resulting child cells are appended to new_cells. +c----------------------------------------------------------------------- + SUBROUTINE subdivide_cell_sub(parent, + $ new_cells, n_new, max_cells, n_k, + $ sl_in, msing_max, coupling_flag, + $ pts_full) + + IMPLICIT NONE + +c --- arguments + TYPE(amr_cell_type), INTENT(IN) :: parent + TYPE(amr_cell_type), INTENT(INOUT) :: new_cells(*) + INTEGER, INTENT(INOUT) :: n_new ! running count of new cells + INTEGER, INTENT(IN) :: max_cells ! capacity of new_cells + INTEGER, INTENT(IN) :: n_k + INTEGER, INTENT(IN) :: msing_max + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + LOGICAL, INTENT(IN) :: coupling_flag + LOGICAL, INTENT(OUT) :: pts_full ! MAX_PTS flag +c --- corner coordinates and D-values from parent + COMPLEX(r8) :: q_bl, q_br, q_tl, q_tr + COMPLEX(r8) :: d_bl, d_br, d_tl, d_tr +c --- midpoint coordinates and D-values (cached via hash) + COMPLEX(r8) :: q_bm, q_tm, q_lm, q_rm, q_mm + COMPLEX(r8) :: d_bm, d_tm, d_lm, d_rm, d_mm + INTEGER :: idx_tmp ! hash-cache index + +c --- extract parent corners (BL=1, BR=2, TL=3, TR=4) + q_bl = parent%Q(1) + q_br = parent%Q(2) + q_tl = parent%Q(3) + q_tr = parent%Q(4) + + d_bl = parent%D(1) + d_br = parent%D(2) + d_tl = parent%D(3) + d_tr = parent%D(4) + +c --- compute 5 midpoint coordinates + q_bm = 0.5d0 * (q_bl + q_br) + q_tm = 0.5d0 * (q_tl + q_tr) + q_lm = 0.5d0 * (q_bl + q_tl) + q_rm = 0.5d0 * (q_br + q_tr) + q_mm = 0.25d0 * (q_bl + q_br + q_tl + q_tr) + +c --- evaluate dispersion at new midpoints (hash-cached) + pts_full = .FALSE. + CALL get_or_compute_v2(q_bm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_bm = D_store(idx_tmp) + CALL get_or_compute_v2(q_tm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_tm = D_store(idx_tmp) + CALL get_or_compute_v2(q_lm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_lm = D_store(idx_tmp) + CALL get_or_compute_v2(q_rm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_rm = D_store(idx_tmp) + CALL get_or_compute_v2(q_mm, idx_tmp, + $ n_k, sl_in, msing_max, + $ coupling_flag, pts_full) + IF (pts_full) RETURN + d_mm = D_store(idx_tmp) + +c --- check space for 4 new cells + IF (n_new + 4 > max_cells) THEN + WRITE(*,*) 'ERROR: Would exceed MAX_CELLS in subdivide' + STOP 'subdivide_cell_sub: MAX_CELLS exceeded' + END IF + +c --- child 1: bottom-left quadrant (BL, BM, LM, MM) + n_new = n_new + 1 + new_cells(n_new)%Q(1) = q_bl + new_cells(n_new)%Q(2) = q_bm + new_cells(n_new)%Q(3) = q_lm + new_cells(n_new)%Q(4) = q_mm + new_cells(n_new)%D(1) = d_bl + new_cells(n_new)%D(2) = d_bm + new_cells(n_new)%D(3) = d_lm + new_cells(n_new)%D(4) = d_mm + new_cells(n_new)%needs_refine = .FALSE. + +c --- child 2: bottom-right quadrant (BM, BR, MM, RM) + n_new = n_new + 1 + new_cells(n_new)%Q(1) = q_bm + new_cells(n_new)%Q(2) = q_br + new_cells(n_new)%Q(3) = q_mm + new_cells(n_new)%Q(4) = q_rm + new_cells(n_new)%D(1) = d_bm + new_cells(n_new)%D(2) = d_br + new_cells(n_new)%D(3) = d_mm + new_cells(n_new)%D(4) = d_rm + new_cells(n_new)%needs_refine = .FALSE. + +c --- child 3: top-left quadrant (LM, MM, TL, TM) + n_new = n_new + 1 + new_cells(n_new)%Q(1) = q_lm + new_cells(n_new)%Q(2) = q_mm + new_cells(n_new)%Q(3) = q_tl + new_cells(n_new)%Q(4) = q_tm + new_cells(n_new)%D(1) = d_lm + new_cells(n_new)%D(2) = d_mm + new_cells(n_new)%D(3) = d_tl + new_cells(n_new)%D(4) = d_tm + new_cells(n_new)%needs_refine = .FALSE. + +c --- child 4: top-right quadrant (MM, RM, TM, TR) + n_new = n_new + 1 + new_cells(n_new)%Q(1) = q_mm + new_cells(n_new)%Q(2) = q_rm + new_cells(n_new)%Q(3) = q_tm + new_cells(n_new)%Q(4) = q_tr + new_cells(n_new)%D(1) = d_mm + new_cells(n_new)%D(2) = d_rm + new_cells(n_new)%D(3) = d_tm + new_cells(n_new)%D(4) = d_tr + new_cells(n_new)%needs_refine = .FALSE. + + RETURN + END SUBROUTINE subdivide_cell_sub + + END MODULE growthrates_mod diff --git a/slayer/gslayer.f b/slayer/gslayer.f index 0a4407c9..1275af91 100644 --- a/slayer/gslayer.f +++ b/slayer/gslayer.f @@ -1,106 +1,188 @@ MODULE gslayer_mod - - USE sglobal_mod, ONLY: out_unit, r8, mu0, m_p, chag, lnLamb, +c----------------------------------------------------------------------- +c gslayer_mod: Single-surface SLAYER driver called by GPEC. +c +c This module contains the gpec_slayer subroutine, which computes +c Delta, torque, and critical field threshold for one rational +c surface. It is the primary interface used by gpec/gpout.f. +c +c Growth-rate scanning, AMR dispersion solvers, and I/O helpers +c have been split out into growthrates_mod (growthrates.f). +c +c Subprograms: +c 1. gpec_slayer - single-surface delta/torque driver +c----------------------------------------------------------------------- + + USE omp_lib + + USE sglobal_mod, ONLY: out_unit,r8,mu0,m_p,chag,lnLamb, $ Q_e,Q_i,pr,pe,c_beta,ds,tau, $ eta,visc,rho_s,lu,omega_e,omega_i, - $ delta_n, - $ Q - USE delta_mod, ONLY: riccati,riccati_out, + $ delta_n,Q, + $ ifac,g_tmp,pi, ! used by AMR + $ tauk,iota_e,D_norm,P_perp,P_tor,delta_eff, ! used by det + $ amr_cell_type,amr_cells,n_amr_cells, ! v2 types + $ Q_store,D_store,n_pts, ! output arrays + $ MAX_PTS,HASH_SZ,HASH_SCALE, ! v1 constants + $ hash_head,hash_next, ! v1 hash + $ MAX_CELLS, ! v2 limit + $ slayer_inputs_type,slayer_outputs_type, + $ deltas_outputs_type, + $ tau_r,dc_tmp,dc_type, + $ sn_str,sm_str + USE delta_mod, ONLY: riccati,riccati_f,riccati_out, $ parflow_flag,PeOhmOnly_flag + USE params_mod + USE layerinputs_mod + USE slayer_netcdf_mod IMPLICIT NONE - + +c --- reconnection regulariser used in psi0 / JxB expressions; +c expose via namelist to make user-configurable. + REAL(r8), PARAMETER :: DELTA_N_PERT = 1.0e-2_r8 + CONTAINS c----------------------------------------------------------------------- c subprogram 1. gpec_slayer. -c run slayer to provide b_crit(ising). -c----------------------------------------------------------------------- -c----------------------------------------------------------------------- -c declarations. +c Single-surface SLAYER driver: compute Delta, torque, and b_crit +c for one rational surface characterised by its plasma profiles. +c +c Steps: +c 1. Derive normalised SLAYER parameters (lu, ds, Q, Q_e, etc.) +c from dimensional inputs. +c 2. Compute baseline Delta via riccati(). +c 3. Scan over a range of Q (rotation) to build a torque +c balance curve, then identify the critical threshold br_th. +c +c TODO: `zeff` and `qval` are INTENT(IN) +c arguments accepted in the interface but never referenced in +c the subroutine body. Remove from signature (breaking API +c change) or add a comment documenting their reserved intent. c----------------------------------------------------------------------- SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, $ omega_i,qval,sval,bt,rs,R0,mu_i,inpr,mms,nns,ascii_flag, $ delta,psi0,jxb,omega_sol,br_th) - REAL(r8),INTENT(IN) :: n_e,t_e,n_i,t_i,omega,omega_e,omega_i, - $ qval,sval,bt,rs,R0,zeff,inpr - INTEGER, INTENT(IN) :: mms,nns,mu_i - LOGICAL, INTENT(IN) :: ascii_flag - COMPLEX(r8),INTENT(OUT) :: delta,psi0 - REAL(r8),INTENT(OUT) :: jxb,omega_sol,br_th - - INTEGER :: i,inum - INTEGER, DIMENSION(1) :: index - - REAL(r8) :: inQ,inQ_e,inQ_i,inpe,inc_beta,inds,intau,inlu - REAL(r8) :: mrs,nrs,rho,b_l,v_a,Qconv,Q0,delta_n_p, - $ lbeta,tau_i,tau_h,tau_r,tau_v - REAL(r8) :: inQ_min,inQ_max,Q_sol - - REAL(r8), DIMENSION(:), ALLOCATABLE :: inQs,iinQs,jxbl,bal - COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: deltal - CHARACTER(3) :: sn,sm +c --- input arguments: dimensional plasma profiles for this surface + REAL(r8),INTENT(IN) :: n_e ! electron density [m^-3] + REAL(r8),INTENT(IN) :: t_e ! electron temperature [eV] + REAL(r8),INTENT(IN) :: n_i ! ion density [m^-3] + REAL(r8),INTENT(IN) :: t_i ! ion temperature [eV] + REAL(r8),INTENT(IN) :: zeff ! (TODO: unused; see header) + REAL(r8),INTENT(IN) :: omega ! plasma rotation frequency + REAL(r8),INTENT(IN) :: omega_e ! electron diamagnetic freq + REAL(r8),INTENT(IN) :: omega_i ! ion diamagnetic freq + REAL(r8),INTENT(IN) :: qval ! (TODO: unused; see header) + REAL(r8),INTENT(IN) :: sval ! magnetic shear + REAL(r8),INTENT(IN) :: bt ! toroidal field [T] + REAL(r8),INTENT(IN) :: rs ! minor radius of surface [m] + REAL(r8),INTENT(IN) :: R0 ! major radius [m] + REAL(r8),INTENT(IN) :: inpr ! Prandtl number + INTEGER, INTENT(IN) :: mms ! poloidal mode number + INTEGER, INTENT(IN) :: nns ! toroidal mode number + INTEGER, INTENT(IN) :: mu_i ! ion mass number (AMU) + LOGICAL, INTENT(IN) :: ascii_flag ! write ASCII torque-balance file +c --- output arguments + COMPLEX(r8),INTENT(OUT) :: delta ! complex Delta + COMPLEX(r8),INTENT(OUT) :: psi0 ! reconnected flux (a.u.) + REAL(r8),INTENT(OUT) :: jxb ! JxB torque (a.u.) + REAL(r8),INTENT(OUT) :: omega_sol ! rotation at torque-balance + REAL(r8),INTENT(OUT) :: br_th ! critical radial field threshold + +c --- loop / scan control + INTEGER :: i ! loop index + INTEGER :: inum ! number of scan points + INTEGER, DIMENSION(1) :: max_idx ! index of max torque balance +c --- local copies of normalised parameters for riccati() + REAL(r8) :: inQ,inQ_e,inQ_i ! normalised frequencies + REAL(r8) :: inpe ! electron Prandtl (set to 0) + REAL(r8) :: inc_beta,inds ! c_beta, ds copies + REAL(r8) :: intau ! tau copy +c --- derived dimensional quantities + REAL(r8) :: mrs,nrs ! real-valued mode numbers + REAL(r8) :: rho ! mass density [kg/m^3] + REAL(r8) :: b_l ! characteristic magnetic field + REAL(r8) :: Qconv ! Q normalisation factor + REAL(r8) :: Q0 ! initial Q (before scan) + REAL(r8) :: delta_n_p ! Delta_n perturbation + REAL(r8) :: lbeta ! local beta + REAL(r8) :: tau_i ! ion collision time + REAL(r8) :: tau_h ! Alfven time across surface + REAL(r8) :: tau_v ! viscous diffusion time +c --- scan workspace + REAL(r8) :: inQ_min,inQ_max ! scan bounds in Q + REAL(r8) :: Q_sol ! Q at torque-balance threshold + REAL(r8), DIMENSION(:), ALLOCATABLE :: inQs ! Q scan array + REAL(r8), DIMENSION(:), ALLOCATABLE :: jxbl ! JxB scan array + REAL(r8), DIMENSION(:), ALLOCATABLE :: bal ! balance scan array + COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: deltal ! delta scan array + CHARACTER(3) :: l_sn,l_sm ! local mode number strings +c --- set delta_mod control flags for this single-surface call parflow_flag=.FALSE. PeOhmOnly_flag=.TRUE. riccati_out=.FALSE. mrs = real(mms,4) nrs = real(nns,4) - - ! String representations of the m and n mode numbers +c----------------------------------------------------------------------- +c build string representations of m and n for file names. +c----------------------------------------------------------------------- IF (nns<10) THEN - WRITE(UNIT=sn,FMT='(I1)') nns - sn=ADJUSTL(sn) + WRITE(UNIT=l_sn,FMT='(I1)') nns + l_sn=ADJUSTL(l_sn) ELSE - WRITE(UNIT=sn,FMT='(I2)') nns + WRITE(UNIT=l_sn,FMT='(I2)') nns ENDIF IF (mms<10) THEN - WRITE(UNIT=sm,FMT='(I1)') mms - sm=ADJUSTL(sm) + WRITE(UNIT=l_sm,FMT='(I1)') mms + l_sm=ADJUSTL(l_sm) ELSEIF (mms<100) THEN - WRITE(UNIT=sm,FMT='(I2)') mms - sm=ADJUSTL(sm) + WRITE(UNIT=l_sm,FMT='(I2)') mms + l_sm=ADJUSTL(l_sm) ELSE - WRITE(UNIT=sm,FMT='(I3)') mms + WRITE(UNIT=l_sm,FMT='(I3)') mms ENDIF - inpe=0.0 ! Waybright added this + inpe=0.0 ! electron Prandtl not used here - tau= t_i/t_e ! ratio of ion to electron temperature - tau_i = 6.6e17*mu_i**0.5*(t_i/1e3)**1.5/(n_e*lnLamb) ! ion colls. - eta= 1.65e-9*lnLamb/(t_e/1e3)**1.5 ! spitzer resistivity (wesson) +c----------------------------------------------------------------------- +c derive normalised SLAYER parameters from dimensional inputs. +c tau, eta, rho, b_l, v_a are intermediate dimensional quantities; +c lu (Lundquist), ds, Qconv, Q, Q_e, Q_i, c_beta are the +c normalised parameters used by riccati(). +c----------------------------------------------------------------------- + tau= t_i/t_e ! ion-to-electron temp ratio + tau_i = 6.6e17*mu_i**0.5*(t_i/1e3)**1.5/(n_e*lnLamb) ! ion coll. time + eta= 1.65e-9*lnLamb/(t_e/1e3)**1.5 ! Spitzer resistivity (Wesson) rho=(mu_i*m_p)*n_e ! mass density - b_l=(nrs/mrs)*rs*sval*bt/R0 ! characteristic magnetic field - v_a=b_l/(mu0*rho)**0.5 ! alfven velocity - rho_s=1.02e-4*(mu_i*t_e)**0.5/bt ! ion Lamour by elec. Temp. + b_l=(nrs/mrs)*nrs*sval*bt/R0 ! characteristic magnetic field + rho_s=1.02e-4*(mu_i*t_e)**0.5/bt ! ion Larmor radius at T_e - tau_h=R0*(mu0*rho)**0.5/(nns*sval*bt) ! alfven time across surface - tau_r=mu0*rs**2.0/eta ! resistive time scale - tau_v=tau_r/inpr ! rho*rs**2.0/visc ! viscous time scale + tau_h=R0*(mu0*rho)**0.5/(nns*sval*bt) ! Alfven transit time + tau_r=mu0*rs**2.0/eta ! resistive diffusion time + tau_v=tau_r/inpr ! viscous diffusion time + visc= rho*rs**2.0/tau_v ! back-calculated viscosity - ! this one must be anomalous. calculated back from pr. - visc= rho*rs**2.0/tau_v + lu=tau_r/tau_h ! Lundquist number + Qconv=lu**(1.0/3.0)*tau_h ! Q normalisation factor - lu=tau_r/tau_h ! Lundquist number - - Qconv=lu**(1.0/3.0)*tau_h ! conversion to Qs based on Cole - - ! note Q depends on Qconv even if omega is fixed. +c --- normalised frequencies Q=Qconv*omega Q_e=-Qconv*omega_e Q_i=-Qconv*omega_i - ! This is the most critical parameter - ds=lu**(1.0/3.0)*rho_s/rs ! conversion based on Cole. + ds=lu**(1.0/3.0)*rho_s/rs ! normalised ion sound radius lbeta=(5.0/3.0)*mu0*n_e*chag*(t_e+t_i)/bt**2.0 - c_beta=(lbeta/(1.0+lbeta))**0.5 + c_beta=(lbeta/(1.0+lbeta))**0.5 ! compressibility parameter - delta_n=lu**(1.0/3.0)/rs ! norm factor for delta primes + delta_n=lu**(1.0/3.0)/rs ! Delta normalisation factor +c --- copy normalised values into local variables for riccati() call inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -109,14 +191,24 @@ SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, intau=tau Q0=Q c----------------------------------------------------------------------- -c calculate basic delta, torque, balance, error fields. +c compute baseline Delta, reconnected flux, and JxB torque. +c delta_n_p uses the module-level DELTA_N_PERT constant; promote +c that constant to a namelist input to make it user-configurable. c----------------------------------------------------------------------- - delta_n_p=1e-2 + delta_n_p = DELTA_N_PERT delta=riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) - psi0=1.0/ABS(delta+delta_n_p) ! a.u. - jxb=-AIMAG(1.0/(delta+delta_n_p)) ! a.u. + psi0=1.0/ABS(delta+delta_n_p) ! reconnected flux (a.u.) + jxb=-AIMAG(1.0/(delta+delta_n_p)) ! JxB torque (a.u.) c----------------------------------------------------------------------- -c find solutions based on simple torque balance. +c torque-balance scan: sweep Q over [inQ_min, inQ_max] and +c compute delta(Q), JxB(Q), and the balance parameter. +c The threshold br_th is sqrt(max(bal)/lu * s^2/2). +c +c TODO: physics bounds from Q0/Q_e (computed +c in the IF/ELSE above) are immediately overridden by the two +c fixed assignments below. To use physics bounds, remove the +c inQ_max=10 / inQ_min=-10 lines. To keep fixed bounds, remove +c the dead IF/ELSE block above. c----------------------------------------------------------------------- IF (Q0>inQ_e) THEN inQ_max=2.0*Q0 @@ -130,11 +222,10 @@ SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, ENDIF ENDIF - ! Scan of rotation - inQ_max=10.0 + inQ_max=10.0 ! TODO: fixed bound; see header inQ_min=-10.0 inum=200 - ALLOCATE(inQs(0:inum),deltal(0:inum),jxbl(0:inum),bal(0:inum)) + ALLOCATE(inQs(0:inum),deltal(0:inum),jxbl(0:inum),bal(0:inum)) DO i=0,inum inQs(i)=inQ_min+(REAL(i)/inum)*(inQ_max-inQ_min) deltal(i)=riccati(inQs(i),inQ_e,inQ_i, @@ -143,10 +234,10 @@ SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, bal(i)=2.0*inpr*(Q0-inQs(i))/jxbl(i) ENDDO - ! Write torque balance curves to file for diagnostic purposes +c --- optionally write torque balance curve to ASCII file IF(ascii_flag)THEN OPEN(UNIT=out_unit,FILE="gpec_slayer_torque_balance_m"// - $ TRIM(sm)//"_n"//TRIM(sn)//".out", + $ TRIM(l_sm)//"_n"//TRIM(l_sn)//".OUT", $ STATUS="UNKNOWN") WRITE(out_unit,'(1x,5(a17))') "inQ","RE(delta)", $ "IM(delta)","jxb","bal" @@ -157,16 +248,13 @@ SUBROUTINE gpec_slayer(n_e,t_e,n_i,t_i,zeff,omega,omega_e, CLOSE(out_unit) ENDIF - ! Identify the threshold from the maximum of the balance parameter - index=MAXLOC(bal) - Q_sol=inQs(index(1)) - omega_sol=inQs(index(1))/Qconv +c --- identify the critical threshold from the maximum balance value + max_idx=MAXLOC(bal) + Q_sol=inQs(max_idx(1)) + omega_sol=inQs(max_idx(1))/Qconv br_th=sqrt(MAXVAL(bal)/lu*(sval**2.0/2.0)) DEALLOCATE(inQs,deltal,jxbl,bal) RETURN END SUBROUTINE gpec_slayer - - END MODULE gslayer_mod - - + END MODULE gslayer_mod \ No newline at end of file diff --git a/slayer/layerinputs.f b/slayer/layerinputs.f new file mode 100644 index 00000000..ffedd819 --- /dev/null +++ b/slayer/layerinputs.f @@ -0,0 +1,581 @@ +c======================================================================= +c MODULE layerinputs_mod +c +c Reads STRIDE NetCDF output and constructs the per-surface +c input arrays required by the SLAYER layer-physics solver. +c +c Subprograms contained: +c 1. read_stride_netcdf_diagonal -- read STRIDE NetCDF, extract +c Deltaprime diagonal, geometry, and equilibrium scalars. +c 2. issurfint -- surface integral by simple +c quadrature (adapted from EQUIL). +c 3. build_inputs -- master driver: reads kinetic +c profiles, evaluates params(), and populates +c slayer_inputs_type for every rational surface. +c======================================================================= + MODULE layerinputs_mod + + USE inputs, ONLY : read_kin, read_equil, kin, chi1 + USE spline_mod, ONLY : spline_alloc, spline_eval, spline_type, + $ spline_dealloc, spline_int, spline_fit + USE sglobal_mod ! SLAYER global scalars, types, constants + USE params_mod ! params() -- compute derived layer params + USE netcdf ! NetCDF Fortran bindings + USE equil_mod, ONLY : equil_read, rzphi, twopi, ro, zo, sq + USE bicube_mod, ONLY : bicube_eval_external, bicube_type + USE slayer_netcdf_mod ! sl_check(), SLAYER NetCDF output + + IMPLICIT NONE + + CONTAINS +c----------------------------------------------------------------------- +c subprogram 1. read_stride_netcdf_diagonal. +c Read the STRIDE NetCDF file and extract the Deltaprime matrix, +c rational-surface geometry (q, psi_n, shear, dgeo, dr), and +c equilibrium scalars (R0, B_t0, psi0, m_psi, n, resm). +c All global-attribute reads use NF90_GLOBAL explicitly; +c variable IDs are obtained via nf90_inq_varid before use. +c----------------------------------------------------------------------- + SUBROUTINE read_stride_netcdf_diagonal(ncfile,msing,dp_mat, + $ Re_dp_diagonal,Im_dp_diagonal,q_rational,psi_n_rational,dgeo, + $ shear,r_o,my_bt0,my_psio,dr_vals,mpsi,nn,resm) + +c --- arguments (all INTENT(OUT) except ncfile) + CHARACTER(512), INTENT(IN) :: ncfile ! path to STRIDE NetCDF + INTEGER, INTENT(OUT) :: msing ! number of singular surfaces + REAL(r8), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) :: dp_mat + ! Deltaprime matrix (msing x msing x 2) + REAL(r8), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: + $ Re_dp_diagonal, ! Re diag(Deltaprime) + $ Im_dp_diagonal ! Im diag(Deltaprime) + REAL(r8), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: + $ q_rational, ! rational-surface q values + $ psi_n_rational, ! normalised psi at each surface + $ shear, ! magnetic shear + $ dgeo ! geometric delta (Shafranov shift) + REAL(r8), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: + $ r_o, ! major radius R0 + $ my_bt0, ! toroidal field B_t0 + $ my_psio, ! poloidal flux psi_0 + $ mpsi, ! poloidal mode number array + $ dr_vals ! radial width dr at each surface + INTEGER, DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: + $ nn, ! toroidal mode number(s) + $ resm ! resonant poloidal mode numbers + +c --- locals + REAL(r8), DIMENSION(:), ALLOCATABLE :: msing_arr ! temp for reading msing attribute + INTEGER(kind=nf90_int) :: ncid, stat ! NetCDF file id / return status + INTEGER(kind=nf90_int) :: r_dim_id, r_dim ! dimension id / length (unused) + INTEGER(kind=nf90_int) :: dp_id, qr_id, pr_id ! variable ids + INTEGER(kind=nf90_int) :: dgeo_id, shear_id ! variable ids + INTEGER(kind=nf90_int) :: resm_id, drr_id ! variable ids + INTEGER(kind=nf90_int), DIMENSION(1) :: start, count ! NetCDF hyperslab + INTEGER :: i ! loop index + INTEGER :: bt0_len, ro_len, psio_len ! attribute lengths + INTEGER :: mpsi_len, msing_len, nn_len, dr_len ! attribute lengths + +c----------------------------------------------------------------------- +c open the STRIDE NetCDF file and read dimension / attribute data. +c----------------------------------------------------------------------- + stat = nf90_open(path=ncfile, mode=NF90_NOWRITE, + $ ncid=ncid) + CALL sl_check(stat) + +c --- read msing (number of singular surfaces) from global attribute + stat = nf90_inquire_attribute(ncid, NF90_GLOBAL, 'msing', + $ len = msing_len) + CALL sl_check(stat) + ALLOCATE(msing_arr(msing_len)) + stat = nf90_get_att(ncid, NF90_GLOBAL, 'msing', msing_arr) + CALL sl_check(stat) + msing = INT(msing_arr(1)) + +c --- allocate output arrays sized by msing + ALLOCATE(Re_dp_diagonal(msing), q_rational(msing), + $ psi_n_rational(msing), shear(msing), dgeo(msing), + $ resm(msing), Im_dp_diagonal(msing), dr_vals(msing)) + ALLOCATE(dp_mat(msing, msing, 2)) + +c --- read lengths of scalar / small-array global attributes + stat = nf90_inquire_attribute(ncid, NF90_GLOBAL, 'ro', + $ len=ro_len) + CALL sl_check(stat) + stat = nf90_inquire_attribute(ncid, NF90_GLOBAL, 'bt0', + $ len=bt0_len) + CALL sl_check(stat) + stat = nf90_inquire_attribute(ncid, NF90_GLOBAL, 'psio', + $ len=psio_len) + CALL sl_check(stat) + stat = nf90_inquire_attribute(ncid, NF90_GLOBAL, 'mpsi', + $ len=mpsi_len) + CALL sl_check(stat) + stat = nf90_inquire_attribute(ncid, NF90_GLOBAL, 'n', + $ len=nn_len) + CALL sl_check(stat) + + ALLOCATE(my_bt0(INT(bt0_len)), r_o(INT(ro_len)), + $ my_psio(INT(psio_len)), + $ mpsi(INT(mpsi_len)), nn(INT(nn_len))) + +c --- obtain NetCDF variable IDs + stat = nf90_inq_varid(ncid, 'Delta_prime', dp_id) + CALL sl_check(stat) + stat = nf90_inq_varid(ncid, 'q_rational', qr_id) + CALL sl_check(stat) + stat = nf90_inq_varid(ncid, 'psi_n_rational', pr_id) + CALL sl_check(stat) + stat = nf90_inq_varid(ncid, 'Delta_geo', dgeo_id) + CALL sl_check(stat) + stat = nf90_inq_varid(ncid, 'shear', shear_id) + CALL sl_check(stat) + stat = nf90_inq_varid(ncid, 'resm', resm_id) + CALL sl_check(stat) + stat = nf90_inq_varid(ncid, 'dr_rational', drr_id) + CALL sl_check(stat) + +c --- read global attributes (equilibrium scalars) + stat = nf90_get_att(ncid, NF90_GLOBAL, 'ro', r_o) + CALL sl_check(stat) + stat = nf90_get_att(ncid, NF90_GLOBAL, 'bt0', my_bt0) + CALL sl_check(stat) + stat = nf90_get_att(ncid, NF90_GLOBAL, 'psio', my_psio) + CALL sl_check(stat) + stat = nf90_get_att(ncid, NF90_GLOBAL, 'mpsi', mpsi) + CALL sl_check(stat) + stat = nf90_get_att(ncid, NF90_GLOBAL, 'n', nn) + CALL sl_check(stat) + +c --- read variable data: Deltaprime matrix and 1-D surface arrays + stat = nf90_get_var(ncid, dp_id, dp_mat, start=(/1,1,1/)) + CALL sl_check(stat) + stat = nf90_get_var(ncid, qr_id, q_rational) + CALL sl_check(stat) + stat = nf90_get_var(ncid, pr_id, psi_n_rational) + CALL sl_check(stat) + stat = nf90_get_var(ncid, dgeo_id, dgeo) + CALL sl_check(stat) + stat = nf90_get_var(ncid, shear_id, shear) + CALL sl_check(stat) + stat = nf90_get_var(ncid, resm_id, resm) + CALL sl_check(stat) + stat = nf90_get_var(ncid, drr_id, dr_vals) + CALL sl_check(stat) + +c --- extract diagonal of the complex Deltaprime matrix + DO i = 1, msing + Re_dp_diagonal(i) = dp_mat(i, i, 1) + Im_dp_diagonal(i) = dp_mat(i, i, 2) + END DO + +c --- close file + stat = nf90_close(ncid) + CALL sl_check(stat) + + END SUBROUTINE read_stride_netcdf_diagonal +c----------------------------------------------------------------------- +c subprogram 2. issurfint. +c Surface integral by simple quadrature, adapted from EQUIL. +c Computes int f(theta) * W(theta) d(theta) where W depends +c on the weight flag `wegt`: +c 0 = jac * |grad psi| +c 1 = R * jac * |grad psi| (R-weighted) +c 2 = jac * |grad psi| / R (1/R-weighted) +c 3 = a * jac * |grad psi| (minor-radius weighted) +c If ave==1, the result is divided by the unweighted surface area. +c +c Geometry is cached via first/fsave/psave to avoid recomputing +c Jacobians when the same surface is queried repeatedly. +c + +c----------------------------------------------------------------------- + FUNCTION issurfint(func,fs,inpsi,wegt,ave, + $ fsave,psave,jacs,delpsi,inr,ina,first) +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- +c --- arguments + INTEGER, INTENT(IN) :: fs ! number of poloidal segments + INTEGER, INTENT(IN) :: wegt ! weight flag (0..3) + INTEGER, INTENT(IN) :: ave ! 1 = return surface average + REAL(r8), INTENT(IN) :: inpsi ! normalised psi of surface + REAL(r8), DIMENSION(0:fs), INTENT(IN) :: func ! integrand array + + LOGICAL, INTENT(INOUT) :: first ! .TRUE. on first call for caching + INTEGER, INTENT(INOUT) :: fsave ! cached fs + REAL(r8), INTENT(INOUT) :: psave ! cached psi + REAL(r8), DIMENSION(0:), INTENT(INOUT) :: jacs ! Jacobian cache + REAL(r8), DIMENSION(0:), INTENT(INOUT) :: delpsi ! |grad psi| cache + REAL(r8), DIMENSION(0:), INTENT(INOUT) :: inr ! R(theta) cache + REAL(r8), DIMENSION(0:), INTENT(INOUT) :: ina ! a(theta) cache +c --- return value + REAL(r8) :: issurfint +c --- locals + INTEGER :: itheta ! poloidal loop index + INTEGER :: ix, iy ! bicube grid hints + REAL(r8) :: rfac, ineta, injac, inarea ! geometry intermediates + REAL(r8), DIMENSION(1,2) :: w ! gradient components + REAL(r8), DIMENSION(0:fs) :: thetas ! normalised theta grid + REAL(r8), DIMENSION(4) :: rzphi_f, rzphi_fx, rzphi_fy + ! bicube_eval_external outputs +c----------------------------------------------------------------------- +c compute / cache geometry if surface changed. +c [bicube_eval_external]: external bicubic interpolation from EQUIL. +c----------------------------------------------------------------------- + issurfint = 0 + inarea = 0 + ix = 0 + iy = 0 + + IF (first .OR. inpsi /= psave .OR. fs /= fsave) THEN + psave = inpsi + fsave = fs + first = .FALSE. + DO itheta = 0, fs + thetas(itheta) = REAL(itheta, r8) / REAL(fs, r8) + ENDDO + DO itheta = 0, fs-1 + CALL bicube_eval_external(rzphi, inpsi, thetas(itheta), 1, + $ ix, iy, rzphi_f, rzphi_fx, rzphi_fy) + rfac = SQRT(rzphi_f(1)) + ineta = twopi * (thetas(itheta) + rzphi_f(2)) + ina(itheta) = rfac + inr(itheta) = ro + rfac * COS(ineta) + injac = rzphi_f(4) + jacs(itheta) = injac +c gradient magnitude: |grad psi| from metric components + w(1,1) = (1 + rzphi_fy(2)) * twopi**2 + $ * rfac * inr(itheta) / injac + w(1,2) = -rzphi_fy(1) * pi * inr(itheta) + $ / (rfac * injac) + delpsi(itheta) = SQRT(w(1,1)**2 + w(1,2)**2) + ENDDO + ENDIF +c----------------------------------------------------------------------- +c perform weighted surface integral based on wegt flag. +c----------------------------------------------------------------------- + IF (wegt == 0) THEN + DO itheta = 0, fs-1 + issurfint = issurfint + $ + jacs(itheta)*delpsi(itheta)*func(itheta)/fs + ENDDO + ELSE IF (wegt == 1) THEN + DO itheta = 0, fs-1 + issurfint = issurfint + $ + inr(itheta)*jacs(itheta)*delpsi(itheta)* + $ func(itheta)/fs + ENDDO + ELSE IF (wegt == 2) THEN + DO itheta = 0, fs-1 + issurfint = issurfint + $ + jacs(itheta)*delpsi(itheta)* + $ func(itheta)/inr(itheta)/fs + ENDDO + ELSE IF (wegt == 3) THEN + DO itheta = 0, fs-1 + issurfint = issurfint + $ + ina(itheta)*jacs(itheta)*delpsi(itheta)* + $ func(itheta)/fs + ENDDO + ELSE + STOP 'ERROR: issurfint wegt must be in [0,1,2,3]' + ENDIF +c----------------------------------------------------------------------- +c optionally normalise by unweighted surface area. +c----------------------------------------------------------------------- + IF (ave == 1) THEN + DO itheta = 0, fs-1 + inarea = inarea + jacs(itheta)*delpsi(itheta)/fs + ENDDO + issurfint = issurfint / inarea + ENDIF +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END FUNCTION issurfint +c----------------------------------------------------------------------- +c subprogram 3. build_inputs. +c Master driver: reads the STRIDE NetCDF file, sets up kinetic +c profiles from the input file, reads the equilibrium, and then +c loops over every singular surface to compute derived layer +c parameters via params() and populate slayer_inputs_type. +c + +c----------------------------------------------------------------------- + SUBROUTINE build_inputs(infile,ncfile,sl_in) +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- +c --- arguments + CHARACTER(512), INTENT(IN) :: infile ! kinetic input file path + CHARACTER(512), INTENT(IN) :: ncfile ! STRIDE NetCDF file path + TYPE(slayer_inputs_type), INTENT(INOUT) :: sl_in +c --- surface-loop control + LOGICAL :: firstsurf ! first-call flag for issurfint + REAL(r8) :: respsi ! normalised psi at current surface + INTEGER :: ising ! loop index over singular surfaces +c --- kinetic profile parameters (for read_kin) + INTEGER :: zi, zimp, mi, mimp ! charge/mass species ids + REAL(r8) :: nfac, tfac, wefac, wpfac ! profile scale factors + REAL(r8) :: e ! elementary charge [C] +c --- mode number workspace + INTEGER :: mms, nns ! poloidal / toroidal mode nums + INTEGER :: mpsi ! poloidal-flux index from attr +c --- local plasma quantities at current surface + REAL(r8) :: n_e, t_e, n_i, t_i ! densities [m^-3], temperatures [eV] + REAL(r8) :: omega, omega_e, omega_i ! toroidal & diamagnetic freqs [rad/s] + REAL(r8) :: my_qval, my_sval ! safety factor, magnetic shear + REAL(r8) :: my_bt, my_rs, R_0 ! toroidal field, minor radius, major radius + REAL(r8) :: zeff ! effective charge from kin%f(9) + REAL(r8) :: dgeo_val ! geometric delta (Shafranov) + REAL(r8) :: mu_i ! ion mass ratio to proton + REAL(r8) :: dr_val ! radial width dr at surface + REAL(r8) :: l_n, l_t ! density / temperature gradient lengths + REAL(r8) :: gammafac ! growth-rate conversion factor + REAL(r8), DIMENSION(3) :: chi_s ! chi_perp, chi_tor, kappa +c --- STRIDE data (from read_stride_netcdf_diagonal) + REAL(r8), DIMENSION(:,:,:), ALLOCATABLE :: dp_mat + REAL(r8), DIMENSION(:), ALLOCATABLE :: Re_dp_diagonal, + $ Im_dp_diagonal, q_rational, psi_n_rational, + $ shear, dgeo, r_o, my_bt0, my_psio, mpsi_arr, + $ dr_vals, dr_arr +c --- local per-surface kinetic arrays (for future NetCDF diagnostic output) + REAL(r8), DIMENSION(:), ALLOCATABLE :: ne_arr, te_arr, + $ ni_arr, ti_arr, zeff_arr, bt_arr, rs_arr, R0_arr, + $ mu_i_arr, omegas_e_arr, omegas_i_arr + INTEGER, DIMENSION(:), ALLOCATABLE :: nn, resm, nns_arr +c --- surface-integral workspace + INTEGER :: msing, i, mthsurf ! surface count, loop idx, theta pts + REAL(r8), DIMENSION(0:512) :: unitfun ! unit function for area integral + INTEGER :: fsave ! cached fs for issurfint + REAL(r8) :: psave ! cached psi for issurfint + REAL(r8), DIMENSION(:), ALLOCATABLE :: jacs, delpsi, rsurf, asurf + REAL(r8) :: a_surf ! flux-surface-averaged minor radius +c --- Fitzpatrick (r-based) shear workspace + REAL(r8) :: a_surf_p, a_surf_m ! a_surf at psiN +/- h + REAL(r8) :: da_dpsiN ! da_surf/dpsiN (Jacobian) + REAL(r8) :: dpsi_h ! finite-diff step for Jacobian + REAL(r8) :: s_fitz ! Fitzpatrick shear r*dq/dr/q +c----------------------------------------------------------------------- +c read STRIDE NetCDF: Deltaprime matrix, geometry, equilibrium scalars. +c----------------------------------------------------------------------- + CALL read_stride_netcdf_diagonal(ncfile,msing,dp_mat, + $ Re_dp_diagonal,Im_dp_diagonal,q_rational, + $ psi_n_rational,dgeo,shear,r_o,my_bt0,my_psio,dr_vals, + $ mpsi_arr,nn,resm) + + mpsi = INT(mpsi_arr(1)) + mthsurf = 512 ! poloidal segments for surface integral + +c --- allocate slayer_inputs_type arrays (one entry per surface) + ALLOCATE(sl_in%qval_arr(msing), sl_in%omegas_arr(msing), + $ sl_in%omegas_e_arr(msing), sl_in%dp_matrix(msing,msing), + $ sl_in%omegas_i_arr(msing), + $ sl_in%Q_e_arr(msing), sl_in%Q_i_arr(msing), + $ sl_in%psi_n_arr(msing), + $ sl_in%Re_dp_arr(msing), sl_in%Im_dp_arr(msing), + $ sl_in%d_crit_arr(msing), sl_in%P_tor_arr(msing), + $ sl_in%P_perp_arr(msing), sl_in%tau_arr(msing), + $ sl_in%D_norm_arr(msing), + $ sl_in%d_beta_arr(msing), sl_in%gammafac_arr(msing), + $ sl_in%c_beta_arr(msing), sl_in%lu_arr(msing), + $ sl_in%Qconv_arr(msing)) + +c --- allocate local kinetic arrays (for future NetCDF diagnostic output) + ALLOCATE(ne_arr(msing), te_arr(msing), ni_arr(msing), + $ ti_arr(msing), zeff_arr(msing), bt_arr(msing), + $ rs_arr(msing), R0_arr(msing), mu_i_arr(msing), + $ nns_arr(msing), dr_arr(msing), + $ omegas_e_arr(msing), omegas_i_arr(msing)) + + ALLOCATE(jacs(0:mthsurf), delpsi(0:mthsurf), + $ rsurf(0:mthsurf), asurf(0:mthsurf)) +c----------------------------------------------------------------------- +c set up kinetic profiles and equilibrium. +c [read_kin]: external, reads kinetic input file into kin spline. +c [equil_read]: external, reads equilibrium into EQUIL module. +c----------------------------------------------------------------------- + zi = 1 ! main-ion charge + zimp = 6 ! impurity charge (carbon) + mi = 2 ! main-ion mass (deuterium) + mimp = 12 ! impurity mass + nfac = 1.0 + tfac = 1.0 + wefac = 1.0 + wpfac = 1.0 + e = 1.6021917e-19 ! elementary charge [C] + chi1 = twopi * my_psio(1) ! total poloidal flux (module-level) + + CALL read_kin(infile,zi,zimp,mi,mimp,nfac, + $ tfac,wefac,wpfac,.false.) + + CALL equil_read(out_unit) + +c store full complex Deltaprime matrix in sl_in + sl_in%dp_matrix(:,:) = CMPLX(dp_mat(:,:,1), dp_mat(:,:,2)) + + dpsi_h = 0.002 ! finite-difference step for Jacobian + +c----------------------------------------------------------------------- +c loop over singular surfaces: evaluate kinetic/equilibrium +c quantities via spline interpolation, compute derived layer +c parameters via params(), and populate sl_in arrays. +c----------------------------------------------------------------------- + DO ising = 1, msing + + respsi = psi_n_rational(ising) ! normalised psi + firstsurf = .TRUE. + unitfun = 1 + +c compute flux-surface-averaged minor radius + a_surf = issurfint(unitfun,mthsurf,respsi,3,1, + $ fsave,psave,jacs,delpsi,rsurf,asurf,firstsurf) + +c compute Jacobian da_surf/dpsiN by central difference + a_surf_p = issurfint(unitfun,mthsurf, + $ MIN(respsi+dpsi_h, REAL(1.0,r8)),3,1, + $ fsave,psave,jacs,delpsi,rsurf,asurf,firstsurf) + a_surf_m = issurfint(unitfun,mthsurf, + $ MAX(respsi-dpsi_h, REAL(0.001,r8)),3,1, + $ fsave,psave,jacs,delpsi,rsurf,asurf,firstsurf) + da_dpsiN = (a_surf_p - a_surf_m) + $ / (MIN(respsi+dpsi_h, REAL(1.0,r8)) + $ - MAX(respsi-dpsi_h, REAL(0.001,r8))) + +c----------------------------------------------------------------------- +c evaluate kinetic splines at this surface. +c [spline_eval]: external, evaluates kin spline at respsi. +c kin%f(1..5) = n_i, n_e, t_i, t_e, omega (SI units) +c kin%f1(1..5) = d/d(psi_n) of the above +c----------------------------------------------------------------------- + CALL spline_eval(kin,respsi,1) + +c diamagnetic frequencies (rad/s) from GPEC kinetic splines. +c These compute the ELECTRON diamagnetic frequency directly. + omega_e = twopi*kin%f(4)*kin%f1(2)/(e*chi1*kin%f(2)) + $ +twopi*kin%f1(4)/(e*chi1) + omega_i = -twopi*kin%f(3)*kin%f1(1)/(e*zi*chi1*kin%f(1)) + $ -twopi*kin%f1(3)/(e*zi*chi1) + sl_in%omegas_e_arr(ising) = omega_e + sl_in%omegas_i_arr(ising) = omega_i + +c extract local plasma quantities from spline + n_e = kin%f(2) + t_e = kin%f(4) / e ! convert J -> eV + n_i = kin%f(1) + t_i = kin%f(3) / e + +c Z_eff: kinetic spline value may be incorrect if ni=ne in gpeckf +c (quasi-neutrality assumption makes Zeff=1). Override to 2.0 +c for deuterium plasma with carbon impurities (matching TJ). +c TODO: fix gpeckf generation to include proper ni for Zeff, +c or read Zeff from a namelist parameter. + zeff = 2.0 + + omega = kin%f(5) + my_qval = q_rational(ising) + my_sval = shear(ising) + dgeo_val = dgeo(ising) + my_bt = my_bt0(1) + my_rs = a_surf + R_0 = r_o(1) + mu_i = 2.0 ! deuterium + dr_val = dr_vals(ising) + +c convert STRIDE shear (psiN-based) to Fitzpatrick shear (r-based). +c STRIDE's shear(ising) = psi_N * d(ln q)/d(psi_N) (psi_N-based). +c Fitzpatrick's rfitzp formula needs s = (r/q) * dq/dr (r-based). +c Chain-rule: s_Fitz = s_psiN * r_s / (psi_N * da/dpsi_N). +c Typical magnitudes: s_psiN ~ 0.3-0.5 mid-plasma, while s_Fitz +c ~ 0.7-1.5 -- r-based shear is usually LARGER than psi_N-based +c because (r/psi_N) > (dr/dpsi_N) inside the mid-radius region. + s_fitz = my_sval * my_rs / (respsi * da_dpsiN) + WRITE(*,'(A,I2,A,F7.4,A,F6.3,A,F8.4,A,F8.4,A,F7.3,A,F7.3)') + $ ' [layerinputs] q=', resm(ising), + $ ' psiN=', respsi, + $ ' r_s=', my_rs, + $ ' da/dpsiN=', da_dpsiN, + $ ' s_psiN=', my_sval, + $ ' -> s_Fitz=', s_fitz + +c transport coefficients from caller-provided arrays. +c guard: arrays may be smaller than msing (e.g. from +c fixed-size namelist); reuse last element if exceeded. + i = MIN(ising, SIZE(sl_in%chi_p_arr)) + chi_s(1) = sl_in%chi_p_arr(i) + i = MIN(ising, SIZE(sl_in%chi_t_arr)) + chi_s(2) = sl_in%chi_t_arr(i) + i = MIN(ising, SIZE(sl_in%kappa_arr)) + chi_s(3) = sl_in%kappa_arr(i) + +c store local kinetic arrays (for future NetCDF diagnostic output) + ne_arr(ising) = n_e + te_arr(ising) = t_e + ni_arr(ising) = n_i + ti_arr(ising) = t_i + zeff_arr(ising) = zeff + bt_arr(ising) = my_bt + rs_arr(ising) = my_rs + R0_arr(ising) = R_0 + mu_i_arr(ising) = mu_i + + mms = resm(ising) + nns = nn(1) + nns_arr(ising) = nn(1) + nr = nn(1) ! module-level toroidal mode number + + l_n = 0.0 + l_t = 0.0 + +c----------------------------------------------------------------------- +c compute derived layer parameters using Fitzpatrick (r-based) +c shear. params() sees s_fitz as `sval`, so lu, tauk, D_norm, +c dc_tmp, etc. are all Fitzpatrick-consistent. Gradient lengths +c are zero here; params() skips its own omega_e/omega_i and we +c set Q_e/Q_i below from the spline-derived frequencies. +c----------------------------------------------------------------------- + CALL params(n_e,t_e,t_i,omega,chi_s,dr_val,dgeo_val, + $ l_n,l_t,my_qval,s_fitz,my_bt,my_rs,R_0,mu_i, + $ zeff,.false.) + +c growth-rate conversion factor: Deltaprime -> gamma + gammafac = (my_rs * Re_dp_diagonal(ising)) / tau_r + +c----------------------------------------------------------------------- +c populate sl_in for this surface from params() globals. +c Q_e/Q_i use spline-derived omega_e/omega_i normalised by +c tauk (= Fitzpatrick S^(1/3) * tau_H, from params()). +c----------------------------------------------------------------------- + sl_in%qval_arr(ising) = INT(my_qval) + sl_in%lu_arr(ising) = lu +c Use the same -tauk prefactor for both so opposite-signed +c spline omega_e, omega_i (electron vs ion diamagnetic +c frequencies) produce opposite-signed Q_e, Q_i. Matches +c params.f convention: Q_e = -Qconv*omega_e, Q_i = -Qconv*omega_i. +c Using +tauk on Q_i cancels the physical sign and drives +c Q_e - Q_i -> 0 for ni=ne, Ti=Te (iota_e blows up). + sl_in%Q_e_arr(ising) = -tauk * omega_e + sl_in%Q_i_arr(ising) = -tauk * omega_i + sl_in%c_beta_arr(ising) = c_beta + sl_in%d_beta_arr(ising) = d_beta + sl_in%D_norm_arr(ising) = D_norm + sl_in%tau_arr(ising) = tau + sl_in%omegas_arr(ising) = omega + sl_in%psi_n_arr(ising) = respsi + sl_in%gammafac_arr(ising)= gammafac + sl_in%Re_dp_arr(ising) = Re_dp_diagonal(ising) + sl_in%Im_dp_arr(ising) = Im_dp_diagonal(ising) + sl_in%d_crit_arr(ising) = dc_tmp + sl_in%P_perp_arr(ising) = P_perp + sl_in%P_tor_arr(ising) = P_tor + sl_in%Qconv_arr(ising) = tauk + ENDDO +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + + END SUBROUTINE build_inputs + + END MODULE layerinputs_mod \ No newline at end of file diff --git a/slayer/params.f b/slayer/params.f index fcf98871..3324165e 100644 --- a/slayer/params.f +++ b/slayer/params.f @@ -1,76 +1,266 @@ +c======================================================================= +c MODULE params_mod +c +c Computes the normalised layer-physics parameters needed by the +c SLAYER dispersion solver from dimensional equilibrium and +c kinetic-profile inputs. +c +c Subprograms contained: +c 1. params -- derive all normalised quantities (Q, ds, c_beta, +c D_norm, P_perp, P_tor, lu, delta_crit, ...) for a +c single rational surface and store them in sglobal_mod. +c======================================================================= MODULE params_mod - USE sglobal_mod + USE sglobal_mod ! SLAYER global scalars, types, constants IMPLICIT NONE CONTAINS c----------------------------------------------------------------------- -c calculate parameters. +c subprogram 1. params. +c Compute all derived layer-physics parameters for a single +c rational surface from dimensional equilibrium / kinetic inputs. +c Results are written to module-level variables in sglobal_mod +c (tau, tau_r, tauk, lu, Q, Q_e, Q_i, ds, c_beta, d_beta, +c D_norm, P_perp, P_tor, delta_n, dc_tmp, eta, visc, rho_s, ...). +c +c Note: `pr` (magnetic Prandtl number) is read from sglobal_mod +c and must be set by the caller (e.g. pr = inpr from namelist). +c If pr == 0, tau_v and visc are set to zero (no viscosity). c----------------------------------------------------------------------- - SUBROUTINE params(n_e,t_e,t_i,omega, + SUBROUTINE params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) - REAL(r8), INTENT(IN) :: n_e,t_e,t_i,omega, - $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff +c --- arguments + REAL(r8), INTENT(IN) :: n_e ! electron density [m^-3] + REAL(r8), INTENT(IN) :: t_e ! electron temperature [eV] + REAL(r8), INTENT(IN) :: t_i ! ion temperature [eV] + REAL(r8), INTENT(IN) :: omega ! toroidal rotation [rad/s] + REAL(r8), INTENT(IN) :: dr_val ! radial width dr at surface + REAL(r8), INTENT(IN) :: dgeo_val ! geometric delta (Shafranov shift factor) + REAL(r8), INTENT(IN) :: l_n ! density gradient length + REAL(r8), INTENT(IN) :: l_t ! temperature gradient length + REAL(r8), INTENT(IN) :: qval ! safety factor + REAL(r8), INTENT(IN) :: sval ! magnetic shear + REAL(r8), INTENT(IN) :: bt ! toroidal field [T] + REAL(r8), INTENT(IN) :: rs ! minor radius [m] + REAL(r8), INTENT(IN) :: R0 ! major radius [m] + REAL(r8), INTENT(IN) :: mu_i ! ion mass ratio to proton + REAL(r8), INTENT(IN) :: zeff ! effective charge + REAL(r8), DIMENSION(3), INTENT(IN) :: chis + ! (1) chi_perp [m^2/s] + ! (2) chi_tor [m^2/s] + ! (3) kappa [m^2/s] + LOGICAL, INTENT(IN) :: params_check ! .TRUE. = print diagnostics - LOGICAL, INTENT(IN) :: params_check +c --- local variables: basic plasma + REAL(r8) :: rho ! mass density [kg/m^3] + REAL(r8) :: b_l ! characteristic magnetic field [T] + REAL(r8) :: v_a ! Alfvén velocity [m/s] + REAL(r8) :: lbeta ! local beta-related quantity +c --- local variables: electron collision time (Braginskii) + REAL(r8) :: tau_ee_num ! numerator of tau_ee formula + REAL(r8) :: tau_ee_denom ! denominator of tau_ee formula + REAL(r8) :: tau_ee ! electron-electron collision time [s] +c --- local variables: parallel conductivity (Spitzer-Härm) + REAL(r8) :: sigma_par_1 ! neoclassical correction factor + REAL(r8) :: sigma_par_2 ! classical conductivity [1/(Ohm*m)] + REAL(r8) :: sigma_par ! parallel conductivity [1/(Ohm*m)] +c --- local variables: timescales + REAL(r8) :: tau_i ! ion collision time [s] + REAL(r8) :: tau_h ! Alfvén transit time [s] + REAL(r8) :: tau_v ! viscous time [s] + REAL(r8) :: tau_tor ! toroidal diffusion time [s] + REAL(r8) :: tau_perp ! perp. diffusion time [s] +c --- local variables: delta_crit iteration + REAL(r8) :: vte ! thermal electron speed [m/s] + REAL(r8) :: chi_par_smfp ! chi_par in short mfp limit + REAL(r8) :: chi_par_lmfp ! chi_par in long mfp limit + REAL(r8) :: chi_par ! effective parallel thermal cond. + REAL(r8) :: Wd ! magnetic island width proxy + REAL(r8) :: Wd_new ! updated Wd for convergence check + INTEGER :: wit ! iteration counter +c --- local variables: unused intermediates + REAL(r8) :: K_val, Csq ! kappa/eta and composite quantity - REAL(r8) :: rho,b_l,v_a,Qconv, - $ lbeta,tau_i,tau_h,tau_r,tau_v +c----------------------------------------------------------------------- +c Coulomb logarithm, basic plasma quantities, and Spitzer +c resistivity. +c----------------------------------------------------------------------- + lnLamb = 24 + 3.0*LOG(10.0) - 0.5*LOG(n_e) + LOG(t_e) - ! mu_i: ion mass ratio to proton - tau= t_i/t_e ! ratio of ion to electron temperature - tau_i = 6.6e17*mu_i**0.5*(t_i/1e3)**1.5/(n_e*lnLamb) ! ion colls. - eta= 1.65e-9*lnLamb/(t_e/1e3)**1.5 ! spitzer resistivity (wesson) - rho=(mu_i*m_p)*n_e ! mass density + tau = t_i / t_e ! T_i / T_e + tau_i = 6.6e17*mu_i**0.5*(t_i/1e3)**1.5 + $ / (n_e*lnLamb) ! ion collision time [s] + eta = 1.65e-9*lnLamb / (t_e/1e3)**1.5 ! Spitzer resistivity (Wesson) + rho = (mu_i*m_p)*n_e ! mass density [kg/m^3] - b_l=(nr/mr)*rs*sval*bt/R0 ! characteristic magnetic field - v_a=b_l/(mu0*rho)**0.5 ! alfven velocity - rho_s=1.02e-4*(mu_i*t_e)**0.5/bt ! ion Lamour by elec. Temp. +c----------------------------------------------------------------------- +c Electron-electron collision time (Braginskii) and Spitzer-Härm +c parallel conductivity. +c----------------------------------------------------------------------- + tau_ee_num = 6.0*SQRT(2.0)*(pi**1.5) + $ *(eps0**2.0)*(m_e**0.5)*(t_e**1.5) + tau_ee_denom = lnLamb*(chag**2.5)*n_e + tau_ee = tau_ee_num / tau_ee_denom - tau_h=R0*(mu0*rho)**0.5/(nn*sval*bt) ! alfven time across surface - tau_r=mu0*rs**2.0/eta ! resistive time scale - tau_v=tau_r/pr !rho*rs**2.0/visc ! viscous time scale - - ! this one must be anomalous. calculated back from pr. - visc= rho*rs**2.0/tau_v - - lu=tau_r/tau_h ! Lundquist number -! pr=tau_r/tau_v ! Prandtl number. Only place needs viscosity. - - omega_e=-t_e/(bt*R0)*(1.0/l_n+1.0/l_t)*qval ! elec. diamag - omega_i=t_i/(bt*R0)*(1.0/l_n+1.0/l_t)*qval ! ion diamag + sigma_par_1 = ( SQRT(2.0) + 13.0*(Zeff/4.0) ) + $ / (Zeff*(SQRT(2.0) + Zeff)) + sigma_par_2 = (n_e * (chag**2.0) * tau_ee) / m_e + sigma_par = sigma_par_1 * sigma_par_2 + +c----------------------------------------------------------------------- +c Characteristic field, Alfvén speed, length scales, and +c fundamental timescales. +c Note: mr, nr (rational m/n) and nn (toroidal mode number) are +c module-level variables set before calling params(). +c----------------------------------------------------------------------- + b_l = (nr/mr)*rs*sval*bt/R0 ! characteristic B [T] + v_a = b_l / (mu0*rho)**0.5 ! Alfvén velocity [m/s] + rho_s = 1.02e-4*(mu_i*t_e)**0.5 / bt ! ion Larmor at T_e [m] + d_i = ( (mu_i*m_p) / (n_e*(chag**2)*mu0) )**0.5 + ! ion skin depth [m] - ! now calculate the main 7 normalized parameters. + tau_h = R0*(mu0*rho)**0.5 / (nn*sval*bt) ! Alfvén time [s] + tau_r = mu0*(rs**2.0)*sigma_par ! resistive time [s] (Fitzpatrick) + IF (pr > 0.0_r8) THEN + tau_v = tau_r / pr ! viscous time [s] + visc = rho*rs**2.0_r8 / tau_v + ELSE + tau_v = 0.0_r8 + visc = 0.0_r8 + END IF - Qconv=lu**(1.0/3.0)*tau_h ! conversion to Qs based on Cole +c Lundquist number + lu = tau_r / tau_h - ! note Q depends on Qconv even if omega is fixed. - Q=Qconv*omega - Q_e=-Qconv*omega_e - Q_i=-Qconv*omega_i +c----------------------------------------------------------------------- +c Diamagnetic frequencies and normalised Q parameters. +c Qconv converts dimensional frequencies to the normalised Q +c used in the SLAYER dispersion relation (Cole scaling). +c----------------------------------------------------------------------- + Qconv = lu**(1.0/3.0) * tau_h ! frequency normalisation (Cole) + tauk = Qconv ! stored in sglobal_mod + Q = Qconv * omega ! normalised rotation frequency + +c Diamagnetic frequencies require gradient lengths; skip when the +c caller does not provide them (l_n <= 0 or l_t <= 0). Callers that +c compute omega_e/omega_i by other means must set Q_e and Q_i +c themselves afterwards. + IF (l_n > 0.0_r8 .AND. l_t > 0.0_r8) THEN + omega_e = -t_e/(bt*R0)*(1.0/l_n + 1.0/l_t)*qval + omega_i = t_i/(bt*R0)*(1.0/l_n + 1.0/l_t)*qval + Q_e = -Qconv * omega_e + Q_i = -Qconv * omega_i + END IF + +c normalised ion Larmor radius (critical stability parameter) + ds = lu**(1.0/3.0) * rho_s / rs + +c----------------------------------------------------------------------- +c Plasma beta and Prandtl-number-like transport ratios. +c----------------------------------------------------------------------- + lbeta = (5.0/3.0)*mu0*n_e*chag*(t_e+t_i) / bt**2.0 + c_beta = (lbeta / (1.0+lbeta))**0.5 + +c kappa-based P_perp (JKP's definition, to implement later if desired) +c K_val = chis(3) / eta +c Csq = c_beta**2.0 + (1.0 - c_beta**2.0)*K_val - ! This is the most critical parameter - ds=lu**(1.0/3.0)*rho_s/rs ! conversion based on Cole. +c IF (ABS(Csq) > 0.0) THEN +c P_perp = Csq +c ELSE +c tau_perp = (rs**2.0) / chis(1) +c END IF - lbeta=(5.0/3.0)*mu0*n_e*chag*(t_e+t_i)/bt**2.0 - c_beta=(lbeta/(1.0+lbeta))**0.5 +c effective perpendicular and toroidal Prandtl numbers + tau_perp = (rs**2.0) / chis(1) + P_perp = tau_r / tau_perp ! perp magnetic Prandtl number - delta_n=lu**(1.0/3.0)/rs ! norm factor for delta primes + tau_tor = (rs**2.0) / chis(2) + P_tor = tau_r / tau_tor ! toroidal magnetic Prandtl number - ! quick diagnostics. +c----------------------------------------------------------------------- +c Normalised beta-related width and Delta norm factor. +c----------------------------------------------------------------------- +c d_beta uses Fitzpatrick's definition (tau' form) + d_beta = c_beta * d_i + D_norm = (d_beta/rs) * lu**(1.0/3.0) + $ * (tau/(1+tau))**(0.5) + + delta_n = lu**(1.0/3.0) / rs ! normalisation for Delta values + +c----------------------------------------------------------------------- +c Critical Deltaprime (dc_tmp) via iterative chi_parallel calculation. +c The island-width Wd is iterated to converge the +c short-mfp / long-mfp interpolation for chi_parallel. +c dc_type (from sglobal_mod) selects the formula: +c 'lar' -- cylindrical (Lutjens) +c 'rfitzp' -- R. Fitzpatrick +c 'toroidal' -- toroidal geometry using dgeo_val +c default -- dc_tmp = 0 +c----------------------------------------------------------------------- + IF (ABS(dr_val) > 0.0) THEN + + vte = SQRT((2.0*(t_e*chag)) / m_e) + chi_par_smfp = (1.581*tau_ee*(vte**2.0)) + $ / (1.0 + 0.2535*Zeff) + + Wd = 0.1 ! initial guess + DO wit = 1, 100 + chi_par_lmfp = (2.0*R0*vte) + $ / (SQRT(pi)*nr*sval*Wd) + chi_par = (chi_par_smfp*chi_par_lmfp) + $ / (chi_par_smfp + chi_par_lmfp) + Wd_new = SQRT(8.0)*((chis(1)/chi_par)**0.25) + $ * (1.0/SQRT((rs/R0)*sval*nr)) + IF (ABS(Wd_new - Wd) / MAX(ABS(Wd), 1.0e-30_r8) + $ < 1.0e-10_r8) THEN + Wd = Wd_new + EXIT + END IF + Wd = Wd_new + END DO + IF (wit > 100) THEN + WRITE(*,*) 'params: Wd iteration failed to converge' + STOP 'params: Wd iteration did not converge' + END IF + + SELECT CASE(dc_type) + CASE('lar') + dc_tmp = 0.5*(-dr_val)*(pi**1.5) + $ *((chi_par/chis(1))**0.25) + $ *( (nr*sval)/(R0*rs) )**0.5 + CASE('rfitzp') + dc_tmp = -(SQRT(2.0)*(pi**(1.5))*dr_val) / Wd + CASE('toroidal') + dc_tmp = 0.5*(-dr_val)*(pi**1.5) + $ *((chi_par/chis(1))**0.25)*dgeo_val + CASE default + dc_tmp = 0.0 + END SELECT + + ELSE + dc_tmp = 0.0 + END IF + +c----------------------------------------------------------------------- +c optional diagnostics (guarded by params_check flag). +c----------------------------------------------------------------------- IF (params_check) THEN - WRITE(*,*)"eta=",eta - WRITE(*,*)"S=",lu - WRITE(*,*)"Q=",Q - WRITE(*,*)"Q_e=",Q_e - WRITE(*,*)"Q_i=",Q_i - WRITE(*,*)"ds=",ds - WRITE(*,*)"c_beta=",c_beta + WRITE(*,*) 'eta = ', eta + WRITE(*,*) 'S = ', lu + WRITE(*,*) 'Q = ', Q + WRITE(*,*) 'Q_e = ', Q_e + WRITE(*,*) 'Q_i = ', Q_i + WRITE(*,*) 'ds = ', ds + WRITE(*,*) 'c_beta = ', c_beta ENDIF - +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- RETURN - END SUBROUTINE params - - END MODULE params_mod + END SUBROUTINE params + + END MODULE params_mod \ No newline at end of file diff --git a/slayer/sglobal.f b/slayer/sglobal.f index 667cd235..bd4d5556 100644 --- a/slayer/sglobal.f +++ b/slayer/sglobal.f @@ -1,25 +1,193 @@ +c======================================================================= +c MODULE sglobal_mod +c +c Global shared state for the SLAYER layer-physics package. +c Contains: +c - physical and mathematical constants +c - module-level scalar variables written by params() and +c consumed by the dispersion solvers (delta, slayer, gslayer) +c - AMR scanner storage (hash-based v1, cell-based v2) +c - derived types: slayer_inputs_type, slayer_outputs_type, +c deltas_outputs_type, result_type, amr_cell_type +c +c TODO: `pr` and `pe` (magnetic Prandtl number and its +c electron analogue) are declared but never explicitly set +c within this module or in params(). They must be initialised +c before params() is called (tau_v = tau_r / pr). Consider +c adding a dedicated init routine or making them INTENT(IN) +c arguments of params(). +c +c======================================================================= MODULE sglobal_mod - USE local_mod, ONLY: r8 + + USE local_mod, ONLY: r8 ! double-precision kind parameter IMPLICIT NONE - INTEGER :: mm,nn - INTEGER :: in_unit,out_unit,out2_unit,out3_unit, - $ bin_unit,bin_2d_unit,input_unit - - REAL(r8) :: mr,nr - REAL(r8) :: Q_e,Q_i,pr,pe,c_beta,ds,tau - REAL(r8) :: eta,visc,rho_s,lu,omega_e,omega_i, - $ delta_n,layfac - COMPLEX(r8) :: Q - - REAL(r8), PARAMETER :: pi=3.1415926535897932385, mu0=4e-7*pi, - $ m_e=9.1094e-31,m_p=1.6726e-27,chag=1.6022e-19, - $ kval=1.3807e-23,lnLamb=17.0 - - ! lnLamb will be updated. - - COMPLEX(r8), PARAMETER :: ifac=(0,1) - - CHARACTER(2) :: sn - + +c----------------------------------------------------------------------- +c I/O unit numbers and trace control. +c----------------------------------------------------------------------- + INTEGER :: in_unit ! standard input unit + INTEGER :: out_unit ! primary output unit + INTEGER :: out2_unit ! secondary output unit + INTEGER :: out3_unit ! tertiary output unit + INTEGER :: bin_unit ! binary output unit + INTEGER :: bin_2d_unit ! 2-D binary output unit + INTEGER :: input_unit ! namelist / input unit + INTEGER :: n_trace ! trace / debug verbosity level + +c----------------------------------------------------------------------- +c Mode numbers (integer and real representations). +c----------------------------------------------------------------------- + INTEGER :: mm ! poloidal mode number (integer) + INTEGER :: nn ! toroidal mode number (integer) + REAL(r8) :: mr ! poloidal mode number (real copy) + REAL(r8) :: nr ! toroidal mode number (real copy) + CHARACTER(2) :: sn_str ! toroidal n as string + CHARACTER(2) :: sm_str ! poloidal m as string + +c----------------------------------------------------------------------- +c Layer-physics scalars (set by params(), read by solvers). +c These are the normalised parameters that enter the SLAYER +c dispersion relation. +c----------------------------------------------------------------------- +c --- temperature / collisionality + REAL(r8) :: tau ! T_i / T_e + REAL(r8) :: eta ! Spitzer resistivity [Ohm*m] + REAL(r8) :: visc ! anomalous viscosity [m^2/s] + REAL(r8) :: lnLamb ! Coulomb logarithm (updated at runtime) +c --- length scales + REAL(r8) :: rho_s ! ion Larmor radius at T_e [m] + REAL(r8) :: d_i ! ion skin depth [m] + REAL(r8) :: d_beta ! beta-weighted ion scale d_beta = c_beta * d_i +c --- timescales + REAL(r8) :: tau_r ! resistive diffusion time [s] + REAL(r8) :: tauk ! Q-conversion factor (= Qconv) +c --- Lundquist and Prandtl numbers + REAL(r8) :: lu ! Lundquist number S = tau_r / tau_h + REAL(r8) :: pr ! magnetic Prandtl number (TODO: see header) + REAL(r8) :: pe ! electron Prandtl number (TODO: see header) + REAL(r8) :: P_perp ! perpendicular magnetic Prandtl number + REAL(r8) :: P_tor ! toroidal magnetic Prandtl number +c --- normalised layer parameters + REAL(r8) :: ds ! normalised ion Larmor radius + REAL(r8) :: c_beta ! compressional beta parameter + REAL(r8) :: D_norm ! normalised beta-weighted ion scale + REAL(r8) :: delta_n ! Delta normalisation factor + REAL(r8) :: Qconv ! frequency normalisation (Cole) +c --- diamagnetic and rotation frequencies + REAL(r8) :: omega_e ! electron diamagnetic frequency [rad/s] + REAL(r8) :: omega_i ! ion diamagnetic frequency [rad/s] + COMPLEX(r8) :: Q ! normalised complex rotation frequency + REAL(r8) :: Q_e ! normalised electron diamagnetic Q + REAL(r8) :: Q_i ! normalised ion diamagnetic Q +c --- stability / Delta_crit + REAL(r8) :: dc_tmp ! computed Delta_crit + REAL(r8) :: delta_eff ! effective Deltaprime shift + CHARACTER(20) :: dc_type ! dc formula selector ('lar','rfitzp','toroidal') +c --- solver workspace / results + COMPLEX(r8) :: g_tmp ! temporary complex growth rate + REAL(r8) :: gamma_fac ! growth-rate conversion factor +c --- miscellaneous + REAL(r8) :: iota_e ! electron iota + REAL(r8) :: layfac ! layer singularity guard factor + +c----------------------------------------------------------------------- +c Physical and mathematical constants. +c----------------------------------------------------------------------- + REAL(r8), PARAMETER :: pi = 3.1415926535897932385d0 + REAL(r8), PARAMETER :: mu0 = 4.0d-7 * pi ! vacuum permeability [H/m] + REAL(r8), PARAMETER :: m_e = 9.1094d-31 ! electron mass [kg] + REAL(r8), PARAMETER :: m_p = 1.6726d-27 ! proton mass [kg] + REAL(r8), PARAMETER :: chag = 1.6021917d-19 ! elementary charge [C] + REAL(r8), PARAMETER :: kval = 1.3807d-23 ! Boltzmann constant [J/K] + REAL(r8), PARAMETER :: eps0 = 8.8542d-12 ! vacuum permittivity [F/m] + COMPLEX(r8), PARAMETER :: ifac = (0.0d0, 1.0d0) ! imaginary unit + +c----------------------------------------------------------------------- +c AMR scanner storage -- hash-based deduplication (v1). +c----------------------------------------------------------------------- + INTEGER, PARAMETER :: MAX_PTS = 1000000 ! max unique eval points + INTEGER, PARAMETER :: HASH_SZ = 1000003 ! hash table size (prime) + REAL(r8), PARAMETER :: HASH_SCALE = 1.0d5 ! Re/Im quantisation scale + INTEGER, ALLOCATABLE :: hash_head(:) ! bucket heads (HASH_SZ) + INTEGER, ALLOCATABLE :: hash_next(:) ! chain pointers (MAX_PTS) + +c----------------------------------------------------------------------- +c AMR scanner storage -- cell-based refinement (v2). +c----------------------------------------------------------------------- + INTEGER, PARAMETER :: MAX_CELLS = 1000000 ! max AMR cells + + TYPE :: amr_cell_type + COMPLEX(r8) :: Q(4) ! corner Q-values (BL, BR, TL, TR) + COMPLEX(r8) :: D(4) ! corner dispersion values + LOGICAL :: needs_refine ! flagged for subdivision + END TYPE amr_cell_type + + TYPE(amr_cell_type), ALLOCATABLE :: amr_cells(:) + INTEGER :: n_amr_cells ! current number of active cells + +c----------------------------------------------------------------------- +c AMR output arrays (shared by v1 and v2). +c----------------------------------------------------------------------- + COMPLEX(r8), ALLOCATABLE :: Q_store(:) ! unique Q-points + COMPLEX(r8), ALLOCATABLE :: D_store(:) ! corresponding D-values + INTEGER :: n_pts ! number of stored points + +c----------------------------------------------------------------------- +c Derived types: solver I/O and scan results. +c----------------------------------------------------------------------- + +c result_type -- torque-scan output bucket + TYPE result_type + REAL(r8), ALLOCATABLE :: inQs(:) ! Re(Q) scan values + REAL(r8), ALLOCATABLE :: iinQs(:) ! Im(Q) scan values + REAL(r8), ALLOCATABLE :: Re_deltas(:) ! Re(Delta) results + REAL(r8), ALLOCATABLE :: Im_deltas(:) ! Im(Delta) results + INTEGER :: count ! number of entries + END TYPE result_type + +c slayer_inputs_type -- per-surface input arrays for SLAYER + TYPE slayer_inputs_type + INTEGER, ALLOCATABLE :: qval_arr(:) ! safety-factor integers + REAL(r8), ALLOCATABLE :: chi_p_arr(:) ! chi_perp [m^2/s] + REAL(r8), ALLOCATABLE :: chi_t_arr(:) ! chi_tor [m^2/s] + REAL(r8), ALLOCATABLE :: kappa_arr(:) ! kappa (thermal cond.) + REAL(r8), ALLOCATABLE :: psi_n_arr(:) ! normalised psi + REAL(r8), ALLOCATABLE :: lu_arr(:) ! Lundquist number + REAL(r8), ALLOCATABLE :: Qconv_arr(:) ! Q-conversion factor + REAL(r8), ALLOCATABLE :: Q_e_arr(:) ! normalised Q_e + REAL(r8), ALLOCATABLE :: Q_i_arr(:) ! normalised Q_i + REAL(r8), ALLOCATABLE :: c_beta_arr(:) ! compressional beta + REAL(r8), ALLOCATABLE :: d_beta_arr(:) ! beta-related width + REAL(r8), ALLOCATABLE :: D_norm_arr(:) ! normalised D + REAL(r8), ALLOCATABLE :: tau_arr(:) ! T_i / T_e + REAL(r8), ALLOCATABLE :: P_perp_arr(:) ! perp Prandtl number + REAL(r8), ALLOCATABLE :: P_tor_arr(:) ! toroidal Prandtl + REAL(r8), ALLOCATABLE :: omegas_arr(:) ! rotation [rad/s] + REAL(r8), ALLOCATABLE :: omegas_e_arr(:) ! omega_e [rad/s] + REAL(r8), ALLOCATABLE :: omegas_i_arr(:) ! omega_i [rad/s] + REAL(r8), ALLOCATABLE :: gammafac_arr(:) ! gamma conversion + REAL(r8), ALLOCATABLE :: Re_dp_arr(:) ! Re(Deltaprime) + REAL(r8), ALLOCATABLE :: Im_dp_arr(:) ! Im(Deltaprime) + REAL(r8), ALLOCATABLE :: d_crit_arr(:) ! Delta_crit + COMPLEX(r8), ALLOCATABLE :: dp_matrix(:,:) ! full Deltaprime matrix + END TYPE slayer_inputs_type + +c slayer_outputs_type -- per-surface solver results + TYPE slayer_outputs_type + COMPLEX(r8), ALLOCATABLE :: dels_db_arr(:) ! Delta from d_beta + COMPLEX(r8), ALLOCATABLE :: gamma_sol_arr(:) ! solved growth rate + COMPLEX(r8), ALLOCATABLE :: gamma_est_arr(:) ! estimated growth rate + REAL(r8), ALLOCATABLE :: br_th_arr(:) ! Br threshold + END TYPE slayer_outputs_type + +c deltas_outputs_type -- scan output (Q vs Delta) + TYPE deltas_outputs_type + REAL(r8), ALLOCATABLE :: inQs(:) ! Re(Q) values + REAL(r8), ALLOCATABLE :: iinQs(:) ! Im(Q) values + REAL(r8), ALLOCATABLE :: real_deltas(:) ! Re(Delta) + REAL(r8), ALLOCATABLE :: imag_deltas(:) ! Im(Delta) + END TYPE deltas_outputs_type + END MODULE sglobal_mod diff --git a/slayer/slayer.f b/slayer/slayer.f index c7a4cd8f..fbc66e43 100644 --- a/slayer/slayer.f +++ b/slayer/slayer.f @@ -1,11 +1,10 @@ c----------------------------------------------------------------------- -c Slab LAYER based on linear drift MHD -c SLAYER: main program -c----------------------------------------------------------------------- -c----------------------------------------------------------------------- -c code organization. -c----------------------------------------------------------------------- -c slayer. +c SLAYER: Slab LAYER linear drift-MHD code. +c Main driver program. +c +c Computes tearing-mode layer quantities (inner layer Delta, +c growth rates, torque balance, field thresholds) from +c slab-geometry drift-MHD matching via a Riccati integration method. c----------------------------------------------------------------------- c----------------------------------------------------------------------- c declarations. @@ -13,125 +12,321 @@ PROGRAM slayer USE sglobal_mod - USE params_mod - USE delta_mod, ONLY: riccati,riccati_out, - $ parflow_flag,PeOhmOnly_flag + USE delta_mod, ONLY: riccati,riccati_f,riccati_del_s, + $ riccati_out,parflow_flag,PeOhmOnly_flag + USE gslayer_mod + USE growthrates_mod + USE layerinputs_mod IMPLICIT NONE - - CHARACTER(128) :: infile - INTEGER :: i,j,k,inum,jnum,knum,inn - INTEGER, DIMENSION(1) :: index - - LOGICAL :: params_flag,QPscan_flag,QPescan_flag,QPscan2_flag, - $ QDscan2_flag,Qbscan_flag,Qscan_flag, - $ onscan_flag,otscan_flag,ntscan_flag,nbtscan_flag, - $ Pe_flag,verbose,ascii_flag,bin_flag,netcdf_flag, - $ bal_flag,stability_flag,riccatiscan_flag,input_flag, - $ params_check - - REAL(r8) :: n_e,t_e,t_i,omega,omega0, - $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff - REAL(r8) :: inQ,inQ_e,inQ_i,inpr,inpe,inc_beta,inds,intau,inlu - REAL(r8) :: psi0,jxb,Q0,Q_sol,br_th - COMPLEX(r8) :: delta,delta_n_p - - REAL(r8) :: inQ_min,inQ_max,j_min,j_max,jpower,k_min,k_max,kpower, - $ Qratio - - INTEGER, DIMENSION(:), ALLOCATABLE :: mms,nns - - REAL(r8), DIMENSION(8) :: inpr_prof - - REAL(r8), DIMENSION(:), ALLOCATABLE :: inQs,iinQs,jxbl,bal, - $ prs,n_es,t_es,t_is,omegas,l_ns,l_ts,qvals,svals, - $ bts,rss,R0s,mu_is,zeffs,Q_soll,br_thl - REAL(r8), DIMENSION(:,:), ALLOCATABLE :: - $ js,ks,psis,jxbs,Q_sols,br_ths +c----------------------------------------------------------------------- +c local scalars — loop indices and counters. +c----------------------------------------------------------------------- + INTEGER :: i,j,k ! general loop indices + INTEGER :: inn ! number of input-file surfaces + INTEGER :: count ! generic counter +c----------------------------------------------------------------------- +c local scalars — numerical resolution. +c----------------------------------------------------------------------- + INTEGER :: inum ! resolution for 1-D scans + INTEGER :: jnum ! resolution for 2-D scan axis 1 + INTEGER :: knum ! resolution for 2-D scan axis 2 + INTEGER :: Q_num ! resolution for stab. scan Re(Q) + INTEGER :: msing_max ! max number of singular surfaces + INTEGER :: n_k ! number of rational surfaces +c----------------------------------------------------------------------- +c local scalars — MAXLOC result holder. +c----------------------------------------------------------------------- + INTEGER, DIMENSION(1) :: iloc ! result of MAXLOC +c----------------------------------------------------------------------- +c control flags — workflow. +c----------------------------------------------------------------------- + LOGICAL :: params_flag ! compute params from kinetic data + LOGICAL :: input_flag ! read multi-surface input file + LOGICAL :: read_eq ! read equilibrium files + LOGICAL :: verbose ! enable progress messages + LOGICAL :: params_check ! print diagnostic output in params +c----------------------------------------------------------------------- +c control flags — physics modes. +c----------------------------------------------------------------------- + LOGICAL :: est_gamma_flag ! estimate growth rate + LOGICAL :: match_gamma_flag ! asymptotically matched gamma + LOGICAL :: Pperp_Ptor_flag ! use Fitzpatrick layer model + LOGICAL :: coupling_flag ! coupled rational surfaces + LOGICAL :: br_th_flag ! Br threshold test scan + LOGICAL :: bal_flag ! torque balance scan + LOGICAL :: stability_flag ! complex-Q delta scan + LOGICAL :: Pe_flag ! include electron pressure +c----------------------------------------------------------------------- +c control flags — parameter-space scans. +c----------------------------------------------------------------------- + LOGICAL :: QPscan_flag ! (Q,P) scan + LOGICAL :: QPescan_flag ! (Q,Pe) scan + LOGICAL :: QPscan2_flag ! (Q,P) scan variant 2 + LOGICAL :: QDscan2_flag ! (Q,D) scan variant 2 + LOGICAL :: Qbscan_flag ! (Q,beta) scan + LOGICAL :: Qscan_flag ! 1-D Q scan + LOGICAL :: onscan_flag ! (omega,n) scan + LOGICAL :: otscan_flag ! (omega,T) scan + LOGICAL :: ntscan_flag ! (n,T) scan + LOGICAL :: nbtscan_flag ! (n,Bt) scan + LOGICAL :: riccatiscan_flag ! Riccati-variable scan + LOGICAL :: stabscan_flag ! stability scan (single surface) + LOGICAL :: coupled_stabscan_flag ! stability scan (coupled) + LOGICAL :: amr_flag ! adaptive mesh refinement scan +c----------------------------------------------------------------------- +c control flags — output format. +c----------------------------------------------------------------------- + LOGICAL :: ascii_flag ! write ASCII output files + LOGICAL :: bin_flag ! write binary output files + LOGICAL :: netcdf_flag ! write NetCDF output files +c----------------------------------------------------------------------- +c local scalars — physical input quantities. +c----------------------------------------------------------------------- + REAL(r8) :: n_e ! electron density [m^-3] + REAL(r8) :: t_e ! electron temperature [eV] + REAL(r8) :: t_i ! ion temperature [eV] + REAL(r8) :: omega ! toroidal rotation [rad/s] + REAL(r8) :: omega0 ! unused + REAL(r8) :: l_n ! density gradient scale length + REAL(r8) :: l_t ! temperature gradient scale length + REAL(r8) :: qval ! safety factor at surface + REAL(r8) :: sval ! magnetic shear at surface + REAL(r8) :: bt ! toroidal field [T] + REAL(r8) :: rs ! minor radius of surface [m] + REAL(r8) :: R0 ! major radius [m] + REAL(r8) :: mu_i ! ion mass number + REAL(r8) :: zeff ! effective charge + REAL(r8) :: dr_val ! radial derivative parameter + REAL(r8) :: dgeo_val ! geometric factor parameter + REAL(r8) :: scan_width ! half-width of complex-Q scan +c----------------------------------------------------------------------- +c local scalars — normalized layer parameters (namelist overrides). +c----------------------------------------------------------------------- + REAL(r8) :: inQ ! normalized ExB rotation freq. + REAL(r8) :: inQ_e ! normalized electron diamagnetic + REAL(r8) :: inQ_i ! normalized ion diamagnetic + REAL(r8) :: inQ_e_ovr(8) ! per-surface Q_e overrides (read_eq) + REAL(r8) :: inQ_i_ovr(8) ! per-surface Q_i overrides (read_eq) + REAL(r8) :: inpr ! normalized pressure gradient + REAL(r8) :: inpe ! normalized electron pressure + REAL(r8) :: inc_beta ! normalized beta + REAL(r8) :: inds ! normalized D (magnetic diffusion) + REAL(r8) :: intau ! normalized tau = T_i/T_e + REAL(r8) :: inlu ! normalized Lundquist number +c----------------------------------------------------------------------- +c local scalars — derived / scratch quantities. +c----------------------------------------------------------------------- + REAL(r8) :: psi0 ! reconnected flux (a.u.) + REAL(r8) :: jxb ! j x B torque (a.u.) + REAL(r8) :: Q0 ! unperturbed rotation frequency + REAL(r8) :: Q_sol ! solved rotation frequency + REAL(r8) :: br_th ! radial field threshold + REAL(r8) :: Qratio ! Q_e/Q ratio for scan2 variants +c----------------------------------------------------------------------- +c local scalars — scan grid helpers. +c----------------------------------------------------------------------- + REAL(r8) :: inQ_min,inQ_max ! rotation scan bounds + REAL(r8) :: j_min,j_max,jpower ! 2-D scan axis 1 + REAL(r8) :: k_min,k_max,kpower ! 2-D scan axis 2 + REAL(r8) :: ing_step ! growth-rate grid step + REAL(r8) :: ing_coarse ! Re(gamma) grid value + REAL(r8) :: iing_coarse ! Im(gamma) grid value +c----------------------------------------------------------------------- +c local scalars — complex quantities. +c----------------------------------------------------------------------- + COMPLEX(r8) :: delta ! layer Delta (tearing index) + COMPLEX(r8) :: delta_n_p ! Deltaprime scale factor + COMPLEX(r8) :: dels_db ! delta_s / d_beta + COMPLEX(r8) :: del_s ! delta_s + COMPLEX(r8) :: ingamma ! initial gamma guess (namelist) + COMPLEX(r8) :: delta_prime ! external Deltaprime (namelist) +c----------------------------------------------------------------------- +c local arrays — transport profile coefficients. +c----------------------------------------------------------------------- + REAL(r8) :: chis(3) ! chi_perp, chi_tor, kappa + REAL(r8), DIMENSION(8) :: chi_p_prof ! chi_perp radial profile + REAL(r8), DIMENSION(8) :: chi_t_prof ! chi_tor radial profile + REAL(r8), DIMENSION(8) :: kappa_prof ! kappa radial profile +c----------------------------------------------------------------------- +c local arrays — multi-surface input-file storage. +c----------------------------------------------------------------------- + INTEGER, DIMENSION(:), ALLOCATABLE :: mms,nns + REAL(r8), DIMENSION(:), ALLOCATABLE :: prs,n_es,t_es,t_is, + $ omegas,l_ns,l_ts,svals,qvals,bts,rss,R0s,mu_is,zeffs, + $ Q_soll,br_thl,pes +c----------------------------------------------------------------------- +c local arrays — scan workspace. +c----------------------------------------------------------------------- + REAL(r8), DIMENSION(:), ALLOCATABLE :: inQs,iinQs + REAL(r8), DIMENSION(:), ALLOCATABLE :: jxbl,bal + REAL(r8), DIMENSION(:,:), ALLOCATABLE :: js,ks,psis,jxbs, + $ Q_sols,br_ths REAL(r8), DIMENSION(:,:,:), ALLOCATABLE :: Q_solss,br_thss - COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: deltal + COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: deltal,outer_deltas COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: deltas +c----------------------------------------------------------------------- +c AMR (adaptive mesh refinement) growth-rate scan. +c----------------------------------------------------------------------- + INTEGER :: AMR_passes ! number of AMR refinement levels + INTEGER :: m_AMR ! effective number of surfaces +c----------------------------------------------------------------------- +c structured input/output types (defined in sglobal_mod). +c----------------------------------------------------------------------- + TYPE(slayer_inputs_type) :: sl_in + TYPE(slayer_outputs_type) :: sl_out + TYPE(deltas_outputs_type), ALLOCATABLE :: all_deltas_out(:) +c----------------------------------------------------------------------- +c file-path strings. +c----------------------------------------------------------------------- + CHARACTER(512) :: infile ! multi-surface input file path + CHARACTER(512) :: ncfile ! NetCDF equilibrium file path - NAMELIST/slayer_input/params_flag,input_flag,infile, - $ mm,nn,n_e,t_e,t_i,omega,l_n,l_t, - $ qval,sval,bt,rs,R0,zeff,mu_i,inQ,inQ_e,inQ_i, - $ inpr,inpr_prof,inpe,inc_beta,inds,intau,inlu,Q0,delta_n_p - NAMELIST/slayer_control/inum,jnum,knum,QPscan_flag,QPscan2_flag, - $ QPescan_flag,QDscan2_flag,Qbscan_flag,Qscan_flag, - $ onscan_flag,otscan_flag,ntscan_flag,nbtscan_flag, - $ layfac,Qratio,parflow_flag,peohmonly_flag,Pe_flag +c----------------------------------------------------------------------- +c namelist groups. +c----------------------------------------------------------------------- + NAMELIST/slayer_input/input_flag,infile, + $ ncfile,params_flag,mm,nn,n_e,t_e,t_i,sval,bt,rs,R0,omega, + $ l_t,l_n,qval,mu_i,zeff,dr_val,dgeo_val,chi_p_prof, + $ chi_t_prof,kappa_prof,inpr,inpe,inQ,inQ_e,inQ_i,inc_beta, + $ inds,intau,Q0,delta_prime,delta_n_p,ingamma, + $ inQ_e_ovr,inQ_i_ovr + NAMELIST/slayer_control/inum,jnum,knum,Q_num,scan_width, + $ AMR_passes,msing_max,dc_type,read_eq,Pperp_Ptor_flag, + $ coupling_flag,QPscan_flag,Qscan_flag,QPescan_flag, + $ Qbscan_flag,onscan_flag,otscan_flag,ntscan_flag, + $ nbtscan_flag,parflow_flag,peohmonly_flag,Pe_flag,layfac NAMELIST/slayer_output/verbose,ascii_flag,bin_flag,netcdf_flag, - $ stability_flag + $ est_gamma_flag,match_gamma_flag,stability_flag, + $ stabscan_flag,coupled_stabscan_flag,amr_flag,br_th_flag, + $ bal_flag NAMELIST/slayer_diagnose/riccati_out,riccatiscan_flag, - $ params_check,bal_flag + $ params_check c----------------------------------------------------------------------- c set initial values. +c defaults are overridden by namelist reads below. c----------------------------------------------------------------------- - mm=2 - nn=1 - mr = real(mm,4) - nr = real(nn,4) - n_e=1e19 - t_e=1e3 - t_i=1e3 - omega=1e4 - l_n=2e-1 - l_t=2e-1 - qval=2.0 - sval=2.0 - bt=1.0 - rs=0.5 - R0=1.0 - mu_i=2.0 - zeff=2.0 - inQ=20.0 ! Q=23.0 for DIII-D example. - inQ_e=2.0 ! Q_e=2.0 for DIII-D example. - inQ_i=-2.6 ! Q_i=-2.6 for DIII-D example. - inpr=1.0 ! - inpr_prof=1.0 - inpe=0.0 - inc_beta=0.7 ! c_beta=0.7 for DIII-D example. - inds=6.0 ! 6.0 for DIII-D example. - intau=1.0 ! 1.0 for DIII-D example. - inlu=1e8 - Q0=4.0 - delta_n_p=(1e-2,1e-2) - inum=400 ! resolution to find error field thresholds. - jnum=500 ! resolution for 2d scan along with Q,omega. - knum=100 ! resolution for 2d scan alont with the other. - in_unit=1 - out_unit=2 - out2_unit=3 - out3_unit=4 - bin_unit=5 - bin_2d_unit=6 - input_unit=7 - QPscan_flag=.FALSE. ! scan (Q,P) space for delta and torque. - QPescan_flag=.FALSE. ! scan (Q,Pe) space for delta and torque. - Qbscan_flag=.FALSE. ! scan (Q,beta) space for delta and torque. - onscan_flag=.FALSE. ! scan (omega,n) space for error fields. - otscan_flag=.FALSE. ! scan (omega,t) space for error fields. - ntscan_flag=.FALSE. ! scan (n,te) space for error fields. - nbtscan_flag=.FALSE. ! scan (n,bt) space for error fields. - layfac=0.02 - Qratio=0.5 - parflow_flag=.FALSE. - PeOhmOnly_flag=.TRUE. - Pe_flag=.FALSE. - params_flag=.TRUE. - input_flag=.FALSE. - infile="input_params.dat" - verbose=.TRUE. - ascii_flag=.TRUE. - bin_flag=.TRUE. - netcdf_flag=.FALSE. - riccati_out=.FALSE. - riccatiscan_flag=.FALSE. - params_check=.FALSE. - bal_flag=.FALSE. - stability_flag=.FALSE. + + ! mode numbers (sglobal_mod: INTEGER mm,nn; REAL mr,nr) + mm = 0 + nn = 0 + mr = 0.0 + nr = 0.0 + + ! kinetic / equilibrium inputs + n_e = 0.0 + t_e = 0.0 + t_i = 0.0 + omega = 0.0 + l_n = 0.0 + l_t = 0.0 + qval = 0.0 + sval = 0.0 + bt = 0.0 + rs = 0.0 + R0 = 0.0 + mu_i = 0.0 + zeff = 0.0 + dr_val = 0.0 + dgeo_val= 0.0 + + ! normalized layer-parameter overrides + inQ = 0.0 + inQ_e = 0.0 + inQ_i = 0.0 + inQ_e_ovr = 0.0 + inQ_i_ovr = 0.0 + inpr = 0.0 + inpe = 0.0 + inc_beta = 0.0 + inds = 0.0 + intau = 0.0 + inlu = 0.0 + Q0 = 0.0 + + ! transport profile coefficients + chi_p_prof = 0.0 + chi_t_prof = 0.0 + kappa_prof = 0.0 + chis = 0.0 + + ! complex namelist inputs + delta_prime = (0.0,0.0) + delta_n_p = (0.0,0.0) + ingamma = (0.0,0.0) + + ! global module scalars (sglobal_mod) + gamma_fac = 0.0 + dc_type = "" + + ! scan resolution defaults + inum = 400 ! 1-D resolution (error-field threshold scans) + jnum = 500 ! 2-D scan axis-1 resolution + knum = 100 ! 2-D scan axis-2 resolution + Q_num = 100 ! stability scan Re(Q) resolution + scan_width = 2.0 ! half-width for complex-Q scans + AMR_passes = 4 ! AMR refinement levels + msing_max = 2 ! max singular surfaces to process + + ! I/O unit numbers (sglobal_mod) + in_unit = 1 + out_unit = 2 + out2_unit = 3 + out3_unit = 4 + bin_unit = 5 + bin_2d_unit= 6 + input_unit = 7 + + ! workflow flags + read_eq = .FALSE. + est_gamma_flag = .FALSE. + match_gamma_flag = .FALSE. + Pperp_Ptor_flag = .FALSE. + coupling_flag = .FALSE. + params_flag = .TRUE. + input_flag = .FALSE. + + ! parameter-space scan flags + QPscan_flag = .FALSE. + QPescan_flag = .FALSE. + Qbscan_flag = .FALSE. + onscan_flag = .FALSE. + otscan_flag = .FALSE. + ntscan_flag = .FALSE. + nbtscan_flag = .FALSE. + + ! physics / model flags + layfac = 0.02 + Qratio = 0.5 + parflow_flag = .FALSE. + PeOhmOnly_flag = .TRUE. + Pe_flag = .FALSE. + + ! file paths + infile = "" + ncfile = "" + + ! output control + verbose = .TRUE. + ascii_flag = .TRUE. + bin_flag = .TRUE. + netcdf_flag = .FALSE. + + ! diagnostic flags + riccati_out = .FALSE. + riccatiscan_flag = .FALSE. + params_check = .FALSE. + + ! remaining physics-mode flags + bal_flag = .FALSE. + stability_flag = .FALSE. + stabscan_flag = .FALSE. + coupled_stabscan_flag= .FALSE. + amr_flag = .FALSE. + br_th_flag = .FALSE. c----------------------------------------------------------------------- c read slayer.in. +c four namelist groups: input, control, output, diagnose. c----------------------------------------------------------------------- IF(verbose) WRITE(*,*)"" IF(verbose) WRITE(*,*)"SLAYER START" @@ -143,18 +338,28 @@ PROGRAM slayer READ(in_unit,NML=slayer_diagnose) CLOSE(UNIT=in_unit) + ! Build toroidal-mode-number string for output filenames. IF (nn<10) THEN - WRITE(UNIT=sn,FMT='(I1)') nn - sn=ADJUSTL(sn) + WRITE(UNIT=sn_str,FMT='(I1)') nn + sn_str=ADJUSTL(sn_str) ELSE - WRITE(UNIT=sn,FMT='(I2)') nn + WRITE(UNIT=sn_str,FMT='(I2)') nn ENDIF c----------------------------------------------------------------------- -c calculate parameters as needed. +c compute normalized layer parameters from kinetic inputs. +c params() (params_mod) converts dimensional plasma profiles into +c the normalized quantities (Q, Q_e, Q_i, c_beta, ds, tau, lu) +c used by the Riccati solver. c----------------------------------------------------------------------- IF (params_flag) THEN - CALL params(n_e,t_e,t_i,omega, + nr = nn ! ensure toroidal mode number is set for Wd iteration + mr = mm ! ensure poloidal mode number is set + chis(1) = chi_p_prof(1) ! chi_perp + chis(2) = chi_t_prof(1) ! chi_tor + chis(3) = kappa_prof(1) ! kappa + CALL params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + ! Copy module-level results into local working variables. inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -163,28 +368,39 @@ PROGRAM slayer intau=tau Q0=Q ELSE - lu=inlu + lu=inlu ! manual Lundquist number when params not computed ENDIF c----------------------------------------------------------------------- -c calculate basic delta, torque, balance, error fields. +c baseline single-surface delta, reconnected flux, & torque. +c skipped when the matched-gamma path is active (it computes +c its own delta internally). c----------------------------------------------------------------------- - delta=riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) - psi0=1.0/ABS(delta+delta_n_p) ! a.u. - jxb=-AIMAG(1.0/(delta+delta_n_p)) ! a.u. - WRITE(*,*)"delta=",delta - WRITE(*,*)"psi0=",psi0 - WRITE(*,*)"jxb=",jxb + IF (.NOT. (match_gamma_flag)) THEN + delta=riccati(inQ,inQ_e,inQ_i,inpr,inc_beta,inds,intau,inpe) + psi0=1.0/ABS(delta+delta_n_p) ! reconnected flux [a.u.] + jxb=-AIMAG(1.0/(delta+delta_n_p)) ! j x B torque [a.u.] + IF (verbose) THEN + WRITE(*,*)"delta=",delta + WRITE(*,*)"psi0=",psi0 + WRITE(*,*)"jxb=",jxb + ENDIF + ENDIF c----------------------------------------------------------------------- -c calculate parameters as needed. +c multi-surface input-file mode. +c reads an external file of (m,n) surfaces with per-surface kinetic +c profiles, computes delta and the error-field threshold for each. c----------------------------------------------------------------------- IF (input_flag) THEN OPEN(UNIT=input_unit,FILE=infile,STATUS="old") READ(input_unit,*)inn - ALLOCATE(mms(inn),nns(inn),prs(inn), - $ n_es(inn),t_es(inn),t_is(inn),omegas(inn), - $ l_ns(inn),l_ts(inn),qvals(inn),svals(inn), - $ bts(inn),rss(inn),R0s(inn),mu_is(inn),zeffs(inn), - $ Q_soll(inn),br_thl(inn)) + ALLOCATE(mms(0:inn-1),nns(0:inn-1),prs(0:inn-1), + $ n_es(0:inn-1),t_es(0:inn-1),t_is(0:inn-1), + $ omegas(0:inn-1), + $ l_ns(0:inn-1),l_ts(0:inn-1),qvals(0:inn-1), + $ svals(0:inn-1), + $ bts(0:inn-1),rss(0:inn-1),R0s(0:inn-1), + $ mu_is(0:inn-1),zeffs(0:inn-1), + $ Q_soll(0:inn-1),br_thl(0:inn-1)) DO k=0,inn-1 READ(input_unit,'(2(1x,I2),14(1x,e12.4))') $ mms(k),nns(k),prs(k), @@ -195,13 +411,13 @@ PROGRAM slayer CLOSE(input_unit) DO k=0,inn-1 - WRITE(*,*)k + WRITE(*,*)k ! surface index mr=REAL(mms(k)) nr=REAL(nns(k)) inpr=prs(k) - CALL params(n_es(k),t_es(k),t_is(k),omegas(k), - $ l_ns(k),l_ts(k),qvals(k),svals(k),bts(k),rss(k),R0s(k), - $ mu_is(k),zeffs(k),params_check) + CALL params(n_es(k),t_es(k),t_is(k),omegas(k),chis,dr_val, + $ dgeo_val,l_ns(k),l_ts(k),qvals(k),svals(k),bts(k), + $ rss(k),R0s(k),mu_is(k),zeffs(k),params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -234,16 +450,16 @@ PROGRAM slayer bal(i)=2.0*inpr*(Q0-inQs(i))/jxbl(i) ENDDO - index=MAXLOC(bal) - Q_soll(k)=inQs(index(1)) + iloc=MAXLOC(bal) + Q_soll(k)=inQs(iloc(1)) br_thl(k)=sqrt(MAXVAL(bal)/lu*(svals(k)**2.0/2.0))*1e4 - WRITE(*,*)"Q_sol=",Q_soll(k) - WRITE(*,*)"br_th=",br_thl(k) + IF (verbose) WRITE(*,*)"Q_sol=",Q_soll(k) + IF (verbose) WRITE(*,*)"br_th=",br_thl(k) DEALLOCATE(inQs,deltal,jxbl,bal) ENDDO OPEN(UNIT=out_unit,FILE="slayer_input_bal_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,(2a17))') "Q_sol","br_th" DO k=0,inn-1 @@ -255,6 +471,613 @@ PROGRAM slayer $ bts,rss,R0s,mu_is,zeffs,Q_soll,br_thl,mms,nns) ENDIF c----------------------------------------------------------------------- +c estimate growth rate via resistive-layer thickness. +c Uses riccati_del_s to get delta_s/d_beta, then scales by +c d_beta to obtain the layer thickness delta_s and the estimated +c gamma. Inputs may come from equilibrium files (read_eq) or +c namelist. Subroutines: build_inputs (layerinputs_mod), +c allocate_inputs, allocate_outputs (growthrates_mod). +c----------------------------------------------------------------------- + IF (est_gamma_flag) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,*)">>> Estimating growth rate" + + IF (read_eq) THEN + ! Read equilibrium files for multi-surface inputs. + ! build_inputs (layerinputs_mod) reads STRIDE NetCDF data. + sl_in%chi_p_arr = chi_p_prof + sl_in%chi_t_arr = chi_t_prof + sl_in%kappa_arr = kappa_prof + + CALL build_inputs(infile,ncfile,sl_in) + + n_k = SIZE(sl_in%qval_arr) + CALL allocate_outputs(n_k,sl_out) + +c Override Q_e/Q_i with namelist values if nonzero. +c inQ_e overrides all surfaces; per-surface arrays take +c precedence (surface k gets inQ_e_ovr(k) if nonzero). + IF (ABS(inQ_e) > 0.0) THEN + WRITE(*,*) 'Overriding all Q_e with inQ_e=',inQ_e + sl_in%Q_e_arr(:) = inQ_e + END IF + IF (ABS(inQ_i) > 0.0) THEN + WRITE(*,*) 'Overriding all Q_i with inQ_i=',inQ_i + sl_in%Q_i_arr(:) = inQ_i + END IF + DO k = 1, MIN(n_k, SIZE(inQ_e_ovr)) + IF (ABS(inQ_e_ovr(k)) > 0.0) THEN + WRITE(*,*) ' Q_e override surface',k,':', + $ sl_in%Q_e_arr(k),' ->',inQ_e_ovr(k) + sl_in%Q_e_arr(k) = inQ_e_ovr(k) + END IF + END DO + DO k = 1, MIN(n_k, SIZE(inQ_i_ovr)) + IF (ABS(inQ_i_ovr(k)) > 0.0) THEN + sl_in%Q_i_arr(k) = inQ_i_ovr(k) + END IF + END DO + + ELSE + ! Single-surface mode: build inputs from namelist. + n_k = 1 + mr = mm + nr = nn + + chis(1) = chi_p_prof(1) ! chi_perp + chis(2) = chi_t_prof(1) ! chi_tor + chis(3) = kappa_prof(1) ! kappa (thermal cond.) + + CALL params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, + $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + + ! Override computed parameters with nonzero namelist values. + IF (ABS(inQ) > 0.0) THEN + Q = inQ ! NAMELIST + END IF + IF (ABS(inQ_e) > 0.0) THEN + Q_e = inQ_e ! NAMELIST + END IF + IF (ABS(inQ_i) > 0.0) THEN + Q_i = inQ_i ! NAMELIST + END IF + IF (inpr > 0.0) THEN + pr = inpr ! NAMELIST + END IF + IF (intau > 0.0) THEN + tau = intau ! NAMELIST + END IF + IF (inds > 0.0) THEN + D_norm = inds ! NAMELIST + END IF + + CALL allocate_inputs(n_k,sl_in) ! growthrates_mod + CALL allocate_outputs(n_k,sl_out) ! growthrates_mod + + sl_in%qval_arr = (/ qval /) + sl_in%omegas_arr = (/ omega /) + !sl_in%Q_arr = (/ Q /) + sl_in%Q_e_arr = (/ Q_e /) + sl_in%Q_i_arr = (/ Q_i /) + sl_in%psi_n_arr = (/ 0.0 /) + sl_in%Re_dp_arr = (/ REAL(delta_prime) /) + sl_in%Im_dp_arr = (/ AIMAG(delta_prime) /) + sl_in%d_crit_arr = (/ dc_tmp /) + sl_in%P_perp_arr = (/ P_perp /) + sl_in%P_tor_arr = (/ P_tor /) + sl_in%tau_arr = (/ tau /) + sl_in%D_norm_arr = (/ D_norm /) + sl_in%d_beta_arr = (/ d_beta /) + sl_in%gammafac_arr = (/ gamma_fac /) + sl_in%c_beta_arr = (/ c_beta /) + sl_in%lu_arr = (/ lu /) + sl_in%Qconv_arr = (/ tauk /) + END IF + + ! Loop over rational surfaces to estimate growth rates. + DO k=1,n_k + WRITE(*,*) + WRITE(*,'(A,I0,A)') 'Calculating growth rate '// + $ 'estimate on q = ', + $ sl_in%qval_arr(k),' rational surface' + + D_norm = sl_in%D_norm_arr(k) + ! First arg is Q_e (electron diamagnetic freq), not Q + ! (ExB freq). This is intentional per riccati_del_s API. + dels_db=riccati_del_s(sl_in%Q_e_arr(k), + $ sl_in%Q_i_arr(k),sl_in%P_perp_arr(k), + $ 5.0*sl_in%D_norm_arr(k)) + + del_s = dels_db * sl_in%d_beta_arr(k) + + sl_out%gamma_est_arr(k) = sl_in%gammafac_arr(k)/del_s + sl_out%dels_db_arr(k) = dels_db + WRITE(*,*) + WRITE(*,'(A,F0.3,A)')'Growth rate estimate = ', + $ REAL(sl_out%gamma_est_arr(k)),' [Hz]' + + ENDDO + + IF (.NOT. (match_gamma_flag)) THEN + sl_out%gamma_sol_arr = (/0./) + CALL output_gamma(est_gamma_flag,m_AMR,sl_in,sl_out, + $ all_deltas_out) + END IF + ENDIF +c----------------------------------------------------------------------- +c asymptotically matched growth rate. +c Matches the inner-layer Delta to the outer-region Delta' to +c find the self-consistent complex growth rate. Supports both +c single-surface and coupled multi-surface (AMR) modes. +c Subroutines: dispersion_AMR_v2, dispersion_det (growthrates_mod), +c riccati_f (delta_mod). +c----------------------------------------------------------------------- + IF (match_gamma_flag) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,*)">>> Calculating asymptotically matched growth rate" + + IF (read_eq) THEN + + IF (.NOT. est_gamma_flag) THEN + sl_in%chi_p_arr = chi_p_prof + sl_in%chi_t_arr = chi_t_prof + sl_in%kappa_arr = kappa_prof + + CALL build_inputs(infile,ncfile,sl_in) + + n_k = SIZE(sl_in%qval_arr) + CALL allocate_outputs(n_k,sl_out) + END IF +c Apply Q_e/Q_i overrides (same logic as est_gamma path) + DO k = 1, MIN(n_k, SIZE(inQ_e_ovr)) + IF (ABS(inQ_e_ovr(k)) > 0.0) THEN + sl_in%Q_e_arr(k) = inQ_e_ovr(k) + END IF + END DO + DO k = 1, MIN(n_k, SIZE(inQ_i_ovr)) + IF (ABS(inQ_i_ovr(k)) > 0.0) THEN + sl_in%Q_i_arr(k) = inQ_i_ovr(k) + END IF + END DO + ELSE + n_k = 1 + + IF (.NOT. est_gamma_flag) THEN + + chis(1) = chi_p_prof(1) ! chi_perp + chis(2) = chi_t_prof(1) ! chi_tor + chis(3) = kappa_prof(1) ! kappa (thermal cond.) + + WRITE(*,*)"chis(1) (chi_perp): ",chis(1) + WRITE(*,*)"chis(2) (chi_tor): ",chis(2) + WRITE(*,*)"chis(3) (kappa): ",chis(3) + + ! Use namelist kinetic inputs instead of equilibrium files + CALL params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, + $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + + ! Override desired normalized parameters + IF (ABS(inQ) > 0.0) THEN + Q = inQ ! NAMELIST + END IF + IF (ABS(inQ_e) > 0.0) THEN + Q_e = inQ_e ! NAMELIST + END IF + IF (ABS(inQ_i) > 0.0) THEN + Q_i = inQ_i ! NAMELIST + END IF + IF (inpr > 0.0) THEN + P_perp = inpr ! NAMELIST + END IF + IF (intau > 0.0) THEN + tau = intau ! NAMELIST + END IF + IF (inds > 0.0) THEN + D_norm = inds ! NAMELIST + END IF + + IF (.NOT. est_gamma_flag) THEN + CALL allocate_inputs(n_k,sl_in) + CALL allocate_outputs(n_k,sl_out) + END IF + + sl_in%qval_arr = (/ qval /) + sl_in%omegas_arr = (/ omega /) + !sl_in%Q_arr = (/ Q /) + sl_in%Q_e_arr = (/ Q_e /) + sl_in%Q_i_arr = (/ Q_i /) + sl_in%psi_n_arr = (/ 0.0 /) + sl_in%Re_dp_arr = (/ REAL(delta_prime) /) + sl_in%Im_dp_arr = (/ AIMAG(delta_prime) /) + sl_in%d_crit_arr = (/ dc_tmp /) + sl_in%P_perp_arr = (/ P_perp /) + sl_in%P_tor_arr = (/ P_tor /) + sl_in%tau_arr = (/ tau /) + sl_in%D_norm_arr = (/ D_norm /) + sl_in%d_beta_arr = (/ d_beta /) + sl_in%gammafac_arr = (/ gamma_fac /) + sl_in%c_beta_arr = (/ c_beta /) + sl_in%lu_arr = (/ lu /) + sl_in%Qconv_arr = (/ tauk /) + + END IF + END IF + +c----------------------------------------------------------------------- +c allocate output arrays for AMR delta storage. +c----------------------------------------------------------------------- + IF (AMR_flag .AND. .NOT. coupling_flag) THEN + ALLOCATE(all_deltas_out(n_k)) + ELSEIF (AMR_flag .AND. coupling_flag) THEN + ALLOCATE(all_deltas_out(1)) + END IF + + IF (AMR_flag) THEN + IF (coupling_flag) THEN + m_AMR = 1 + ELSE + m_AMR = MIN(n_k,msing_max) + END IF + END IF + + WRITE(*,*),"Rational q domain: ",sl_in%qval_arr +c----------------------------------------------------------------------- +c loop over rational surfaces to find matched growth rates. +c----------------------------------------------------------------------- + DO k=1,MIN(n_k,msing_max) + WRITE(*,*) + WRITE(*,'(A,I0,A)') 'Calculating growth rate on q = ', + $ sl_in%qval_arr(k),' rational surface:' + + ! Load per-surface parameters into module-level scalars. + Q_e = sl_in%Q_e_arr(k) + Q_i = sl_in%Q_i_arr(k) + P_perp = sl_in%P_perp_arr(k) + P_tor = sl_in%P_tor_arr(k) + tau = sl_in%tau_arr(k) + D_norm = sl_in%D_norm_arr(k) + c_beta = sl_in%c_beta_arr(k) + tauk = sl_in%Qconv_arr(k) + iota_e = Q_e / (Q_e - Q_i) + + WRITE(*,*)"Q_e: ",Q_e + WRITE(*,*)"Q_i: ",Q_i + WRITE(*,*)"P_perp: ",P_perp + WRITE(*,*)"P_tor: ",P_tor + WRITE(*,*)"tau: ",tau + WRITE(*,*)"D_norm: ",D_norm + WRITE(*,*)"tauk: ",tauk + WRITE(*,*)"iota_e: ",iota_e + WRITE(*,*)"Delta_prime: ",sl_in%Re_dp_arr(k) + WRITE(*,*)"Delta_crit: ",sl_in%d_crit_arr(k) + + ! Calculate (Deltaprime - D_crit)/S^1/3 + delta_eff = (sl_in%Re_dp_arr(k) - + $ sl_in%d_crit_arr(k))/(sl_in%lu_arr(k)**(1.0/3.0)) + pe = 0.0 + + ! Placeholder: gamma_sol_arr filled for external root-finder. + sl_out%gamma_sol_arr(k) = 0.0 + +c----------------------------------------------------------------------- +c uncoupled AMR scan (one surface at a time). +c dispersion_AMR_v2 (growthrates_mod) populates Q_store, D_store. +c----------------------------------------------------------------------- + IF (AMR_flag .AND. .NOT. coupling_flag) THEN + + WRITE(*,'(A,I0,A)') 'Calling uncoupled AMR scan on q = ', + $ sl_in%qval_arr(k),' rational surface:' + CALL dispersion_AMR_v2(n_k,sl_in,msing_max,scan_width, + $ Q_num,AMR_passes,coupling_flag) + + ! Re-allocate output arrays for this surface. + IF (ALLOCATED(all_deltas_out(k)%inQs)) + $ DEALLOCATE(all_deltas_out(k)%inQs) + IF (ALLOCATED(all_deltas_out(k)%iinQs)) + $ DEALLOCATE(all_deltas_out(k)%iinQs) + IF (ALLOCATED(all_deltas_out(k)%real_deltas)) + $ DEALLOCATE(all_deltas_out(k)%real_deltas) + IF (ALLOCATED(all_deltas_out(k)%imag_deltas)) + $ DEALLOCATE(all_deltas_out(k)%imag_deltas) + + ALLOCATE(all_deltas_out(k)%inQs(n_pts), + $ all_deltas_out(k)%iinQs(n_pts)) + ALLOCATE(all_deltas_out(k)%real_deltas(n_pts), + $ all_deltas_out(k)%imag_deltas(n_pts)) + + ! Flatten unique AMR points into 1-D output arrays. + DO i = 1, n_pts + all_deltas_out(k)%inQs(i) = REAL(Q_store(i)) + all_deltas_out(k)%iinQs(i) = -AIMAG(Q_store(i)) + + all_deltas_out(k)%real_deltas(i) = REAL(D_store(i)) + all_deltas_out(k)%imag_deltas(i) = AIMAG(D_store(i)) + END DO + + ! Clean up temporary AMR memory. + DEALLOCATE(Q_store, D_store) + + WRITE(*,'(A,I2,A,I7,A,2ES14.6)') + $ ' Surface', k, ': n_pts=', + $ SIZE(all_deltas_out(k)%real_deltas), + $ ' out_chksum=', + $ SUM(all_deltas_out(k)%real_deltas), + $ SUM(all_deltas_out(k)%imag_deltas) + + END IF + +c----------------------------------------------------------------------- +c single-surface stability scan on [Re(Q), Im(Q)] grid. +c Uses riccati_f() (new Fitzpatrick TJ-like formalism) or +c. riccati() (orig. SLAYER/Waelbroeck). +c----------------------------------------------------------------------- + IF ((stabscan_flag)) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,'(A,F0.1)')' >>> Running [Re(Q),'// + $ 'Im(Q)] scan with Q width = ', + $ scan_width + + ing_step = (2.0*scan_width) / (Q_num - 1) + count = 0 + + ALLOCATE(inQs(1:(Q_num+1)),iinQs(1:Q_num)) + ALLOCATE(deltas(1:(Q_num+1),1:Q_num)) + + DO i = 1, (Q_num+1) + DO j = 1, Q_num + ing_coarse = -scan_width + (i - 1) * ing_step ! added 0.5* + iing_coarse = -scan_width + (j - 1) * ing_step + ! Evaluate riccati function + g_tmp = CMPLX(ing_coarse,iing_coarse) + IF (Pperp_Ptor_flag) THEN + delta=riccati_f() + ELSE + delta=riccati(iing_coarse,Q_e,Q_i,P_perp, + $ c_beta,D_norm,tau,pe, + $ iinQ=ing_coarse) + END IF + inQs(i) = ing_coarse + IF (Pperp_Ptor_flag) THEN + iinQs(j) = iing_coarse + ELSE + iinQs(j) = -iing_coarse + END IF + deltas(i,j) = delta + ENDDO + ENDDO + + + IF (k<10) THEN + WRITE(UNIT=sm_str,FMT='(I1)') sl_in%qval_arr(k) + sm_str=ADJUSTL(sm_str) + ELSE + WRITE(UNIT=sm_str,FMT='(I2)') sl_in%qval_arr(k) + ENDIF + + OPEN(UNIT=out_unit,FILE="slayer_stability_n"// + $ TRIM(sn_str)//"m"//TRIM(sm_str)//".out", + $ STATUS="UNKNOWN") + WRITE(out_unit,'(1x,4(a17))') "RE(Q)", + $ "IM(Q)","RE(delta)","IM(delta)" + DO i=1,Q_num+1 + DO j=1,Q_num + WRITE(out_unit,'(1x,4(es17.8e3))') + $ inQs(i),iinQs(j), + $ REAL(deltas(i,j)),AIMAG(deltas(i,j)) + ENDDO + ENDDO + CLOSE(out_unit) + + DEALLOCATE(inQs,iinQs,deltas) + ENDIF + ENDDO + + IF (.NOT. (est_gamma_flag)) THEN + sl_in%d_beta_arr = (/ 0. /) + sl_out%dels_db_arr = (/ 0. /) + END IF + +c----------------------------------------------------------------------- +c coupled AMR scan (all surfaces simultaneously). +c dispersion_AMR_v2 (growthrates_mod) with coupling_flag = .TRUE. +c----------------------------------------------------------------------- + IF (AMR_flag .AND. coupling_flag) THEN + + CALL dispersion_AMR_v2(n_k,sl_in,msing_max,scan_width, + $ Q_num,AMR_passes,coupling_flag) + + ! Re-allocate output arrays. + IF (ALLOCATED(all_deltas_out(1)%inQs)) + $ DEALLOCATE(all_deltas_out(1)%inQs) + IF (ALLOCATED(all_deltas_out(1)%iinQs)) + $ DEALLOCATE(all_deltas_out(1)%iinQs) + IF (ALLOCATED(all_deltas_out(1)%real_deltas)) + $ DEALLOCATE(all_deltas_out(1)%real_deltas) + IF (ALLOCATED(all_deltas_out(1)%imag_deltas)) + $ DEALLOCATE(all_deltas_out(1)%imag_deltas) + + ALLOCATE(all_deltas_out(1)%inQs(n_pts), + $ all_deltas_out(1)%iinQs(n_pts)) + ALLOCATE(all_deltas_out(1)%real_deltas(n_pts), + $ all_deltas_out(1)%imag_deltas(n_pts)) + + ! Flatten unique AMR points into 1-D output arrays. + DO i = 1, n_pts + all_deltas_out(1)%inQs(i) = REAL(Q_store(i)) + + all_deltas_out(1)%iinQs(i) = -AIMAG(Q_store(i)) ! verify this sign convention + + all_deltas_out(1)%real_deltas(i) = REAL(D_store(i)) + all_deltas_out(1)%imag_deltas(i) = AIMAG(D_store(i)) + END DO + + ! Clean up temporary AMR memory. + DEALLOCATE(Q_store, D_store) + + END IF + +c----------------------------------------------------------------------- +c coupled-surface stability scan on [Re(Q), Im(Q)] grid. +c Uses dispersion_det (growthrates_mod) for the full dispersion +c determinant including inter-surface coupling. +c----------------------------------------------------------------------- + IF (coupled_stabscan_flag) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,'(A,F0.1)')' >>> Running [Re(Q),'// + $ 'Im(Q)] determinant scan with radius = ', + $ scan_width + + ing_step = (2.0 * scan_width) / (Q_num - 1) + count = 0 + + !IF (.NOT. stabscan_flag) THEN + ALLOCATE(inQs(1:(Q_num+1)),iinQs(1:Q_num)) + ALLOCATE(deltas(1:(Q_num+1),1:Q_num)) + !END IF + + inQs=0.0; iinQs=0.0; deltas=(0.0,0.0) + + DO i = 1, (Q_num+1) + DO j = 1, Q_num + ing_coarse = -scan_width + (i - 1) * ing_step + iing_coarse = -scan_width + (j - 1) * ing_step + inQs(i) = ing_coarse + iinQs(j) = iing_coarse + + ! Evaluate determinant + g_tmp = CMPLX(ing_coarse,iing_coarse) + deltas(i,j)=dispersion_det(g_tmp,n_k,sl_in,msing_max) + ENDDO + ENDDO + + OPEN(UNIT=out_unit,FILE="slayer_determinants_n"// + $ TRIM(sn_str)//".out", STATUS="UNKNOWN") + WRITE(out_unit,'(1x,4(a17))') "RE(Q)", + $ "IM(Q)","RE(det)","IM(det)" + DO i=1,Q_num+1 + DO j=1,Q_num + WRITE(out_unit,'(1x,4(es17.8e3))') + $ inQs(i),iinQs(j), + $ REAL(deltas(i,j)),AIMAG(deltas(i,j)) + ENDDO + ENDDO + CLOSE(out_unit) + + DEALLOCATE(inQs,iinQs,deltas) + END IF + + CALL output_gamma(est_gamma_flag,m_AMR,sl_in,sl_out, + $ all_deltas_out) + stop + ENDIF ! match_gamma_flag +c----------------------------------------------------------------------- +c Br threshold test scan (analytic). +c For testing & verification only. Scans rotation to find the +c critical radial-field threshold from a simple torque balance. +c----------------------------------------------------------------------- + IF (br_th_flag) THEN + WRITE(*,*)"------------------------------------------" + WRITE(*,*)">>> Computing Br threshold" + + IF (read_eq) THEN + sl_in%chi_p_arr = chi_p_prof + sl_in%chi_t_arr = chi_t_prof + sl_in%kappa_arr = kappa_prof + CALL build_inputs(infile,ncfile,sl_in) + n_k = SIZE(sl_in%qval_arr) + CALL allocate_outputs(n_k,sl_out) + ELSE + n_k = 1 + + chis(1) = chi_p_prof(1) + chis(2) = chi_t_prof(1) + chis(3) = kappa_prof(1) + + CALL params(n_e,t_e,t_i,omega,chis,dr_val,dgeo_val, + $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + + inQ=Q + inQ_e=Q_e + inQ_i=Q_i + inc_beta=c_beta + inds=ds + intau=tau + + CALL allocate_inputs(n_k,sl_in) + CALL allocate_outputs(n_k,sl_out) + + sl_in%qval_arr = (/ qval /) + sl_in%omegas_arr = (/ omega /) + sl_in%Q_e_arr = (/ Q_e /) + sl_in%Q_i_arr = (/ Q_i /) + sl_in%psi_n_arr = (/ 0.0 /) + sl_in%Re_dp_arr = (/ 0.0 /) + sl_in%Im_dp_arr = (/ 0.0 /) + sl_in%d_crit_arr = (/ 0.0 /) + sl_in%P_perp_arr = (/ P_perp /) + sl_in%P_tor_arr = (/ P_tor /) + sl_in%tau_arr = (/ tau /) + sl_in%D_norm_arr = (/ D_norm /) + sl_in%d_beta_arr = (/ d_beta /) + sl_in%gammafac_arr = (/ gamma_fac /) + sl_in%c_beta_arr = (/ c_beta /) + sl_in%lu_arr = (/ lu /) + sl_in%Qconv_arr = (/ tauk /) + END IF +c----------------------------------------------------------------------- +c loop over rational surfaces to compute Br threshold. +c----------------------------------------------------------------------- + delta_n_p = 1e-2 + inum = 200 + inQ_max = 10.0 + inQ_min = -10.0 + + DO k=1,n_k + WRITE(*,*) + WRITE(*,'(A,I0,A)') 'Computing Br threshold on q = ', + $ sl_in%qval_arr(k),' rational surface' + + Q_e = sl_in%Q_e_arr(k) + Q_i = sl_in%Q_i_arr(k) + inQ_e = Q_e + inQ_i = Q_i + inpr = sl_in%P_perp_arr(k) + inc_beta = sl_in%c_beta_arr(k) + inds = sl_in%D_norm_arr(k) + intau = sl_in%tau_arr(k) + Q0 = sl_in%Q_e_arr(k) + + ALLOCATE(inQs(0:inum),deltal(0:inum), + $ jxbl(0:inum),bal(0:inum)) + DO i=0,inum + inQs(i)=inQ_min+(REAL(i)/inum)*(inQ_max-inQ_min) + deltal(i)=riccati(inQs(i),inQ_e,inQ_i, + $ inpr,inc_beta,inds,intau,inpe) + jxbl(i)=-AIMAG(1.0/(deltal(i)+delta_n_p)) + bal(i)=2.0*inpr*(Q0-inQs(i))/jxbl(i) + ENDDO + + iloc=MAXLOC(bal) + Q_sol=inQs(iloc(1)) + br_th=SQRT(MAXVAL(bal)/sl_in%lu_arr(k) + $ *(sval**2.0/2.0)) + + sl_out%br_th_arr(k) = br_th + sl_out%gamma_sol_arr(k) = 0.0 + sl_out%gamma_est_arr(k) = 0.0 + sl_out%dels_db_arr(k) = 0.0 + + WRITE(*,'(A,ES12.4)') ' br_th = ', br_th + DEALLOCATE(inQs,deltal,jxbl,bal) + ENDDO +c----------------------------------------------------------------------- +c write output. +c----------------------------------------------------------------------- + CALL output_gamma(est_gamma_flag,m_AMR,sl_in,sl_out, + $ all_deltas_out) + STOP + ENDIF ! br_th_flag +c----------------------------------------------------------------------- c find solutions based on simple torque balance. c----------------------------------------------------------------------- IF (bal_flag)THEN @@ -288,7 +1111,7 @@ PROGRAM slayer ! write components of torque balance IF(ascii_flag)THEN OPEN(UNIT=out_unit,FILE="slayer_bal_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,5(a17))') "inQ","RE(delta)", $ "IM(delta)","jxb","bal" @@ -299,8 +1122,8 @@ PROGRAM slayer CLOSE(out_unit) ENDIF - index=MAXLOC(bal) - Q_sol=inQs(index(1)) + iloc=MAXLOC(bal) + Q_sol=inQs(iloc(1)) br_th=sqrt(MAXVAL(bal)/lu*(sval**2.0/2.0))*1e4 WRITE(*,*)"Q_sol=",Q_sol WRITE(*,*)"br_th=",br_th @@ -324,10 +1147,10 @@ PROGRAM slayer $ inds,intau,inpe,iinQ=iinQs(j)) ENDDO ENDDO - + IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_stability_n"// - $ TRIM(sn)//".out", STATUS="UNKNOWN") + $ TRIM(sn_str)//".out", STATUS="UNKNOWN") WRITE(out_unit,'(1x,4(a17))') "RE(Q)", $ "IM(Q)","RE(delta)","IM(delta)" DO i=0,inum @@ -364,7 +1187,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_riccatiscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,5(a17))') "x","yphs","yamp", $ "RE(delta)","IM(delta)" DO j=0,jnum @@ -406,7 +1229,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_QPescan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Q","Pe","RE(delta)", $ "IM(delta)","psi","jxb" DO j=0,jnum @@ -421,7 +1244,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_QPescan_n' - $ //TRIM(sn)//'.bin', + $ //TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -461,7 +1284,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_QPscan"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Q","Pr","RE(delta)", $ "IM(delta)","psi","jxb" DO j=0,jnum @@ -476,7 +1299,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE="slayer_QPscan_" - $ //TRIM(sn)//".bin", + $ //TRIM(sn_str)//".bin", $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -510,7 +1333,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_Qscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Q","Pr","RE(delta)", $ "IM(delta)","psi","jxb" DO j=0,jnum @@ -552,7 +1375,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_QPscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Q","Pr","RE(delta)", $ "IM(delta)","psi","jxb" DO j=0,jnum @@ -567,7 +1390,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE="slayer_QPscan_n"// - $ TRIM(sn)//".bin", + $ TRIM(sn_str)//".bin", $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -620,7 +1443,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_QDscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -651,8 +1474,9 @@ PROGRAM slayer DO k=0,knum ks(j,k)=k_min+(k_max-k_min)*(REAL(k)/knum) - CALL params(n_e*ks(j,k),t_e,t_i,omega*js(j,k), - $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + CALL params(n_e*ks(j,k),t_e,t_i,omega*js(j,k),chis, + $ dr_val,dgeo_val,l_n,l_t,qval,sval,bt,rs, + $ R0,mu_i,zeff,params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -683,8 +1507,8 @@ PROGRAM slayer jxb=-AIMAG(1.0/(delta+delta_n_p)) bal(i)=2.0*inpr*(Q0-inQs(i))/jxb ENDDO - index=MAXLOC(bal) - Q_sols(j,k)=inQs(index(1)) + iloc=MAXLOC(bal) + Q_sols(j,k)=inQs(iloc(1)) br_ths(j,k)=sqrt(MAXVAL(bal)/lu)*1e4 WRITE(*,*)"br_ths=",br_ths(j,k) ENDDO @@ -692,7 +1516,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_onscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Omega","Density", $ "Omega_i","Omega_e","Omega_sol","Field_Threshold" DO j=0,jnum @@ -709,7 +1533,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_onscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -739,8 +1563,8 @@ PROGRAM slayer ks(j,k)=k_min+(k_max-k_min)*(REAL(k)/knum) CALL params(n_e,t_e*ks(j,k),t_i*ks(j,k), - $ omega*js(j,k),l_n,l_t,qval,sval,bt, - $ rs,R0,mu_i,zeff,params_check) + $ omega*js(j,k),chis,dr_val,dgeo_val,l_n,l_t,qval, + $ sval,bt,rs,R0,mu_i,zeff,params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -769,8 +1593,8 @@ PROGRAM slayer jxb=-AIMAG(1.0/(delta+delta_n_p)) bal(i)=2.0*inpr*(Q0-inQs(i))/jxb ENDDO - index=MAXLOC(bal) - Q_sols(j,k)=inQs(index(1)) + iloc=MAXLOC(bal) + Q_sols(j,k)=inQs(iloc(1)) br_ths(j,k)=sqrt(MAXVAL(bal)/lu)*1e4 WRITE(*,*)"t_e=",t_e*ks(j,k),"br_ths=",br_ths(j,k) ENDDO @@ -778,7 +1602,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_otscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Omega","Temperature", $ "Omega_i","Omega_e","Omega_sol","Field_Threshold" DO j=0,jnum @@ -795,7 +1619,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_otscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -825,7 +1649,8 @@ PROGRAM slayer ks(j,k)=k_min+(k_max-k_min)*(REAL(k)/knum) CALL params(n_e*ks(j,k),t_e*js(j,k),t_i*js(j,k),omega, - $ l_n,l_t,qval,sval,bt,rs,R0,mu_i,zeff,params_check) + $ chis,dr_val,dgeo_val,l_n,l_t,qval,sval,bt,rs,R0, + $ mu_i,zeff,params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -854,8 +1679,8 @@ PROGRAM slayer jxb=-AIMAG(1.0/(delta+delta_n_p)) bal(i)=2.0*inpr*(Q0-inQs(i))/jxb ENDDO - index=MAXLOC(bal) - Q_sols(j,k)=inQs(index(1)) + iloc=MAXLOC(bal) + Q_sols(j,k)=inQs(iloc(1)) br_ths(j,k)=sqrt(MAXVAL(bal)/lu)*1e4 WRITE(*,*)"br_ths=",br_ths(j,k) ENDDO @@ -863,7 +1688,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_ntscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,6(a17))') "Temperature","Density", $ "Omega_i","Omega_e","Omega_sol","Field_Threshold" DO j=0,jnum @@ -879,7 +1704,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_ntscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum @@ -909,9 +1734,9 @@ PROGRAM slayer ks(j,k)=k_min+(k_max-k_min)*(REAL(k)/knum) - CALL params(n_e*ks(j,k),t_e,t_i,omega, - $ l_n,l_t,qval,sval,bt*js(j,k),rs,R0,mu_i,zeff, - $ params_check) + CALL params(n_e*ks(j,k),t_e,t_i,omega,chis,dr_val, + $ dgeo_val,l_n,l_t,qval,sval,bt*js(j,k),rs,R0,mu_i, + $ zeff,params_check) inQ=Q inQ_e=Q_e inQ_i=Q_i @@ -941,8 +1766,8 @@ PROGRAM slayer jxb=-AIMAG(1.0/(delta+delta_n_p)) bal(i)=2.0*inpr*(Q0-inQs(i))/jxb ENDDO - index=MAXLOC(bal) - Q_sols(j,k)=inQs(index(1)) + iloc=MAXLOC(bal) + Q_sols(j,k)=inQs(iloc(1)) br_ths(j,k)=sqrt(MAXVAL(bal)/lu)*1e4 WRITE(*,*)"br_ths=",br_ths(j,k) ENDDO @@ -950,7 +1775,7 @@ PROGRAM slayer IF (ascii_flag) THEN OPEN(UNIT=out_unit,FILE="slayer_nbtscan_n"// - $ TRIM(sn)//".out",STATUS="UNKNOWN") + $ TRIM(sn_str)//".out",STATUS="UNKNOWN") WRITE(out_unit,'(1x,4(a17))') "Bt","Density", $ "Omega_sol","Field_Threshold" DO j=0,jnum @@ -966,7 +1791,7 @@ PROGRAM slayer IF (bin_flag) THEN OPEN(UNIT=bin_2d_unit,FILE='slayer_nbtscan_n'// - $ TRIM(sn)//'.bin', + $ TRIM(sn_str)//'.bin', $ STATUS='UNKNOWN',POSITION='REWIND',FORM='UNFORMATTED') WRITE(bin_2d_unit)1,0 WRITE(bin_2d_unit)jnum,knum diff --git a/slayer/slayer_netcdf.f b/slayer/slayer_netcdf.f new file mode 100644 index 00000000..82e9bd91 --- /dev/null +++ b/slayer/slayer_netcdf.f @@ -0,0 +1,408 @@ +c======================================================================= +c file slayer_netcdf.f +c Writes SLAYER solver results to a NetCDF output file. +c +c The per-surface scalar inputs (Lundquist number, Q-normalisation, +c Prandtl numbers, …) and the solver outputs (growth rates, Delta +c values) are stored in a single NetCDF-3/64-bit-offset file named +c slayer_output_n.nc. +c +c Ragged AMR scan data (variable number of evaluation points per +c surface) are zero-padded into rectangular arrays before writing. +c +c (Resolved: FLAG 1 -- early RETURN for msing == 0. +c FLAG 2 -- buffers now use fill_val; _FillValue attributes added. +c FLAG 3 -- Q_id removed (unused). +c FLAG 4 -- c_b_id removed (unused). +c FLAG 5 -- stale unused-variable comment removed. +c FLAG 6 -- version from INCLUDE "version.inc". +c FLAG 7 -- local sn_local replaces global sn_str.) +c======================================================================= +c----------------------------------------------------------------------- +c code organisation. +c----------------------------------------------------------------------- +c 0. slayer_netcdf_mod -- module declarations +c 1. sl_check -- NetCDF status checker +c 2. slayer_netcdf_out -- main output routine +c----------------------------------------------------------------------- +c +c----------------------------------------------------------------------- +c subprogram 0. slayer_netcdf_mod. +c Module wrapper — imports sglobal_mod (shared types and globals) +c and the NetCDF Fortran-90 API. +c----------------------------------------------------------------------- + MODULE slayer_netcdf_mod + + USE sglobal_mod + USE netcdf + USE ieee_arithmetic, ONLY: ieee_value, ieee_quiet_nan + + IMPLICIT NONE + + CONTAINS +c +c----------------------------------------------------------------------- +c subprogram 1. sl_check. +c Assert that a NetCDF operation succeeded; abort with a message +c if it did not. +c----------------------------------------------------------------------- + SUBROUTINE sl_check(stat) +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + INTEGER, INTENT(IN) :: stat ! return code from any nf90_* call +c----------------------------------------------------------------------- +c check status and abort on error. +c----------------------------------------------------------------------- + IF (stat /= nf90_noerr) THEN + PRINT *, TRIM(nf90_strerror(stat)) + STOP "ERROR: failed to write/read netcdf file" + ENDIF + + RETURN + END SUBROUTINE sl_check +c +c----------------------------------------------------------------------- +c subprogram 2. slayer_netcdf_out. +c Write per-surface SLAYER inputs and solver outputs to the +c NetCDF file slayer_output_n.nc . +c +c Arguments: +c msing -- number of rational surfaces +c m_AMR -- number of AMR-scanned surfaces +c est_gamma_flag -- .TRUE. to include estimated growth rates +c sl_in -- slayer_inputs_type (per-surface inputs) +c sl_out -- slayer_outputs_type (solver results) +c all_deltas_out -- array(m_AMR) of deltas_outputs_type (AMR +c scan results, potentially ragged) +c----------------------------------------------------------------------- + SUBROUTINE slayer_netcdf_out(msing, m_AMR, est_gamma_flag, + $ sl_in, sl_out, all_deltas_out) +c----------------------------------------------------------------------- +c declarations -- subroutine arguments. +c----------------------------------------------------------------------- + INTEGER, INTENT(IN) :: msing ! number of rational surfaces + INTEGER, INTENT(IN) :: m_AMR ! number of AMR surfaces + LOGICAL, INTENT(IN) :: est_gamma_flag ! include estimated gammas? + + TYPE(slayer_inputs_type), INTENT(IN) :: sl_in + TYPE(slayer_outputs_type), INTENT(IN) :: sl_out + TYPE(deltas_outputs_type), INTENT(IN) :: all_deltas_out(m_AMR) + +c----------------------------------------------------------------------- +c declarations -- NetCDF file and dimension IDs. +c----------------------------------------------------------------------- + INTEGER :: ncid ! NetCDF file ID + INTEGER :: qsing_dim ! dim: rational surfaces (msing) + INTEGER :: nAMR_dim ! dim: AMR surfaces (m_AMR) + INTEGER :: i_dim ! dim: Re/Im component (2) + INTEGER :: dim_pts_id ! dim: max AMR eval points + +c----------------------------------------------------------------------- +c declarations -- NetCDF variable IDs. +c Each *_id holds the handle returned by nf90_def_var and is +c later passed to the matching nf90_put_var call. +c----------------------------------------------------------------------- + INTEGER :: qsing_id ! "r" — surface index + INTEGER :: qr_id ! "q_rational" — safety factor + INTEGER :: omegas_id ! "omegas" — rotation freq + INTEGER :: qc_id ! "tau_k" — Q-conversion + INTEGER :: Q_e_id ! "Q_e" — norm Q_e + INTEGER :: Q_i_id ! "Q_i" — norm Q_i + INTEGER :: S_id ! "S" — Lundquist + INTEGER :: pr_id ! "psi_n_rational" — norm psi + INTEGER :: p_perp_id ! "P_perp" + INTEGER :: p_tor_id ! "P_tor" + INTEGER :: Dnorm_id ! "D" — normalised D + INTEGER :: dpp_id ! "Delta_prime_rational" (complex) + INTEGER :: dc_id ! "Delta_crit_rational" + INTEGER :: dels_db_id ! "delta_s_d_b" (complex) + INTEGER :: d_b_id ! "d_beta" + INTEGER :: gs_id ! "growth rate" (complex) + INTEGER :: ge_id ! "est. growth rate" (complex) + INTEGER :: br_th_id ! "br_th" — Br threshold + +c AMR variable IDs + INTEGER :: var_q_id ! "Q_AMR" — scan Q-points + INTEGER :: var_d_id ! "Deltas_AMR" — scan Delta values + INTEGER :: var_npts_id ! "n_amr_pts" — points per surface + +c----------------------------------------------------------------------- +c declarations -- AMR rectangular-buffer workspace. +c The ragged per-surface scan data are padded into fixed-size +c rectangular arrays before writing to NetCDF. +c----------------------------------------------------------------------- + INTEGER :: max_pts_all ! max points across surfaces + INTEGER :: s ! surface loop index + INTEGER :: n_curr ! points on current surface + REAL(r8), ALLOCATABLE :: buffer_q(:,:,:) ! (pts, surf, Re/Im) + REAL(r8), ALLOCATABLE :: buffer_d(:,:,:) ! (pts, surf, Re/Im) + INTEGER, ALLOCATABLE :: n_pts_arr(:) ! points per surface + REAL(r8) :: fill_val ! NaN padding for ragged arrays + +c----------------------------------------------------------------------- +c declarations -- miscellaneous locals. +c----------------------------------------------------------------------- + CHARACTER(64) :: ncfile ! output file name + CHARACTER(2) :: sn_local ! local n-string for filename + LOGICAL, PARAMETER :: debug_flag = .FALSE. ! verbose trace + INCLUDE "version.inc" + +c----------------------------------------------------------------------- +c build the output filename from the toroidal mode number. +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, "Called slayer_netcdf_out" + + IF (nn < 10) THEN + WRITE(UNIT=sn_local, FMT='(I1)') nn + sn_local = ADJUSTL(sn_local) + ELSE + WRITE(UNIT=sn_local, FMT='(I2)') nn + ENDIF + ncfile = "slayer_output_n"//TRIM(sn_local)//".nc" + IF (debug_flag) PRINT *, ncfile + +c----------------------------------------------------------------------- +c create the NetCDF file (clobber any existing file). +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, " - Creating netcdf file" + CALL sl_check( nf90_create(ncfile, + $ cmode=OR(NF90_CLOBBER, NF90_64BIT_OFFSET), ncid=ncid) ) + +c----------------------------------------------------------------------- +c reform ragged AMR Delta outputs into rectangular arrays. +c +c Each surface may have a different number of AMR scan points. +c We find the maximum, allocate rectangular buffers of that size, +c and copy in the per-surface data. Unused trailing slots are +c filled with fill_val (-9.99E33). +c----------------------------------------------------------------------- + +c step 1: find the maximum AMR grid size across all surfaces. + max_pts_all = 0 + IF (ALLOCATED(all_deltas_out(1)%inQs)) THEN + DO s = 1, m_AMR + max_pts_all = MAX(max_pts_all, + $ SIZE(all_deltas_out(s)%inQs)) + END DO + END IF + +c step 2: allocate rectangular buffers (pts × surfaces × Re/Im). + ALLOCATE(buffer_q(max_pts_all, m_AMR, 2)) + ALLOCATE(buffer_d(max_pts_all, m_AMR, 2)) + ALLOCATE(n_pts_arr(m_AMR)) + + fill_val = ieee_value(1.0d0, ieee_quiet_nan) + buffer_q = fill_val ! NaN padding for ragged arrays + buffer_d = fill_val + n_pts_arr = 0 + +c step 3: flatten the ragged data into the buffers. + IF (ALLOCATED(all_deltas_out(1)%inQs)) THEN + DO s = 1, m_AMR + n_curr = SIZE(all_deltas_out(s)%inQs) + n_pts_arr(s) = n_curr + +c Q-coordinate (Re and Im parts) + buffer_q(1:n_curr, s, 1) = + $ all_deltas_out(s)%inQs(1:n_curr) + buffer_q(1:n_curr, s, 2) = + $ all_deltas_out(s)%iinQs(1:n_curr) + +c Delta result (Re and Im parts) + buffer_d(1:n_curr, s, 1) = + $ all_deltas_out(s)%real_deltas(1:n_curr) + buffer_d(1:n_curr, s, 2) = + $ all_deltas_out(s)%imag_deltas(1:n_curr) + END DO + END IF + +c----------------------------------------------------------------------- +c define global file attributes. +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, " - Defining netcdf globals" + CALL sl_check( nf90_put_att(ncid, nf90_global, + $ "title", "SLAYER outputs") ) + CALL sl_check( nf90_put_att(ncid, nf90_global, + $ "version", version) ) + +c----------------------------------------------------------------------- +c define dimensions and per-surface NetCDF variables. +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, " - Defining dimensions in netcdf" + WRITE(*,*) ">>> Writing results to NetCDF output file" + + IF (msing == 0) THEN + WRITE(*,*) "WARNING: msing == 0, skipping NetCDF output" + DEALLOCATE(buffer_q, buffer_d, n_pts_arr) + CALL sl_check( nf90_close(ncid) ) + RETURN + END IF + +c -- core dimensions -- + CALL sl_check( nf90_def_dim(ncid, "r", msing, qsing_dim) ) + CALL sl_check( nf90_def_dim(ncid, "r_AMR", m_AMR, nAMR_dim) ) + CALL sl_check( nf90_def_dim(ncid, "i", 2, i_dim) ) + +c -- scalar per-surface variables -- + CALL sl_check( nf90_def_var(ncid, "r", nf90_int, + $ qsing_dim, qsing_id) ) + CALL sl_check( nf90_def_var(ncid, "q_rational", + $ nf90_double, qsing_dim, qr_id) ) + CALL sl_check( nf90_def_var(ncid, "omegas", nf90_double, + $ qsing_dim, omegas_id) ) + CALL sl_check( nf90_def_var(ncid, "tau_k", nf90_double, + $ qsing_dim, qc_id) ) + CALL sl_check( nf90_def_var(ncid, "Q_e", nf90_double, + $ qsing_dim, Q_e_id) ) + CALL sl_check( nf90_def_var(ncid, "Q_i", nf90_double, + $ qsing_dim, Q_i_id) ) + CALL sl_check( nf90_def_var(ncid, "S", nf90_double, + $ qsing_dim, S_id) ) + CALL sl_check( nf90_def_var(ncid, "psi_n_rational", + $ nf90_double, qsing_dim, pr_id) ) + CALL sl_check( nf90_def_var(ncid, "P_perp", nf90_double, + $ qsing_dim, p_perp_id) ) + CALL sl_check( nf90_def_var(ncid, "P_tor", nf90_double, + $ qsing_dim, p_tor_id) ) + +c----------------------------------------------------------------------- +c define additional variables (D, Delta', growth rates). +c----------------------------------------------------------------------- + CALL sl_check( nf90_def_var(ncid, "D", nf90_double, + $ qsing_dim, Dnorm_id) ) + CALL sl_check( nf90_def_var(ncid, "Delta_prime_rational", + $ nf90_double, (/qsing_dim, i_dim/), dpp_id) ) + CALL sl_check( nf90_def_var(ncid, "Delta_crit_rational", + $ nf90_double, qsing_dim, dc_id) ) + + IF (est_gamma_flag) THEN + CALL sl_check( nf90_def_var(ncid, "delta_s_d_b", + $ nf90_double, (/qsing_dim, i_dim/), dels_db_id) ) + CALL sl_check( nf90_def_var(ncid, "d_beta", nf90_double, + $ qsing_dim, d_b_id) ) + CALL sl_check( nf90_def_var(ncid, "est. growth rate", + $ nf90_double, (/qsing_dim, i_dim/), ge_id) ) + END IF + + CALL sl_check( nf90_def_var(ncid, "growth rate", + $ nf90_double, (/qsing_dim, i_dim/), gs_id) ) + + IF (ALLOCATED(sl_out%br_th_arr)) THEN + CALL sl_check( nf90_def_var(ncid, "br_th", + $ nf90_double, qsing_dim, br_th_id) ) + END IF + +c----------------------------------------------------------------------- +c define AMR scan dimensions and variables. +c----------------------------------------------------------------------- + CALL sl_check( nf90_def_dim(ncid, "amr_pts", + $ max_pts_all, dim_pts_id) ) + CALL sl_check( nf90_def_var(ncid, "n_amr_pts", NF90_INT, + $ (/nAMR_dim/), var_npts_id) ) + CALL sl_check( nf90_def_var(ncid, "Q_AMR", NF90_DOUBLE, + $ (/dim_pts_id, nAMR_dim, i_dim/), var_q_id) ) + CALL sl_check( nf90_def_var(ncid, "Deltas_AMR", NF90_DOUBLE, + $ (/dim_pts_id, nAMR_dim, i_dim/), var_d_id) ) + +c set _FillValue attribute on AMR arrays for proper padding. + CALL sl_check( nf90_put_att(ncid, var_q_id, + $ "_FillValue", fill_val) ) + CALL sl_check( nf90_put_att(ncid, var_d_id, + $ "_FillValue", fill_val) ) + +c----------------------------------------------------------------------- +c end NetCDF define mode. +c----------------------------------------------------------------------- + CALL sl_check( nf90_enddef(ncid) ) + +c----------------------------------------------------------------------- +c write per-surface scalar variables. +c----------------------------------------------------------------------- + CALL sl_check( nf90_put_var(ncid, qsing_id, + $ sl_in%qval_arr) ) + CALL sl_check( nf90_put_var(ncid, qr_id, + $ sl_in%qval_arr) ) + CALL sl_check( nf90_put_var(ncid, pr_id, + $ sl_in%psi_n_arr) ) + CALL sl_check( nf90_put_var(ncid, omegas_id, + $ sl_in%omegas_arr)) + CALL sl_check( nf90_put_var(ncid, S_id, + $ sl_in%lu_arr) ) + CALL sl_check( nf90_put_var(ncid, qc_id, + $ sl_in%Qconv_arr) ) + CALL sl_check( nf90_put_var(ncid, Q_e_id, + $ sl_in%Q_e_arr) ) + CALL sl_check( nf90_put_var(ncid, Q_i_id, + $ sl_in%Q_i_arr) ) + CALL sl_check( nf90_put_var(ncid, p_perp_id, + $ sl_in%P_perp_arr)) + CALL sl_check( nf90_put_var(ncid, p_tor_id, + $ sl_in%P_tor_arr) ) + CALL sl_check( nf90_put_var(ncid, Dnorm_id, + $ sl_in%D_norm_arr)) + +c----------------------------------------------------------------------- +c write complex Delta' as a RESHAPE'd (msing, 2) real array. +c----------------------------------------------------------------------- + CALL sl_check( nf90_put_var(ncid, dpp_id, + $ RESHAPE( (/sl_in%Re_dp_arr, sl_in%Im_dp_arr/), + $ (/msing, 2/) )) ) + CALL sl_check( nf90_put_var(ncid, dc_id, + $ sl_in%d_crit_arr) ) + +c----------------------------------------------------------------------- +c write estimated growth-rate outputs (only when requested). +c----------------------------------------------------------------------- + IF (est_gamma_flag) THEN + CALL sl_check( nf90_put_var(ncid, dels_db_id, + $ RESHAPE( (/REAL(sl_out%dels_db_arr), + $ AIMAG(sl_out%dels_db_arr)/), + $ (/msing, 2/) )) ) + + CALL sl_check( nf90_put_var(ncid, d_b_id, + $ sl_in%d_beta_arr) ) + + CALL sl_check( nf90_put_var(ncid, ge_id, + $ RESHAPE( (/REAL(sl_out%gamma_est_arr), + $ AIMAG(sl_out%gamma_est_arr)/), + $ (/msing, 2/) )) ) + END IF + +c----------------------------------------------------------------------- +c write solved growth rate (always present). +c----------------------------------------------------------------------- + CALL sl_check( nf90_put_var(ncid, gs_id, + $ RESHAPE( (/REAL(sl_out%gamma_sol_arr), + $ AIMAG(sl_out%gamma_sol_arr)/), + $ (/msing, 2/) )) ) + + IF (ALLOCATED(sl_out%br_th_arr)) THEN + CALL sl_check( nf90_put_var(ncid, br_th_id, + $ sl_out%br_th_arr) ) + END IF + +c----------------------------------------------------------------------- +c write AMR scan arrays. +c----------------------------------------------------------------------- + CALL sl_check( nf90_put_var(ncid, var_npts_id, n_pts_arr) ) + CALL sl_check( nf90_put_var(ncid, var_q_id, buffer_q) ) + CALL sl_check( nf90_put_var(ncid, var_d_id, buffer_d) ) + +c----------------------------------------------------------------------- +c deallocate AMR rectangular buffers. +c----------------------------------------------------------------------- + DEALLOCATE(buffer_q, buffer_d, n_pts_arr) + +c----------------------------------------------------------------------- +c close the NetCDF file. +c----------------------------------------------------------------------- + IF (debug_flag) PRINT *, " - Closing netcdf file" + CALL sl_check( nf90_close(ncid) ) + +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE slayer_netcdf_out + END MODULE slayer_netcdf_mod \ No newline at end of file diff --git a/stride/free.f b/stride/free.f index c7cd1f34..043007b1 100755 --- a/stride/free.f +++ b/stride/free.f @@ -166,7 +166,6 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,op_netcdf_out) INTEGER :: ipert,jpert,isol,info,lwork INTEGER, DIMENSION(mpert) :: m INTEGER, DIMENSION(1) :: imax - REAL(r8) :: v1 REAL(r8), DIMENSION(mpert) :: ep,ev,et REAL(r8), DIMENSION(3*mpert-1) :: rwork COMPLEX(r8) :: phase,norm @@ -174,6 +173,15 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,op_netcdf_out) COMPLEX(r8), DIMENSION(mpert,mpert) :: wt,wpt,wvt COMPLEX(r8), DIMENSION(mpert,mpert) :: nmat,smat CHARACTER(24), DIMENSION(mpert) :: message + INTEGER :: ipsi,itheta + REAL(r8), DIMENSION(sq%mx+1) :: ln_q + REAL(r8), DIMENSION(msing) :: dgeo,shr + TYPE(spline_type) :: psi_t,avg_dpsi_spl,avg_bsq_spl,v_spl, + $ vol_spl,shr_spl + REAL(r8) :: bsq,chi1,dpsisq,myeta,jac,psifac,q,q1,respsi, + $ rfac,v1,v21,v22,v23,v33,al,Lam,mytheta,myr,V_s + REAL(r8), DIMENSION(:), POINTER :: avg + TYPE(spline_type), TARGET :: fspl c----------------------------------------------------------------------- c write formats. c----------------------------------------------------------------------- @@ -293,15 +301,155 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,op_netcdf_out) WRITE(out_unit,90)(isol,ep(isol),ev(isol),isol=1,mpert) WRITE(out_unit,80) c----------------------------------------------------------------------- +c compute toroidal Delta_crit +c----------------------------------------------------------------------- + ! Prepare toroidal flux spline + CALL spline_alloc(psi_t,SIZE(sq%fsi(:, 4))-1,1) + psi_t%xs=sq%xs(:) + psi_t%fs(:,1)=sq%fsi(:,4)*twopi*psio ! Un-normalize toroidal flux + CALL spline_fit(psi_t,"extrap") + + ! Prepare geometric splines + CALL spline_alloc(avg_dpsi_spl,SIZE(sq%xs(:))-1,1) + avg_dpsi_spl%xs=sq%xs(:) + CALL spline_alloc(avg_bsq_spl,SIZE(sq%xs(:))-1,1) + avg_bsq_spl%xs=sq%xs(:) + CALL spline_alloc(v_spl,SIZE(sq%xs(:))-1,1) + v_spl%xs=sq%xs(:) + ! Cumulative volume V(psi_N) [m^3] = integral of dV/dpsi from axis. + CALL spline_alloc(vol_spl,SIZE(sq%xs(:))-1,1) + vol_spl%xs=sq%xs(:) + vol_spl%fs(:,1)=sq%fsi(:,3) + CALL spline_fit(vol_spl,"extrap") + + ! Prepare shear spline + CALL spline_alloc(shr_spl,SIZE(sq%xs(:))-1,1) + shr_spl%xs=sq%xs(:) + ln_q=LOG(sq%fs(:,4)) + shr_spl%fs(:,1)=ln_q + + CALL spline_fit(shr_spl,"extrap") + + CALL spline_alloc(fspl,mtheta,3) + fspl%xs=rzphi%ys + + DO ipsi=0,mpsi + psifac=sq%xs(ipsi) + v1=sq%fs(ipsi,3) + q=sq%fs(ipsi,4) + q1=sq%fs1(ipsi,4) + chi1=twopi*psio +c----------------------------------------------------------------------- +c evaluate coordinates and jacobian. +c----------------------------------------------------------------------- + DO itheta=0,mtheta + CALL bicube_eval(rzphi,rzphi%xs(ipsi),rzphi%ys(itheta),1) + mytheta=rzphi%ys(itheta) + rfac=SQRT(rzphi%f(1)) + myeta=twopi*(mytheta+rzphi%f(2)) + myr=ro+rfac*COS(myeta) + jac=rzphi%f(4) +c----------------------------------------------------------------------- +c evaluate other local quantities. +c----------------------------------------------------------------------- + v21=rzphi%fy(1)/(2*rfac*jac) + v22=(1+rzphi%fy(2))*twopi*rfac/jac + v23=rzphi%fy(3)*myr/jac + v33=twopi*myr/jac + bsq=chi1**2*(v21**2+v22**2+(v23+q*v33)**2) + dpsisq=(twopi*myr)**2*(v21**2+v22**2) +c----------------------------------------------------------------------- +c evaluate integrands. +c----------------------------------------------------------------------- + fspl%fs(itheta,1)=dpsisq*(v1**2) ! Converting to \nabla V + fspl%fs(itheta,2)=bsq + fspl%fs(itheta,3)=v1 + fspl%fs(itheta,:)=fspl%fs(itheta,:)*(jac/v1) + ENDDO +c----------------------------------------------------------------------- +c integrate quantities with respect to theta. +c----------------------------------------------------------------------- + CALL spline_fit(fspl,"periodic") + CALL spline_int(fspl) + avg => fspl%fsi(mtheta,:) + avg_dpsi_spl%fs(ipsi,1)=avg(1) + avg_bsq_spl%fs(ipsi,1)=avg(2) + v_spl%fs(ipsi,1)=avg(3) + ENDDO + CALL spline_dealloc(fspl) + + CALL spline_fit(avg_dpsi_spl,"extrap") + CALL spline_fit(avg_bsq_spl,"extrap") + CALL spline_fit(v_spl,"extrap") + + DO ising=1,msing + respsi=sing(ising)%psifac + + ! Evaluate splines on rational surface + CALL spline_eval(sq,respsi,1) + CALL spline_eval(psi_t,respsi,1) + CALL spline_eval(avg_dpsi_spl,respsi,1) + CALL spline_eval(avg_bsq_spl,respsi,1) + CALL spline_eval(v_spl,respsi,1) + CALL spline_eval(vol_spl,respsi,1) + CALL spline_eval(shr_spl,respsi,1) + +c Connor et al. 2015 (PPCF 57 065001) eq. (59): +c Delta_crit = (pi^{3/2}/2) (chi_par/chi_perp)^{1/4} V_s +c * [alpha^2 Lambda^2/(<|grad V|^2>)]^{1/4} +c * (-D_R) +c Identify dgeo = V_s * [alpha^2 Lambda^2 +c /(<|grad V|^2>)]^{1/4}. +c +c Connor's alpha has dimension m^3/Wb: the phase 2 pi i n u / q +c must be dimensionless with u = psi_pol'(V) * theta (Hamada), +c which has units Wb/m^3, so alpha = (2 pi n / q) * v1_SI with +c v1_SI = dV/dpsi_pol [m^3/Wb]. +c Lambda = -q'(V) * (psi_pol'(V))^2 (SI: Wb^2/m^9) +c = -psio^2 * (dq/dpsi_N) / (dV/dpsi_N)^3 +c where psio is total poloidal flux at edge (Wb). +c V_s = V(psi_N_res) (m^3) +c STRIDE's sq uses psi_N (psifac); convert via psi_pol = psi_N * psio. +c sq%f(3) = dV / d(psi_N) [m^3] +c sq%f1(4) = dq / d(psi_N) [dimensionless] +c psio = chi1/(2 pi) [Wb] (total poloidal flux) +c v1_SI = sq%f(3)/psio [m^3/Wb] + V_s = vol_spl%f(1) + q = sq%f(4) + + ! alpha = (2 pi n / q) * v1_SI [m^3/Wb] + al = twopi * nn / q * sq%f(3) / psio + + ! Lambda = -psio^2 * (dq/dpsi_N) / (dV/dpsi_N)^3 [Wb^2/m^9] + Lam = -(psio**2.0) * sq%f1(4) / (sq%f(3)**3.0) + + shr(ising) = respsi * shr_spl%f1(1) +c Final 2 sqrt(2 pi q) factor bridges Connor's Hamada-0-to-1 +c convention (eq. 59) to the Fitzpatrick/PEST3 r_s-normalized +c convention used by STRIDE's Delta_prime (PEST3 matrix) and +c SLAYER's rfitzp dc_type. In the LAR limit this reduces +c dgeo -> sqrt(ns/(R r_s)), matching Connor eq. (61) times r_s. + dgeo(ising) = 2.0d0*SQRT(twopi*q) * V_s + $ * ( (((al**2.0)*(Lam**2.0))/ + $ (avg_bsq_spl%f(1)*avg_dpsi_spl%f(1)))**0.25 ) + ENDDO + ! Deallocate the new cumulative volume spline + CALL spline_dealloc(vol_spl) +c----------------------------------------------------------------------- c optionally write netcdf file. c----------------------------------------------------------------------- IF(present(op_netcdf_out))THEN IF(op_netcdf_out) CALL stride_netcdf_out(wp,wv,wt,ep,ev,et, - $ delta_prime_mat,plasma1,vacuum1,total1) + $ delta_prime_mat,plasma1,vacuum1,total1,shr,dgeo) ENDIF c----------------------------------------------------------------------- c deallocate c----------------------------------------------------------------------- + CALL spline_dealloc(psi_t) + CALL spline_dealloc(avg_dpsi_spl) + CALL spline_dealloc(avg_bsq_spl) + CALL spline_dealloc(v_spl) + CALL spline_dealloc(shr_spl) CALL stride_dealloc c----------------------------------------------------------------------- c terminate. diff --git a/stride/ode.F b/stride/ode.F index 4d9a8281..36f0e7d5 100644 --- a/stride/ode.F +++ b/stride/ode.F @@ -113,7 +113,7 @@ SUBROUTINE ode_run INTEGER, DIMENSION(2*mpert) :: ipiv COMPLEX(r8), DIMENSION(2*mpert,2*mpert) :: uFMInv COMPLEX(r8), DIMENSION(2*mpert) :: uwork - COMPLEX(r8), DIMENSION(2*mpert,2*mpert) :: identityMat + COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: identityMat !NOTE #1: INTEGER, DIMENSION(:), POINTER, PRIVATE :: iwork ! REAL(r8), DIMENSION(:), POINTER, PRIVATE :: rwork @@ -153,8 +153,6 @@ SUBROUTINE ode_run $ ABS(nn*sing(iS)%q1) scalc(iS)%singEdgesLR(2) = sing(iS)%psifac + singfac_min/ $ ABS(nn*sing(iS)%q1) - WRITE(*,'(1x,i5,2(es11.3))') iS,scalc(iS)%singEdgesLR(1), - $ scalc(iS)%singEdgesLR(2) !This finds the index of the singular column scalc(iS)%sing_col = NINT(nn*sing(iS)%q)-mlow+1 @@ -187,6 +185,7 @@ SUBROUTINE ode_run ALLOCATE(rwork(lrw), zwork(lzw), iwork(liw), $ atol(2*mpert,2*mpert)) + ALLOCATE(identityMat(2*mpert,2*mpert)) identityMat = 0.0_r8 DO i = 1,2*mpert identityMat(i,i) = 1.0_r8 @@ -435,6 +434,7 @@ SUBROUTINE ode_run CALL ZVODE1(sing_derFM,neq,uFM,startPsi,endPsi, $ itol,rtol,atol,itask,istate,iopt,zwork,lzw, $ rwork,lrw,iwork,liw,ode_nojac,mf,rpar,ipar) + IF(istate < 0) EXIT ENDDO IF ( grid_packing == "naive" ) THEN CALL SYSTEM_CLOCK(COUNT=fTime) @@ -462,6 +462,7 @@ SUBROUTINE ode_run c----------------------------------------------------------------------- c terminate. c----------------------------------------------------------------------- + DEALLOCATE(identityMat) RETURN END SUBROUTINE ode_run c----------------------------------------------------------------------- diff --git a/stride/stride.F b/stride/stride.F index e2b135e9..d0a0ebe5 100644 --- a/stride/stride.F +++ b/stride/stride.F @@ -358,7 +358,8 @@ PROGRAM stride total1=0_r8 ALLOCATE(mx0(mpert,mpert),vx0(mpert)) CALL stride_netcdf_out(mx0,mx0,mx0,vx0,vx0,vx0, - $ delta_prime_mat,plasma1,vacuum1,total1) + $ delta_prime_mat,plasma1,vacuum1,total1, + $ SPREAD(0._r8,1,msing),SPREAD(0._r8,1,msing)) CALL stride_dealloc ENDIF IF(mat_flag .OR. ode_flag)DEALLOCATE(amat,bmat,cmat,ipiva,jmat) diff --git a/stride/stride_netcdf.f b/stride/stride_netcdf.f index f1f29ac0..0d7d8781 100644 --- a/stride/stride_netcdf.f +++ b/stride/stride_netcdf.f @@ -48,23 +48,30 @@ END SUBROUTINE check c----------------------------------------------------------------------- c declarations. c----------------------------------------------------------------------- - SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) + SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp, + $ pl1,va1,to1,shr,dgeo) REAL(r8), DIMENSION(mpert), INTENT(IN) :: epi,evi,eti COMPLEX(r8), DIMENSION(mpert,mpert), INTENT(IN) :: wp,wv,wt COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: dp REAL(r8), INTENT(IN) :: pl1,va1,to1 + REAL(r8), DIMENSION(msing), INTENT(IN) :: shr,dgeo INTEGER :: i, ncid, $ i_dim, m_dim, mo_dim, p_dim, i_id, m_id, mo_id, p_id, $ f_id, q_id, dv_id, mu_id, di_id, dr_id, ca_id, $ wp_id, wpv_id, wv_id, wvv_id, wt_id, wtv_id, $ r_dim, rp_dim, l_dim, lp_dim, r_id, rp_id, l_id, lp_id, - $ pr_id, qr_id, dp_id, ap_id, bp_id, gp_id, dpp_id + $ pr_id, qr_id, dp_id, ap_id, bp_id, gp_id, dpp_id, + $ shear_id,resm_id,drr_id,dgeo_id COMPLEX(r8), DIMENSION(mpert) :: ep,ev,et CHARACTER(2) :: sn CHARACTER(64) :: ncfile + REAL(r8) :: respsi + REAL(r8), DIMENSION(msing) :: dr_rationals + INTEGER, DIMENSION(msing) :: resm + INTEGER :: ising,jsing COMPLEX(r8), DIMENSION(msing,msing) :: ap,bp,gammap,deltap @@ -85,6 +92,13 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) ep = CMPLX(epi, 0.0) ev = CMPLX(evi, 0.0) et = CMPLX(eti, 0.0) + + ! evaluate resistive interchange parameter on rational surfaces + DO i=1,msing + respsi = sing(i)%psifac + CALL spline_eval(locstab,respsi,0) + dr_rationals(i)=locstab%f(2)/respsi + END DO c----------------------------------------------------------------------- c open files c----------------------------------------------------------------------- @@ -109,6 +123,7 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) CALL check( nf90_put_att(ncid,nf90_global,'mhigh', mhigh)) CALL check( nf90_put_att(ncid,nf90_global,'mpert', mpert)) CALL check( nf90_put_att(ncid,nf90_global,'mband', mband)) + CALL check( nf90_put_att(ncid,nf90_global,'msing', msing)) CALL check( nf90_put_att(ncid,nf90_global,'psilow', psilow)) CALL check( nf90_put_att(ncid,nf90_global,'amean', amean)) CALL check( nf90_put_att(ncid,nf90_global,'rmean', rmean)) @@ -158,6 +173,9 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) CALL check( nf90_def_dim(ncid, "psi_n", sq%mx+1, p_dim) ) CALL check( nf90_def_var(ncid, "psi_n", nf90_double, p_dim, p_id)) IF(msing>0)THEN + DO i=1,msing + resm(i)=NINT(sing(i)%q*nn) + ENDDO CALL check( nf90_def_dim(ncid,"lr_index",2*msing,l_dim) ) CALL check( nf90_def_var(ncid,"lr_index",nf90_int,l_dim,l_id)) CALL check( nf90_put_att(ncid,l_dim,"long_name", @@ -182,7 +200,15 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) CALL check( nf90_def_var(ncid,"q_rational",nf90_double, $ r_dim,qr_id) ) CALL check( nf90_put_att(ncid,qr_id,"long_name", - $ "Safety Factor at Rational Surfaces") ) + $ "Safety Factor at Rational Surfaces") ) + CALL check( nf90_def_var(ncid, "shear", nf90_double, r_dim, + $ shear_id) ) + CALL check( nf90_def_var(ncid, "Delta_geo", nf90_double, r_dim, + $ dgeo_id) ) + CALL check( nf90_def_var(ncid, "resm", nf90_int, r_dim, + $ resm_id) ) + CALL check( nf90_def_var(ncid, "dr_rational", nf90_double, + $ r_dim, drr_id) ) ENDIF ! define variables IF(debug_flag) PRINT *," - Defining variables in netcdf" @@ -246,6 +272,10 @@ SUBROUTINE stride_netcdf_out(wp,wv,wt,epi,evi,eti,dp,pl1,va1,to1) $ i=1,msing)/)) ) CALL check( nf90_put_var(ncid,qr_id, (/(sing(i)%q, $ i=1,msing)/)) ) + CALL check( nf90_put_var(ncid,shear_id, shr) ) + CALL check( nf90_put_var(ncid,resm_id, resm) ) + CALL check( nf90_put_var(ncid,dgeo_id, dgeo) ) + CALL check( nf90_put_var(ncid,drr_id, dr_rationals)) ENDIF IF(debug_flag) PRINT *," - Putting profile variables in netcdf"