module geometry implicit none public private :: bgrad, root, fluxavg real, allocatable, dimension(:) :: grho, theta, bmag, gradpar, & cvdrift, cvdrift0, gbdrift, gbdrift0, gds2, gds21, gds22, & jacob, job, jobp, iob2p real :: rhoc, rmaj, r_geo, shift, dbetadrho real :: qinp, shat, akappa, akappri, tri, tripri, dpressdrho real :: delrho, rmin, rmax real :: s_hat_input, p_prime_input, invLp_input, beta_prime_input, & alpha_input, dp_mult real :: cd_psi0, cd_delpsi, cd_c0, Z_eff integer :: nperiod integer :: itor, iflux, irho real :: dp_new, di_new real, private :: psi_0, psi_a real :: B_T0, avgrmid, dvdrhon, surfarea, grho1n, grho2n, drhodpsin real, dimension(3) :: rpval real :: rpmin, rpmax, ak0 integer :: isym, ismooth, k1, k2, big integer :: eqinit = 1 logical :: gen_eq, vmom_eq, efit_eq, ppl_eq, local_eq, in_nt, writelots, equal_arc integer :: bishop public :: beta_a_fun, eikcoefs, geofax, iofrho, pbarofrho, beta_of_rho, & qfun, rpofrho, rcenter, init_theta, nth_get, f_trap, hirsh_bs !procedures integer, private :: ntgrid, ntheth, nth, nrgrid, ntheta character*80 :: eqfile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Input: !!! rhoc :: minor radius of surface of interest (see irho below) !!! itor :: default one. Choose zero to get "s-alpha" equilibria (see iflux=0) !!! iflux :: chooses mode of operation. !!! case (0) :: no equilibrium. Miller parameterization of local equilibrium !!! iflux = 0 forces local_eq = .true. and irho = 2 !!! case (1) :: use numerical equilibrium. 4 choices: !!! ppl_eq = .true. :: use PPL style NetCDF equilibrium (Menard) !!! gen_eq = .true. :: use GA style NetCDF equilibrium (TOQ) !!! -------------------> [note: if either of the above == .true., !!! set eqfile = input file name] !!! vmom_eq = .true. :: use VMOMS equilibrium !!! efit_eq = .true. :: use EFIT equilibrium !!! case (2) :: running inside nt without actual equilibria !!! iflux = 2 forces local_eq = .true. and irho = 2 !!! case (10) :: used by nt to initialize some variables. Obsolescent. !!! irho :: choose flux surface label !!! case (1) :: sqrt(toroidal flux)/sqrt(toroidal flux of LCFS) !!! case (2) :: diameter/(diameter of LCFS). recommended !!! case (3) :: poloidal flux/(poloidal flux of LCFS) !!! !!! equal_arc :: .true. makes grad_par coefficient a constant. !!! !!! !! New feature: Given an actual equilibrium, allow s_hat and d pressure/d rho !!! !! to vary self-consistently (a la Greene and Chance, Bishop). We use the !!! !! Bishop formalism [C. M. Bishop, et al., NF, Vol. 24, 1579 (1984).] !!! !!! bishop :: uses Bishop relations to generate coefficients. !!! case(0) :: do not use Bishop relations !!! case(1) :: use Bishop relations with actual equilibrium shat, p' !!! case(3) :: use Bishop relations with new shat, L_p from below !!! s_hat_input :: new magnetic shear, rho/q*dq/drho !!! invLp_input :: new -1/pressure * d pressure/d rho !!! case(4) :: use Bishop relations with new shat, beta' from below !!! s_hat_input :: new magnetic shear, rho/q*dq/drho !!! beta_prime_input :: new d beta/d rho !!! case(5) :: use Bishop relations with new shat, alpha from below !!! s_hat_input :: new magnetic shear, rho/q*dq/drho !!! alpha_input :: new alpha = -q**2 * rmaj * d pressure/d rho !!! case(6) :: use Bishop relations with equilibrium shat, beta' from below !!! beta_prime_input :: new d beta/d rho !!! case(7) :: use Bishop relations with equilibrium shat, beta' from below !!! dp_mult :: d pressure/d rho = equilibrium gradient * dp_mult !!! (In each case, the definition of rho is determined by irho above) !!! !!! nperiod :: number of cells of width 2 pi in theta. !!! !!! rmaj :: major radius of magnetic axis (in units of aminor) if iflux = 1 !!! or major radius of center of flux surface of interest (units of aminor) otherwise !!! r_geo :: major radius of last closed flux surface !!! shift :: d rmaj/drho. (Typically < 0.) !!! dIdrho :: d I/drho where I = R B_T !!! qinp :: safety factor q at surface of interest !!! shat :: d q/drho at surface of interest !!! akappa :: elongation (k) of surface of interest !!! akappri :: dk/drho !!! tri :: triangularity of surface of interest !!! tripri :: d tri/drho !!! dpressdrho :: dpressure /drho = 0.5 dbeta /drho at surface of interest. !!! !!! delrho :: numerical parameter. Should be "small enough". Typically 0.001 ok. !!! rmin :: minimum minor radius at which items are evaluated !!! rmax :: should be equal to unity !!! !!! ismooth :: 1 or more -> smooth resulting coefficients. Rarely used. !!! ak0, k1, k2 :: smoothing parameters. Rarely used. !!! !!! isym :: 1 -> assume up-down symmetric equilibrium. !!! in_nt :: .true. if running inside of nt. !!! writelots :: .true. for more output !!! !!! Output: !!! theta :: the theta grid of the results. The final grid has !!! pi - 2*pi*nperiod <= theta <= 2*pi*nperiod - pi !!! bmag :: |B| !!! gradpar :: coefficient of b_hat.grad operator !!! grho :: grad(rho) !!! cvdrift :: coefficient of v_par**2/Omega b_hat x (b_hat.grad b_hat).grad !!! cvdrift0 :: part of total curvature drift operator proportional to theta_0 !!! gbdrift :: coefficient of mu/Omega b_hat x (grad B).grad operator !!! gbdrift0 :: part of total grad B drift operator proportional to theta_0 !!! gds2 :: part of grad_perp**2 operator independent of theta_0 !!! gds21 :: part of grad_perp**2 operator proportional to theta_0 !!! gds22 :: part of grad_perp**2 operator proportional to theta_0**2 !!! jacob :: Jacobian !!! !!! Use: !!! Include this module, set input variables appropriately and call eikcoefs. !!! !!! NOTE: If vmom_eq, gen_eq, and ppl_eq are false, you should call init_theta !!! before eikcoefs to set the number of theta gridpoints per 2 pi. !!! !!! I do not recommend calling other routines, but a few are available for !!! diagnostic purposes. !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! contains subroutine eikcoefs use veq, only: vmomin, veq_init use geq, only: eqin, geq_init use peq, only: peqin => eqin, peq_init use eeq, only: efitin, eeq_init => efit_init use leq, only: leqin, dpdrhofun implicit none ! Local variables: real, allocatable, dimension(:) :: smodel, dsdrp, grho2, seik, seik1, seik2 real, allocatable, dimension(:) :: dsdthet, dsdthet1, dsdthet2, grho1 real, allocatable, dimension(:) :: gdsdum1, gdsdum2, gdsdum3 real, allocatable, dimension(:) :: th_bish, Rpol, ltheta real, allocatable, dimension(:) :: rgrid, rgrid1, rgrid2, Bpolmag, Bmod, dSdl real, allocatable, dimension(:) :: rmajor, ans, ds, arcl real, allocatable, dimension(:,:) :: thgrad, rpgrad, crpgrad, grads real, allocatable, dimension(:,:) :: bvector, gradrptot, gradstot real, allocatable, dimension(:) :: a_bish, b_bish, c_bish, d_bish real, allocatable, dimension(:) :: da_bish, db_bish, dc_bish, dd_bish, trip real :: rtmp, rho, rho1, rho2, rp, rp1, rp2, drhodrp, dqdrho real :: dpsidrp, dpsidrho, dqdrp, dsdrptot, drhodpsi, dqdpsi real :: qval, qval1, qval2, delrp, akpri, shift0, dum real :: pi, bi, Ltot, rinv real :: a_b, b_b, c_b, d_b real :: s_hat_new, s_hat, dp, di, pressure, tmp, ddp, ddi character*1 :: char integer :: i, j, k, itot, ndum, nthg, n, nrg ! compute the initial constants pi=2.*acos(0.) if(iflux == 0 .or. iflux == 2) then if(.not. local_eq) write(*,*) 'Forcing local_eq = true' local_eq = .true. if(irho /= 2) write(*,*) 'Forcing irho = 2' irho = 2 if(bishop == 0) then write(*,*) 'Forcing bishop = 1' bishop = 1 endif endif call check(vmom_eq, gen_eq, efit_eq, ppl_eq, local_eq) if(.not. vmom_eq .and. .not. gen_eq .and. .not. ppl_eq & .and. .not. allocated(theta)) then write(*,*) 'You should call init_theta to specify the ' write(*,*) 'number of theta grid points per 2 pi ' write(*,*) 'before calling eikcoefs.' write(*,*) write(*,*) 'Proceeding with ntheta = 32.' ntheta = 32 call init_theta(ntheta) endif if(iflux /= 1 .and. iflux /= 10) then if(iflux == 2) call radial call leqin(rmaj, R_geo, akappa, akappri, tri, tripri, rhoc, delrho, shift, & qinp, s_hat_input, dpressdrho, ntgrid) if(.not.allocated(gds22)) call alloc_module_arrays(ntgrid) call alloc_local_arrays(ntgrid) endif nrgrid = 1 select case (iflux) case (0) avgrmid=1. case (1, 10) if(vmom_eq) then call vmomin( psi_0, psi_a, rmaj, B_T0, avgrmid, eqinit, in_nt, nthg, nrg) nrgrid = nrg call tdef(nthg) else if(gen_eq) then call eqin(eqfile, psi_0, psi_a, rmaj, B_T0, avgrmid, eqinit, in_nt, nthg, nrg) nrgrid = nrg call tdef(nthg) else if(ppl_eq) then call peqin(eqfile, psi_0, psi_a, rmaj, B_T0, avgrmid, eqinit, in_nt, nthg, nrg) nrgrid = nrg call tdef(nthg) else if(efit_eq) then if(big <= 0) big = 8 call efitin(eqfile, psi_0, psi_a, rmaj, B_T0, avgrmid, eqinit, big) endif if(iflux == 10) return if(.not.allocated(gds22)) call alloc_module_arrays(ntgrid) call alloc_local_arrays(ntgrid) case (2) continue end select if(writelots) then write(11,*) 'All lengths normalized to the avg midplane minor radius' if(irho /= 2) then write(11,*) 'except rho, which always goes from 0 to 1.' write(11,*) 'The rho normalization is determined by irho.' endif write(11,*) 'avgrmid = ',avgrmid,' meters' if(iflux == 1) write(11,*) 'R_mag = ',rmaj if(iflux == 1) write(11,*) 'B_T0 = ',B_T0 endif rho=rhoc rho1 = rhoc + delrho rho2 = rhoc - delrho if(rho1 > 1. .or. rho2 < 0.) then write(*,*) 'Reduce delrho' stop endif if(iflux == 1) then rpmin = psi_0 rpmax = psi_a rp = rpofrho(rho) rp1 = rpofrho(rho1) rp2 = rpofrho(rho2) else rpmin = 0. rpmax = 1.0 rp=rho rp1=rho1 rp2=rho2 endif ! if(iflux == 0 .and. R_geo /= rmaj) then ! write(6,*) ! write(6,*) 'You have set R_geo not equal to Rmaj.' ! write(6,*) 'This is unusual! Make sure you mean it.' ! write(6,*) ! endif if(iflux == 1) R_geo=rcenter(rpmax) if(efit_eq) then call rtg(rgrid, rp) call rtg(rgrid1, rp1) call rtg(rgrid2, rp2) else rgrid = rp rgrid1 = rp1 rgrid2 = rp2 endif if(writelots) then write(11,*) 'rp of rho = ',rp write(11,*) 'Rcenter(rho) = ',rcenter(rp) write(11,*) 'R_geo = ',R_geo,' or ',btori(rgrid(0), 0.) endif if(R_geo > Rmaj + 1.e-5) then write(6,*) write(6,*) 'Warning: the center of the LCFS is outboard of the ' write(6,*) 'center of the reference surface. This is probably wrong.' write(6,*) endif if(writelots) then do i= -nth,nth write(99,*) 'i=',i,' thet= ',theta(i),' rgrid= ',rgrid(i) ! write(99,*) theta(i),Rpos(rgrid(i),theta(i)),Zpos(rgrid(i),theta(i)) enddo endif call drho_drp(rp, drhodrp) if(nperiod > 1) call periodic_copy(theta, 2.*pi) if(writelots) write(11,*) 'drhodrp=',drhodrp ! ! should test whether this is a new equilibrium ! if(gen_eq) call geq_init if(ppl_eq) call peq_init if(vmom_eq) call veq_init if(efit_eq) call eeq_init call rmajortgrid(rgrid, theta, rmajor) ! compute gradient of rp char='P' if(bishop == 0) then call grad(rgrid, theta, rpgrad, char, dum, nth, ntgrid) else call bgrad(rgrid, theta, rpgrad, char, dum, nth, ntgrid) endif if(bishop /= 0) then call loftheta(rgrid, theta, ltheta) call grad(rgrid, theta, crpgrad, char, dum, nth, ntgrid) call th_bishop(crpgrad, th_bish, nth) call B_pol(rgrid, theta, crpgrad, Bpolmag, nth) call B_mod(rgrid, theta, Bpolmag, Bmod, nth) call R_pol(theta, th_bish, ltheta, Rpol, nth) if(nperiod > 1) then call periodic_copy(rmajor, 0.) call periodic_copy(Bpolmag, 0.) call periodic_copy(Bmod, 0.) call periodic_copy(th_bish, 0.) call periodic_copy(Rpol, 0.) Ltot = ltheta(nth)-ltheta(-nth) call periodic_copy(ltheta, Ltot) endif endif ! compute coordinate gradients call thetagrad(rgrid,thgrad) ! compute eikonal S if(iflux == 0 .or. iflux == 2) qval=qfun(pbarofrho(rho)) call eikonal(rgrid, rpgrad, thgrad, qval, seik, dsdthet, dpsidrp) if(writelots) write(11,*) 'q= ',qval if(iflux == 0 .or. iflux == 2) then bi = btori(rgrid(0), theta(0)) drhodpsin = drhodrp/dpsidrp do i=-nth,nth rinv = invRfun(rgrid(i), theta(i)) Bpolmag(i) = sqrt(rpgrad(i,1)**2+rpgrad(i,2)**2)*rinv*abs(dpsidrp) Bmod(i) = sqrt( (bi*rinv)**2 + bpolmag(i)**2) enddo if(nperiod > 1) then call periodic_copy(Bpolmag, 0.) call periodic_copy(Bmod, 0.) endif endif ! do i=-nth,nth ! write(43,*) theta(i), Bpolmag(i), rpos(rgrid(i), theta(i)), zpos(rgrid(i), theta(i)), seik(i) ! enddo ! ! ****************** ! Debug: ! Set itor=1 and the first two columns should be the ! same (better agreement results from higher ntheta) ! Set itor=0 and the first and third columns should ! be equal. All of this, of course, with iflux=0. ! ! This amounts to checking the expression for ! nu for axisymmetric, low beta equilibrium: ! nu = - q r/R sin(theta). ! ! call rmajortgrid(rgrid, theta, rmajor) ! do i=-nth,nth ! smodel(i)=-qval*(theta(i)-rgrid(i)*sin(theta(i))/rmajor(i)) ! write(*,*) seik(i),smodel(i),-qval*theta(i) ! enddo ! ****************** if(writelots) write(11,*) 'dpsidrp= ',dpsidrp ! ! compute derivatives of eikonal if(bishop /= 0) then drhodpsi = drhodrp / dpsidrp n=ntgrid if(.not. allocated(a_bish)) then allocate(a_bish(-n:n), b_bish(-n:n), c_bish(-n:n), d_bish(-n:n), & da_bish(-n:n), db_bish(-n:n), dc_bish(-n:n), dd_bish(-n:n), trip(-n:n)) endif bi = btori(rgrid(0),theta(0)) da_bish = 1/rmajor**2 + (bi/(Bpolmag*rmajor**2))**2 db_bish = bi/(Bpolmag*rmajor)**2 dc_bish = bi*(sin(th_bish) + rmajor/Rpol)/(-Bpolmag*rmajor**4) call tripprod2dtgrid(rpgrad, thgrad, rgrid, trip) if (nperiod > 1) call periodic_copy(trip, 0.) if(iflux == 0 .or. iflux == 2) trip = trip * dpsidrp da_bish = da_bish/trip db_bish = db_bish/trip dc_bish = dc_bish/trip call integrate(da_bish, theta, a_bish, ntgrid) call integrate(db_bish, theta, b_bish, ntgrid) call integrate(dc_bish, theta, c_bish, ntgrid) a_b = a_bish(nth)-a_bish(-nth) b_b = b_bish(nth)-b_bish(-nth) c_b = c_bish(nth)-c_bish(-nth) a_bish = - a_bish * Bpolmag * rmajor b_bish = - b_bish * Bpolmag * rmajor c_bish = - c_bish * Bpolmag * rmajor ! find shat, pressure' and report them if(iflux == 0 .or. iflux == 2) then dp = dpdrhofun(rhoc)*drhodpsi else dp = dpfun(rgrid(0),theta(0)) endif di = dbtori(rgrid(0),theta(0)) tmp = rho/qval/(2.*pi)/drhodpsi s_hat = tmp*(a_b*di+ b_b*dp + c_b*2) pressure = pfun(rgrid(0),0.) write(11,*) write(11,*) 'Quantity, equilibrium value, value used' select case (bishop) case (2) s_hat_new = s_hat_input dp_new = p_prime_input di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b write(11,*) 'p_prime = ',dp,', ',dp_new if(writelots) write(*,*) 'p_prime = ',dp,', ',dp_new case (3) s_hat_new = s_hat_input dp_new = -invLp_input*pressure*drhodpsi di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b write(11,*) '1/L_p = ',-dp/drhodpsi/pressure,', ',-dp_new/drhodpsi/pressure if(writelots) write(*,*) '1/L_p = ',-dp/drhodpsi/pressure,', ',-dp_new/drhodpsi/pressure case (4) s_hat_new = s_hat_input dp_new = 0.5*beta_prime_input*drhodpsi di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b write(11,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi if(writelots) write(*,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi case (5) s_hat_new = s_hat_input dp_new = -alpha_input/qval**2/rmaj*drhodpsi di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b write(11,*) 'alpha = ', -qval**2*rmaj*dp/drhodpsi,', ', & -qval**2*rmaj*dp_new/drhodpsi if(writelots) write(*,*) 'alpha = ', -qval**2*rmaj*dp/drhodpsi,', ', & -qval**2*rmaj*dp_new/drhodpsi case (6) s_hat_new = s_hat dp_new = 0.5*beta_prime_input*drhodpsi di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b write(11,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi if(writelots) write(*,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi case (7) s_hat_new = s_hat dp_new = dp * dp_mult di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b write(11,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi if(writelots) write(*,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi case (8) s_hat_new = s_hat_input dp_new = dp * dp_mult di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b write(11,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi if(writelots) write(*,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi case (9) s_hat_new = s_hat_input dp_new = 0.5*beta_prime_input*drhodpsi di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b write(11,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi if(writelots) write(*,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi case default if(iflux /= 1) then s_hat_new = s_hat_input dp_new = dp di_new = (s_hat_new/tmp -2*c_b -dp_new*b_b) / a_b else s_hat_new = s_hat dp_new = dp di_new = di endif write(11,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi if(writelots) write(*,*) 'd beta/d rho = ',2.*dp/drhodpsi,', ',2.*dp_new/drhodpsi end select write(11,*) 's_hat ',s_hat,', ',s_hat_new if(writelots) write(*,*) 's_hat ',s_hat,', ',s_hat_new write(11,*) shat = s_hat_new dqdrp = s_hat_new*qval/rho*drhodrp dbetadrho = 2*dp_new/drhodpsi call gradl(ltheta, seik, dSdl, 2*pi*qval, nth) if(nperiod > 1) call periodic_copy(dSdl, 0.) grads(:,1) = a_bish*di_new + b_bish*dp_new + c_bish*2 grads(:,2) = dSdl ! High-n kink (peeling) mode terms: ! ! One could do a kind of s-alpha model here, with p'' and I'', but ! I will look at that later. ! ! We have dI/dPsi, dp/dPsi. For kink-type terms, we will need a little more: ! (remember, bpolmag -> -bpolmag b/c of different definition in Bishop) ! Coefficient of I'' term: da_bish = 1. ! Coefficient of p'' term: db_bish = bi/Bmod**2 ! Remaining part: dc_bish = di_new*dp_new/Bmod**2 & + 2.*bi*dp_new/bmod**4 & * (-Bpolmag/(rmajor*Rpol) + dp_new & -bi**2*sin(th_bish)/(-Bpolmag*rmajor**4)) ! find (J_par/B)', pressure'' and report them if(iflux == 0 .or. iflux == 2) then ! dp = dpdrhofun(rhoc)*drhodpsi write(*,*) 'error: ddp not calculated' else ddp = (dpfun(rgrid1(0),0.)-dpfun(rgrid2(0),0.)) & /(psi(rgrid1(0),0.)-psi(rgrid2(0),0.)) endif ! This returns nonsense for iflux /= 1, but it is ok b/c ! it should be recalculated in that case ddi = (dbtori(rgrid1(0),0.)-dbtori(rgrid2(0),0.)) & /(psi(rgrid1(0),0.)-psi(rgrid2(0),0.)) ! da_bish = da_bish/trip ! db_bish = db_bish/trip ! dc_bish = dc_bish/trip ! call integrate(da_bish, theta, a_bish, ntgrid) ! call integrate(db_bish, theta, b_bish, ntgrid) ! call integrate(dc_bish, theta, c_bish, ntgrid) ! a_b = a_bish(nth)-a_bish(-nth) ! b_b = b_bish(nth)-b_bish(-nth) ! c_b = c_bish(nth)-c_bish(-nth) job = di_new + bi*dp_new/bmod**2 jobp = (da_bish*ddi + db_bish*ddp + dc_bish)/drhodpsi iob2p = dc_bish/drhodpsi/dp_new ! Connor peeling mode stability condition if(writelots .and. .not. in_nt .and. nperiod == 1) then dd_bish = - job*(bi/(Bpolmag*rmajor**2))**2/trip call integrate(dd_bish, theta, d_bish, ntgrid) d_b = d_bish(nth)-d_bish(-nth) dqdpsi = s_hat_new*qval/rho*drhodpsi if(rho > 0.9) then write(*,*) write(*,*) 'This should be evaluated at the edge:' write(*,*) 'Peeling mode stability' write(*,*) 'sqrt(1 - 4 D_M) > RHS' write(*,*) 'RHS = ',1.+2./(2.*pi*dqdpsi)*d_b write(*,*) endif endif else ! rp derivative delrp=delrho/drhodrp rp1=rp-delrp rp2=rp+delrp if (iflux == 0 .or. iflux == 2) then dqdrho=dqdrhofun(rho) qval1=qval-dqdrho*delrho qval2=qval+dqdrho*delrho endif call seikon(rp1, qval1, seik1, dsdthet1, rgrid1, ltheta, bpolmag, dum) call seikon(rp2, qval2, seik2, dsdthet2, rgrid2, ltheta, bpolmag, dum) do i=-nth,nth dsdrp(i)=(seik2(i)-seik1(i))/(2.*delrp) enddo if(iflux == 1) then dqdrho=(qval2-qval1)/(2.*delrho) shat=dqdrho*rho/qval endif if(writelots) then write(11,*) 'dqdrho=',dqdrho,' qval1,2= ',qval1,qval2 write(11,*) 's_hat= ',shat write (11,*) 'i,rgrid,rgrid1,rgrid2:' do i=-nth,nth write (11,*) i,rgrid(i),rgrid1(i),rgrid2(i) enddo endif rpval(1)=rp rpval(2)=rp1 rpval(3)=rp2 if(iflux == 0 .or. iflux == 2) then dbetadrho = 2.*dpdrhofun(rhoc) else dbetadrho = 2.*dpfun(rgrid(0),theta(0))/drhodrp endif ! ! compute gradient of S dqdrp=dqdrho*drhodrp do k=-nperiod+1,nperiod-1 do j=1,2 do i=-nth,nth itot=i+k*ntheta dsdrptot=dsdrp(i)-k*2*pi*dqdrp grads(itot,j)=rpgrad(i,j)*dsdrptot+thgrad(i,j)*dsdthet(i) enddo enddo enddo endif if(writelots) then write(11,*) 'i,theta,grads1,2= ' do i=-ntgrid,ntgrid write(11,*) i,theta(i),grads(i,1),grads(i,2) enddo endif ! compute total grad S including toroidal part call gradstottgrid(rmajor, grads, gradstot) ! compute the magnitude squared of grad S dpsidrho=dpsidrp/drhodrp if(iflux == 1) drhodpsin = drhodrp/dpsidrp call dottgridf(gradstot, gradstot, gdsdum1) do k=-nperiod+1,nperiod-1 do i=-nth,nth itot=i+k*ntheta gradrptot(itot,1)=rpgrad(i,1) gradrptot(itot,2)=rpgrad(i,2) gradrptot(itot,3)=0. enddo enddo call dottgridf(gradstot, gradrptot, gdsdum2) call dottgridf(gradrptot, gradrptot, gdsdum3) gds2 = gdsdum1 *dpsidrho**2 gds21 = gdsdum2*dqdrp *dpsidrho**2 gds22 = gdsdum3*dqdrp*dqdrp*dpsidrho**2 ! compute magnetic field call bvectortgrid(rgrid, theta, nth, rpgrad, dpsidrp, bvector) ! compute curvature and grad b drift terms call drift(rgrid, rp, bvector, gradstot, gradrptot, & dqdrp, dpsidrho, drhodrp, Bmod, Bpolmag, Rpol, th_bish, ltheta) ! compute the magnetic field along theta ! call bmagtgrid(rgrid, bmagtg) if(ismooth >= 1) then call smoothie(2*nth+1, Bmod(-nth:nth), Bmod(-nth:nth), k1, k2, ak0) Bmod(nth)=Bmod(-nth) endif ! if(itor /= 0) then ! call periodic_copy(bmagtg, 0.) ! bmag = bmagtg ! if(nperiod > 1 ) call periodic_copy(Bmod, 0.) bmag = Bmod ! else ! do k= -nperiod+1,nperiod-1 ! do i=-nth,nth ! itot= i+k*ntheta ! bmag(itot)= 1./(1.+rgrid(i)*cos(theta(i))/rmaj) ! enddo ! enddo ! endif if(isym == 1) call sym(bmag, 0, ntgrid) ! compute the gradparallel coefficient bi = btori(rgrid(0),theta(0)) do k=-nperiod+1,nperiod-1 do i=-nth,nth itot=i+k*ntheta gradpar(itot) = -bi/(rmajor(i)**2*bmag(i)*dsdthet(i)) enddo enddo if(isym == 1) call sym(gradpar, 0, ntgrid) if(equal_arc) then call arclength (ntheta, nperiod, gradpar, arcl) theta = arcl endif ! compute |grad rho|: do k=-nperiod+1,nperiod-1 do i=-nth,nth itot=i+k*ntheta grho(itot)=sqrt(rpgrad(i,1)**2+rpgrad(i,2)**2)*abs(drhodrp) enddo enddo if(isym == 1) call sym(grho, 0, ntgrid) ! compute grad rho*Jacobian of the transformation to the new coordinates: ! Axisymmetry assumed to get this (nu indep. of phi): do i=-nth,nth jacob(i)=abs(dpsidrho/(gradpar(i)*bmag(i))) ds(i)=grho(i)*jacob(i) ! write(*,*) theta(i),' Jac_eik= ',jacob(i) enddo ! compute surface area from A=Int(R sqrt(r**2+drdtheta**2) dtheta dphi) ! does not work with equal_arc grid ! surfarea=surfareafun(rgrid) ! write(*,*) 'surface area= ',surfarea,' avgrmid**2' ! compute the surface area from A=Int(J |grad rho| dtheta dalpha) call integrate(ds, theta, ans, nth) surfarea=2.*pi*(ans(nth)-ans(-nth)) ! write(*,*) 'surface area= ',surfarea,' avgrmid**2' ! compute dV/drhon = 2 pi Int(J dtheta) call integrate(jacob, theta, ans, nth) dvdrhon=2.*pi*(ans(nth)-ans(-nth)) ! write(*,*) 'dV/drhon= ',dvdrhon ! write(*,*) '< |grad rho| >= ',surfarea/dvdrhon ! write(*,*) ! compute < |grad rho|**2> do i=-nth,nth grho2(i)=grho(i)**2*jacob(i) enddo call integrate(grho2, theta, ans, nth) grho2n=2.*pi*(ans(nth)-ans(-nth)) if(isym == 1) then call sym(gds2, 0, ntgrid) call sym(gds21, 1, ntgrid) call sym(gds22, 0, ntgrid) endif ! write(*,*) hirsh_bs(bmag(-nth:nth), 2., 0.5)*bi*dp_new write(25,*) rhoc, f_trap(bmag(-nth:nth)) if(eqinit >= 1) eqinit=0 call dealloc_local_arrays contains subroutine alloc_local_arrays(n) integer n allocate(smodel (-n:n), & dsdrp (-n:n), & grho2 (-n:n), & grho1 (-n:n), & seik (-n:n), & seik1 (-n:n), & seik2 (-n:n), & dsdthet (-n:n), & dsdthet1 (-n:n), & dsdthet2 (-n:n), & gdsdum1 (-n:n), & gdsdum2 (-n:n), & gdsdum3 (-n:n), & th_bish (-n:n), & Rpol (-n:n), & ltheta (-n:n), & rgrid (-n:n), & rgrid1 (-n:n), & rgrid2 (-n:n), & Bpolmag (-n:n), & Bmod (-n:n), & dSdl (-n:n), & rmajor (-n:n), & ans (-n:n), & ds (-n:n), & arcl (-n:n) ) allocate(thgrad (-n:n,2), & rpgrad (-n:n,2), & crpgrad (-n:n,2), & grads (-n:n,2)) allocate(bvector(-n:n,3), & gradrptot (-n:n,3), & gradstot (-n:n,3)) end subroutine alloc_local_arrays subroutine dealloc_local_arrays deallocate(smodel, & dsdrp , & grho2 , & grho1 , & seik , & seik1 , & seik2 , & dsdthet , & dsdthet1 , & dsdthet2 , & gdsdum1 , & gdsdum2 , & gdsdum3 , & th_bish , & Rpol , & ltheta , & rgrid , & rgrid1 , & rgrid2 , & Bpolmag , & Bmod , & dSdl , & rmajor , & ans , & ds , & arcl ) deallocate(thgrad, & rpgrad , & crpgrad , & grads ) deallocate(bvector, & gradrptot , & gradstot ) end subroutine dealloc_local_arrays end subroutine eikcoefs function surfareafun(rgrid) real surfareafun real, dimension(:) :: rgrid real, dimension(-ntgrid:ntgrid) :: ds, drdth, ans real pi integer i write(*,*) 'surfareafun not generalized yet? needs to be checked.' pi=2.*acos(0.) drdth(-nth)=(rgrid(nth-1)-rgrid(-nth+1))/(theta(nth-1)-theta(-nth+1)) do i=-nth+1,nth-1 drdth(i)=(rgrid(i+1)-rgrid(i-1))/(theta(i+1)-theta(i-1)) enddo drdth(nth)=(rgrid(-nth+1)-rgrid(nth-1))/(theta(-nth+1)-theta(nth-1)) do i=-nth,nth ds(i)=2.*pi/invRfun(rgrid(i),theta(i))*sqrt(rgrid(i)**2+drdth(i)**2) enddo call integrate(ds, theta, ans, nth) surfareafun=ans(nth)-ans(-nth) end function surfareafun subroutine bvectortgrid(rgrid, theta, nth, rpgrad, dpsidrp, bvector) integer, intent (in) :: nth real, dimension(-ntgrid: ), intent (in) :: rgrid, theta real, dimension(-ntgrid:,:), intent (in) :: rpgrad real, dimension(-ntgrid:,:), intent (out) :: bvector real, intent (in) :: dpsidrp real :: bi real, dimension(-ntgrid:ntgrid) :: Rinv integer :: i bi = btori(rgrid(0), theta(0)) do i = -nth, nth Rinv(i) = invRfun(rgrid(i), theta(i)) enddo do i=-nth,nth bvector(i,1) =-dpsidrp*rpgrad(i,2) bvector(i,2) = dpsidrp*rpgrad(i,1) bvector(i,3) = bi enddo do i=-nth,nth bvector(i,1) = bvector(i,1) * Rinv(i) bvector(i,2) = bvector(i,2) * Rinv(i) bvector(i,3) = bvector(i,3) * Rinv(i) enddo return end subroutine bvectortgrid subroutine gradstottgrid (rmajor, grads, gradstot) real, dimension(-ntgrid:), intent (in) :: rmajor real, dimension(-ntgrid:, :), intent (in) :: grads real, dimension(-ntgrid:, :), intent (out) :: gradstot integer :: i, k, itot do k=-nperiod+1,nperiod-1 do i=-nth,nth itot=i+k*ntheta gradstot(itot,1)=grads(itot,1) gradstot(itot,2)=grads(itot,2) if(itor == 0) then gradstot(itot,3)=0. else gradstot(itot,3)=1./rmajor(i) endif enddo enddo end subroutine gradstottgrid subroutine crosstgrid(a, b, c) real, dimension(-ntgrid:, :), intent (in) :: a, b real, dimension(-ntgrid:, :), intent (out) :: c integer :: i do i = -nth, nth c(i,1)=a(i,2)*b(i,3)-b(i,2)*a(i,3) c(i,2)=a(i,3)*b(i,1)-b(i,3)*a(i,1) c(i,3)=a(i,1)*b(i,2)-b(i,1)*a(i,2) enddo end subroutine crosstgrid subroutine dottgrid(a, b, c) real, dimension (-ntgrid:, :) :: a, b real, dimension (-ntgrid:) :: c integer :: i, k, itot do k = -nperiod+1, nperiod-1 do i = -nth, nth itot = i + k*ntheta c(itot) = a(i,1)*b(itot,1) + a(i,2)*b(itot,2) + a(i,3)*b(itot,3) enddo enddo end subroutine dottgrid subroutine dottgridf(a, b, c) real a(-ntgrid:,:),b(-ntgrid:,:),c(-ntgrid:) integer i c = a(:,1)*b(:,1) + a(:,2)*b(:,2) + a(:,3)*b(:,3) end subroutine dottgridf subroutine bmagtgrid(rgrid, bmagtg) real, dimension (-ntgrid:), intent (in) :: rgrid real, dimension (-ntgrid:), intent (out) :: bmagtg integer i do i=-nth,nth bmagtg(i)=bmagfun(rgrid(i),theta(i)) enddo if(ismooth >= 1) then call smoothie(2*nth+1,bmagtg(-nth:nth),bmagtg(-nth:nth),k1,k2,ak0) bmagtg(nth)=bmagtg(-nth) endif end subroutine bmagtgrid real function bmagfun(r, thet) real, intent (in) :: r, thet real :: bt if(itor == 0) then bmagfun = abs(btori(r, thet)*invRfun(r, thet)) else if(iflux == 1 .and. (gen_eq .or. vmom_eq)) then bmagfun = bmodfun(r, thet) ! needs to be upgraded!!! else bt = btori(r, thet)*invRfun(r, thet) write(*,*) 'error: bmagfun not defined.' stop ! bmagfun = sqrt(bt**2+bpmagfun(r, thet)**2) endif end function bmagfun function Rpos(r, thet) use geq, only: geq_R => Rpos use peq, only: peq_R => Rpos use leq, only: leq_R => Rpos use eeq, only: eeq_R => Rpos use veq, only: veq_R => Rpos real, intent (in) :: r, thet real :: rpos if(gen_eq) Rpos = geq_R (r, thet) if(ppl_eq) Rpos = peq_R (r, thet) if(vmom_eq) Rpos = veq_R (r, thet) if(efit_eq) Rpos = eeq_R (r, thet) if(local_eq)Rpos = leq_R (r, thet) end function Rpos function Zpos(r, thet) use geq, only: geq_Z => Zpos use peq, only: peq_Z => Zpos use leq, only: leq_Z => Zpos use eeq, only: eeq_Z => Zpos use veq, only: veq_Z => Zpos real, intent (in) :: r, thet real :: Zpos if(gen_eq) Zpos = geq_Z (r, thet) if(ppl_eq) Zpos = peq_Z (r, thet) if(vmom_eq) Zpos = veq_Z (r, thet) if(efit_eq) Zpos = eeq_Z (r, thet) if(local_eq)Zpos = leq_Z (r, thet) end function Zpos function invRfun(r, thet) use geq, only: geq_invR => invR use peq, only: peq_invR => invR use leq, only: leq_invR => invR use eeq, only: eeq_invR => invR use veq, only: veq_invR => invR real, intent (in) :: r, thet real :: invRfun integer :: i ! if(itor == 0) then ! invRfun=1./rmaj ! else if(gen_eq) invRfun = geq_invR (r, thet) if(ppl_eq) invRfun = peq_invR (r, thet) if(vmom_eq) invRfun = veq_invR (r, thet) if(efit_eq) invRfun = eeq_invR (r, thet) if(local_eq)invRfun = leq_invR (r, thet) ! endif end function invRfun subroutine drift(rgrid, rp, bvector, gradstot, & gradrptot, dqdrp, dpsidrho, drhodrp, Bmod, & Bpolmag, Rpol, th_bish, ltheta) real, dimension (-ntgrid:), intent (in) :: rgrid, Bmod, Bpolmag, Rpol, th_bish, ltheta real, dimension (-ntgrid:, :), intent (in) :: bvector, gradstot, gradrptot real, intent (in) :: rp, dqdrp, dpsidrho, drhodrp real, dimension (-ntgrid:ntgrid, 3) :: bgradtot, pgradtot, igradtot, & dummy, dummy1, curve real, dimension (-ntgrid:ntgrid, 2) :: pgrad, igrad, bgrad1, bgrad2, bgrad3 real, dimension (-ntgrid:ntgrid) :: gbdrift1, gbdrift2, cvdrift1, cvdrift2, & dummy2, dell real, dimension (2*ntgrid + 1) :: dumdum1, dumdum2 real :: pressure, err, val, val0, bt, dum, r, t, bi character*1 char integer :: i, k, itot, ndum ndum = 2*nth + 1 ! if (itor == 0) then ! gbdrift = (-sin(theta)*gradstot(:,1)-cos(theta)*gradstot(:,2))/rmaj ! gbdrift=2.*dpsidrho*gbdrift ! gbdrift0=(-sin(theta)*gradrptot(:,1)-cos(theta)*gradrptot(:,2))/rmaj ! gbdrift0=2.*dpsidrho*gbdrift0*dqdrp ! cvdrift=gbdrift ! cvdrift0=gbdrift0 ! ! else char='B' if(bishop >= 1) then call bishop_gradB(rgrid, Bmod, Bpolmag, Rpol, th_bish, ltheta, bgrad1) else call grad(rgrid, theta, bgrad1, char, dum, nth, ntgrid) endif if(ismooth >= 2) then do i=-nth,nth dumdum1(nth+i+1)=bgrad1(i,1) dumdum2(nth+i+1)=bgrad1(i,2) enddo call smoothie(ndum,dumdum1(1:ndum),dumdum1(1:ndum),k1,k2,ak0) call smoothie(ndum,dumdum2(1:ndum),dumdum2(1:ndum),k1,k2,ak0) do i=-nth,nth bgrad1(i,1)=dumdum1(i+nth+1) bgrad1(i,2)=dumdum2(i+nth+1) enddo endif do i=-nth,nth bgradtot(i,1)=bgrad1(i,1) bgradtot(i,2)=bgrad1(i,2) bgradtot(i,3)=0. enddo ! create curvature char='R' select case (bishop) case (0) call grad(rgrid, theta, pgrad, char, rp, nth, ntgrid) case (1) pgradtot = gradrptot*dpsidrho*drhodrp*dpfun(rgrid(0), 0.) case default pgradtot = gradrptot*dpsidrho*drhodrp*dp_new end select if(bishop < 1) then do i=-nth,nth pgradtot(i,1) = pgrad(i,1) pgradtot(i,2) = pgrad(i,2) pgradtot(i,3) = 0. enddo endif do i=-nth,nth curve(i,1)=bgradtot(i,1)/bmod(i)+pgradtot(i,1)/bmod(i)**2 curve(i,2)=bgradtot(i,2)/bmod(i)+pgradtot(i,2)/bmod(i)**2 curve(i,3)=0. enddo ! write(*,*) 'curve 1, 2 ',curve(0,1), curve(0,2) call crosstgrid(bvector,bgradtot,dummy1) call dottgrid(dummy1,gradstot,gbdrift1) call dottgrid(dummy1,gradrptot,gbdrift2) call crosstgrid(bvector,curve,dummy) call dottgrid(dummy,gradstot,cvdrift1) call dottgrid(dummy,gradrptot,cvdrift2) do k=-nperiod+1,nperiod-1 do i=-nth,nth itot=i+k*ntheta gbdrift(itot) =2.*dpsidrho*gbdrift1(itot)/bmod(i)**3 gbdrift0(itot)=2.*dpsidrho*gbdrift2(itot)*dqdrp/bmod(i)**3 cvdrift(itot) =2.*dpsidrho*cvdrift1(itot)/bmod(i)**2 cvdrift0(itot)=2.*dpsidrho*cvdrift2(itot)*dqdrp/bmod(i)**2 enddo enddo ! endif if(isym == 1) then call sym(gbdrift, 0, ntgrid) call sym(cvdrift, 0, ntgrid) call sym(gbdrift0, 1, ntgrid) call sym(cvdrift0, 1, ntgrid) endif end subroutine drift subroutine sym(a, isign, ntgrid) integer i, isign, ntgrid real, dimension(-ntgrid:), intent (in out) :: a if(isign == 0) then do i=1,ntgrid a(i)=0.5*(a(i)+a(-i)) a(-i)=a(i) enddo else do i=1,ntgrid a(i)=0.5*(a(i)-a(-i)) a(i)=-a(-i) enddo a(0)=0. endif end subroutine sym subroutine seikon(rp, qval, seik, dsdthet, rgrid, ltheta, bpolmag, dpsidrp) real :: rp, qval real, dimension(-ntgrid:) :: seik, dsdthet, rgrid, ltheta, bpolmag real, dimension(-ntgrid:ntgrid, 2) :: thgrad, rpgrad real :: dpsidrp, dum character*1 char rgrid = rp char='P' if (bishop == 0) then call grad(rgrid, theta, rpgrad, char, dum, nth, ntgrid) else call bgrad(rgrid, theta, rpgrad, char, dum, nth, ntgrid) end if call thetagrad(rgrid, thgrad) call eikonal(rgrid, rpgrad, thgrad, qval, seik, dsdthet, dpsidrp) end subroutine seikon subroutine eikonal(rgrid, rpgrad, thgrad, qval, seik, dsdthet, dpsidrp) real, dimension(-ntgrid:ntgrid) :: trip, seik, dsdthet, rgrid real, dimension(-ntgrid:, :) :: rpgrad, thgrad real :: qval, dpsidrp, pi, bi integer :: i pi=2.*acos(0.) call tripprod2dtgrid(rpgrad, thgrad, rgrid, trip) bi = btori(rgrid(0),theta(0)) do i=-nth,nth dsdthet(i)=-bi/trip(i)*invRfun(rgrid(i),theta(i))**2 enddo if(nperiod > 1) call periodic_copy(dsdthet, 0.) call integrate(dsdthet, theta, seik, ntgrid) if(iflux == 0 .or. iflux == 2) then dpsidrp=-(seik(nth)-seik(-nth))/(2.*pi*qval) do i=-nth,nth seik(i)=seik(i)/dpsidrp dsdthet(i)=dsdthet(i)/dpsidrp enddo else dpsidrp=1. qval=-(seik(nth)-seik(-nth))/(2.*pi) endif end subroutine eikonal subroutine thetagrad(rgrid, thgrad) real, dimension(-ntgrid:), intent (in) :: rgrid real, dimension(-ntgrid:, :), intent (out) :: thgrad real :: dum character*1 char integer :: i char='T' if(bishop == 0) then if(efit_eq) then write(*,*) 'error in thetagrad' stop endif call grad(rgrid, theta, thgrad, char, dum, nth, ntgrid) else call bgrad(rgrid, theta, thgrad, char, dum, nth, ntgrid) end if if(nperiod > 1) then call periodic_copy(thgrad(-ntgrid:ntgrid,1),0.) call periodic_copy(thgrad(-ntgrid:ntgrid,2),0.) endif end subroutine thetagrad function diameter(rp) use geq, only: geq_diameter => diameter, geq_init_diameter => initialize_diameter use peq, only: peq_diameter => diameter, peq_init_diameter => initialize_diameter use veq, only: veq_diameter => diameter, veq_init_diameter => initialize_diameter use eeq, only: eeq_diameter => diameter integer :: i, initd = 1 real :: rp, pi, diameter pi=2.*acos(0.) if(rp <= rpmin .and. .not. efit_eq) then diameter = 0. return endif if(eqinit == 1) initd = 1 if(iflux /= 1) then diameter=2.*rp else if (gen_eq) i=geq_init_diameter(initd) if (ppl_eq) i=peq_init_diameter(initd) if (vmom_eq) i=veq_init_diameter(initd) initd=0 if(gen_eq) diameter = geq_diameter(rp) if(ppl_eq) diameter = peq_diameter(rp) if(vmom_eq) diameter = veq_diameter(rp) if(efit_eq) diameter = eeq_diameter(rp) endif end function diameter function psibarfun (rp) real, intent (in) :: rp real :: psibarfun psibarfun=min(1.,max(0.,(rp-psi_0)/(psi_a-psi_0))) end function psibarfun function rpfun(r, thet) real, intent (in) :: r, thet real :: rpfun rpfun=psi(r, thet) end function rpfun function rpofrho(rho) real :: rpofrho real :: a, b, xerrbi, xerrsec, fval, soln, rho integer :: nsolv,ier a=rpmin b=rpmax xerrbi=1.e-8 xerrsec=1.e-8 nsolv=1 if(irho == 1) then fval=rho*rho*phi(b) call root(phi, fval, a, b, xerrbi, xerrsec, nsolv, ier, soln) elseif(irho == 2) then fval=rho*diameter(b) call root(diameter, fval, a, b, xerrbi, xerrsec, nsolv, ier, soln) elseif(irho == 3) then fval=rho*psibarfun(b) call root(psibarfun, fval, a, b, xerrbi, xerrsec, nsolv, ier, soln) endif rpofrho=soln if(ier > 0) write(11,*) 'error in rpofrho,rho=',rho end function rpofrho subroutine tripprod2dtgrid(x, y, rgrid, val) real, dimension(-ntgrid:), intent (in) :: rgrid real, dimension(-ntgrid:), intent (out) :: val real, dimension(-ntgrid:, :), intent(in) :: x, y integer :: i ! factor of 1/R from grad zeta do i=-nth,nth val(i)=(x(i,1)*y(i,2)-x(i,2)*y(i,1))*invRfun(rgrid(i),theta(i)) enddo end subroutine tripprod2dtgrid function btori(r, thet) use veq, only: veq_btori => btori , vm_init_btori => initialize_btori use geq, only: geq_btori => btori , geq_init_btori => initialize_btori use peq, only: peq_btori => btori , ppl_init_btori => initialize_btori use eeq, only: eeq_btori => btori, efit_init_btori => initialize_btori use leq, only: leq_btori => btori real :: btori real, intent (in) :: r, thet real :: pbar, f real, save :: r_last, theta_last, I_last integer :: i, initb = 1 ! ! In the present code, most calls to this routine have the same r, thet, so: ! if(r == r_last .and. thet == theta_last .and. eqinit /= 1) then btori = I_last return endif if(eqinit == 1) initb = 1 if(iflux == 1) then pbar=min(1.,max(0.,(rpfun(r,thet)-psi_0)/(psi_a-psi_0))) if(vmom_eq) then i=vm_init_btori(initb) f=veq_btori(pbar) elseif (gen_eq) then i=geq_init_btori(initb) f=geq_btori(pbar) elseif (ppl_eq) then i=ppl_init_btori(initb) f=peq_btori(pbar) elseif (efit_eq) then i = efit_init_btori(initb) f=eeq_btori(pbar) endif if(initb == 1) initb = 0 btori=f else btori=leq_btori(pbar) endif r_last = r theta_last = thet I_last = btori end function btori function dbtori(r, thet) ! returns dI/dpsi use veq, only: veq_dbtori => dbtori , vm_init_dbtori => initialize_dbtori use geq, only: geq_dbtori => dbtori , initialize_dbtori use peq, only: peq_dbtori => dbtori , ppl_init_dbtori => initialize_dbtori use leq, only: leq_dbtori => dbtori use eeq, only: eeq_dbtori => dbtori, efit_init_dbtori => initialize_dbtori real :: dbtori real, intent (in) :: r, thet real :: pbar, f integer :: i, initdb = 1 if(eqinit == 1) initdb = 1 if(iflux == 1) then pbar=min(1.,max(0.,(rpfun(r,thet)-psi_0)/(psi_a-psi_0))) if(vmom_eq) then i=vm_init_dbtori(initdb) f=veq_dbtori(pbar) elseif (gen_eq) then i=initialize_dbtori(initdb) f=geq_dbtori(pbar) elseif (ppl_eq) then i=ppl_init_dbtori(initdb) f=peq_dbtori(pbar) elseif (efit_eq) then i = efit_init_dbtori(initdb) f=eeq_dbtori(pbar) endif if(initdb == 1) initdb = 0 dbtori=f else dbtori=leq_dbtori(pbar) endif end function dbtori function iofrho(rho) use geq, only: geq_iofpbar => btori use peq, only: peq_iofpbar => btori use veq, only: veq_iofpbar => btori use eeq, only: eeq_iofpbar => btori use leq, only: leq_i => btori real :: iofrho, f real, intent (in) :: rho if(vmom_eq) iofrho = veq_iofpbar(pbarofrho(rho)) if (gen_eq) iofrho = geq_iofpbar(pbarofrho(rho)) if (ppl_eq) iofrho = peq_iofpbar(pbarofrho(rho)) if(efit_eq) iofrho = eeq_iofpbar(pbarofrho(rho)) if(local_eq) iofrho = leq_i(rho) iofrho = abs(iofrho) end function iofrho subroutine integrate(arg, grid, ans, n) real, dimension(-ntgrid:), intent (in) :: arg, grid real, dimension(-ntgrid:), intent (out) :: ans integer :: i, n ans=0. do i=1,n ans(i)=ans(i-1)+0.5*(grid(i)-grid(i-1))*(arg(i)+arg(i-1)) enddo do i=-1,-n,-1 ans(i)=ans(i+1)+0.5*(grid(i)-grid(i+1))*(arg(i)+arg(i+1)) enddo end subroutine integrate subroutine rmajortgrid(rgrid, theta, rmajor) real, dimension(-ntgrid:), intent(in) :: rgrid, theta real, dimension(-ntgrid:), intent(out) :: rmajor integer :: i do i=-nth,nth rmajor(i)=Rpos(rgrid(i),theta(i)) enddo return end subroutine rmajortgrid function qfun(pbar) use geq, only: geq_qfun => qfun, geq_init_q => initialize_q use peq, only: peq_qfun => qfun, ppl_init_q => initialize_q use veq, only: veq_qfun => qfun, vm_init_q => initialize_q use eeq, only: eeq_qfun => qfun, efit_init_q => initialize_q use leq, only: leq_qfun => qfun integer :: i, initq = 1 real :: pbar, f, qfun ! if(in_nt .and. iflux == 1) then ! write(*,*) 'qfun not update in nt.' !stop !endif if(iflux == 0 .or. iflux == 2) then qfun=leq_qfun(pbar) return endif if(eqinit ==1 ) initq = 1 if(vmom_eq) then i = vm_init_q(initq) qfun = veq_qfun(pbar) elseif (gen_eq) then i = geq_init_q(initq) qfun = geq_qfun(pbar) elseif (ppl_eq) then i = ppl_init_q(initq) qfun = peq_qfun(pbar) elseif (efit_eq) then i = efit_init_q(initq) qfun = eeq_qfun(pbar) endif if(initq == 1) initq = 0 end function qfun function pfun(r,thet) use veq, only: veq_pfun => pfun, vm_init_pressure => initialize_pressure use geq, only: geq_pfun => pfun, geq_init_pressure => initialize_pressure use peq, only: peq_pfun => pfun, ppl_init_pressure => initialize_pressure use eeq, only: eeq_pfun => pfun, efit_init_pressure => initialize_pressure use leq, only: leq_pfun => pfun real, intent (in) :: r, thet real :: pfun, f real pbar integer :: i, initp = 1 if(iflux /= 1) then pfun=leq_pfun(0.) return endif f=0. pbar=min(1.,max(0.,(psi(r,thet)-psi_0)/(psi_a-psi_0))) if(eqinit == 1) initp = 1 if(vmom_eq) then i=vm_init_pressure(initp) pfun = veq_pfun(pbar) write(*,*) 'check units of p! pfun' elseif (gen_eq) then i=geq_init_pressure(initp) pfun = geq_pfun(pbar) elseif (ppl_eq) then i=ppl_init_pressure(initp) pfun = peq_pfun(pbar) elseif (efit_eq) then i = efit_init_pressure(initp) pfun = eeq_pfun(pbar) endif if(initp == 1) initp = 0 end function pfun function dpfun(r, thet) use veq, only: veq_dpfun => dpfun, vm_init_dpressure => initialize_dpressure use geq, only: geq_dpfun => dpfun, geq_init_dpressure => initialize_dpressure use peq, only: peq_dpfun => dpfun, ppl_init_dpressure => initialize_dpressure use eeq, only: eeq_dpfun => dpfun, efit_init_dpressure => initialize_dpressure use leq, only: leq_dpfun => dpfun real, intent (in) :: r, thet real :: dpfun, f real pbar integer :: i, initdp = 1 if(iflux /= 1) then dpfun=leq_dpfun(0.) return endif f=0. pbar=min(1.,max(0.,(psi(r,thet)-psi_0)/(psi_a-psi_0))) if(eqinit == 1) initdp = 1 if(vmom_eq) i = vm_init_dpressure(initdp) if(gen_eq) i = geq_init_dpressure(initdp) if(ppl_eq) i = ppl_init_dpressure(initdp) if(efit_eq) i = efit_init_dpressure(initdp) initdp=0 if(vmom_eq) dpfun = veq_dpfun(pbar) if(gen_eq) dpfun = geq_dpfun(pbar) if(ppl_eq) dpfun = peq_dpfun(pbar) if(efit_eq) dpfun = eeq_dpfun(pbar) end function dpfun function beta_of_rho(rho) use geq, only: geq_beta => betafun, geq_init_beta => initialize_beta use peq, only: peq_beta => betafun, ppl_init_beta => initialize_beta use veq, only: veq_beta => betafun, vm_init_beta => initialize_beta use eeq, only: eeq_beta => betafun, efit_init_beta => initialize_beta real, intent(in) :: rho real :: beta_of_rho, f integer :: i, initbeta = 1 if(iflux == 0 .or. iflux == 2) then beta_of_rho = 0. return endif if(in_nt .and. iflux == 1) then if(eqinit == 1) initbeta = 1 endif if(vmom_eq) then i = vm_init_beta(initbeta) f = veq_beta(pbarofrho(rho)) elseif (gen_eq) then i = geq_init_beta(initbeta) f = geq_beta(pbarofrho(rho)) elseif (ppl_eq) then i = ppl_init_beta(initbeta) f = peq_beta(pbarofrho(rho)) elseif(efit_eq) then i = efit_init_beta(initbeta) f = eeq_beta(pbarofrho(rho)) endif initbeta = 0 beta_of_rho=f end function beta_of_rho function beta_a_fun(r,thet) use geq, only: geq_beta => betafun, geq_init_beta => initialize_beta use peq, only: peq_beta => betafun, ppl_init_beta => initialize_beta use veq, only: veq_beta => betafun, vm_init_beta => initialize_beta use eeq, only: eeq_beta => betafun, efit_init_beta => initialize_beta real, intent (in) :: r, thet real :: beta_a_fun real :: pbar, f integer :: i, initbeta = 1 if(iflux == 0 .or. iflux == 2) then beta_a_fun=0. return endif if(in_nt .and. iflux == 1) then if(eqinit == 1) initbeta = 1 endif f=0. pbar=min(1.,max(0.,(psi(r,thet)-psi_0)/(psi_a-psi_0))) if(vmom_eq) then i = vm_init_beta(initbeta) f = veq_beta(pbar) elseif (gen_eq) then i = geq_init_beta(initbeta) f = geq_beta(pbar) elseif (ppl_eq) then i = ppl_init_beta(initbeta) f = peq_beta(pbar) elseif(efit_eq) then i = efit_init_beta(initbeta) f = eeq_beta(pbar) endif initbeta = 0 beta_a_fun=f end function beta_a_fun function pbarofrho(rho) real, intent (in) :: rho real :: pbarofrho ! punt for now; not used if(iflux /= 1) then ! write(*,*) 'error in pbarofrho!' pbarofrho=0. return endif pbarofrho=min(1.,max(0.,(rpofrho(rho)-psi_0)/(psi_a-psi_0))) end function pbarofrho function dqdrhofun(rho) real, intent (in) :: rho real dqdrhofun dqdrhofun=shat*qfun(pbarofrho(rho))/(rhoc) end function dqdrhofun subroutine root(f,fval,a,b,xerrbi,xerrsec,nsolv,ier,soln) real, external :: f real :: fval,a,b,a1,b1,f1,f2,f3,trial,xerrbi,xerrsec,soln,aold integer :: i,ier,nsolv,niter,isolv ier=0 a1=a b1=b f1=f(a1)-fval f2=f(b1)-fval if(xerrbi < 0.) goto 1000 if(f1*f2 > 0.) then write(11,*) 'f1 and f2 have same sign in bisection routine' write(11,*) 'a1,f1,b1,f2=',a1,f1,b1,f2 write(11,*) 'fval=',fval ier=1 goto 1000 endif niter=int(alog(abs(b-a)/xerrbi)/alog(2.))+1 do i=1,niter trial=0.5*(a1+b1) f3=f(trial)-fval if(f3*f1 > 0.) then a1=trial f1=f3 else b1=trial f2=f3 endif ! write(11,*) 'i,a1,f1,b1,f2 ',i,a1,f1,b1,f2 enddo 1000 continue if( abs(f1) > abs(f2) ) then f3=f1 f1=f2 f2=f3 aold=a1 a1=b1 b1=aold endif ! start secant method ! write(11,*) 'a1,f1,b1,f2',a1,f1,b1,f2 isolv=0 do i=1,10 aold=a1 f3=f2-f1 if (abs(f3) < 1.e-11) f3=1.e-11 a1=a1-f1*(b1-a1)/f3 f2=f1 b1=aold ! write(11,*) 'a1,f1,b1,f2',a1,f1,b1,f2 if(abs(a1-b1) < xerrsec) isolv=isolv+1 if (isolv >= nsolv) goto 9000 f1=f(a1)-fval enddo 9000 soln=a1 end subroutine root function phi(rp) integer, parameter :: nimax = 200 real :: phi real, intent (in) :: rp real :: pbar, pb(nimax), dpb integer :: i, ni ni=200 if(ni > nimax) then write(*,*) 'Increase nimax in phi function' stop endif pbar=(rp-psi_0)/(psi_a-psi_0) if(.not.between(pbar,0.,1.)) then write(*,*) 'Defining Phi outside of plasma' write(*,*) 'rp, psi_0, psi_a = ',rp,psi_0,psi_a endif phi=0. do i=1,ni pb(i)=pbar*float(i-1)/float(ni-1) enddo dpb=pbar/float(ni-1) do i=1,ni-1 ! fixed by MB phi=phi+0.5*dpb*(qfun(pb(i+1))+qfun(pb(i))) enddo phi=phi*(psi_a-psi_0) ! fixed by MB end function phi function psi(r, thet) use geq, only: geq_psi => psi, geq_init_psi => initialize_psi use peq, only: peq_psi => psi, ppl_init_psi => initialize_psi use veq, only: veq_psi => psi, vm_init_psi => initialize_psi use eeq, only: eeq_psi => psi use leq, only: leq_psi => psi real :: psi real, intent (in) :: r, thet integer :: init = 1, i if(r == 0.) then psi=psi_0 return endif if(r == 1. .and. thet == 0.) then psi=psi_a return endif ! ! Reinitialize if equilibrium has changed. ! if(eqinit == 1) init = 1 if(vmom_eq) then i = vm_init_psi(init) psi = veq_psi(r, thet) elseif (gen_eq) then i = geq_init_psi(init) psi = geq_psi(r, thet) elseif (ppl_eq) then i = ppl_init_psi(init) psi = peq_psi(r, thet) else if(efit_eq) then psi = eeq_psi(r, thet) else if(local_eq) then psi = leq_psi(r, thet) endif init = 0 end function psi function between(x,x1,x2) logical :: between real :: x, x1, x2 if( ((x1 <= x).and.(x <= x2)) .or. ((x1 >= x).and.(x >= x2)) ) then between = .true. else between = .false. end if end function between subroutine geofax(rho, t, e, d) real, intent (in) :: rho real, intent (out) :: t, e, d real :: a, b real, dimension(-ntgrid:ntgrid) :: rgrid real :: rp, hmaxu, rmax integer :: i rp=rpofrho(rho) if(efit_eq) then call rtg(rgrid, rp) else rgrid = rp endif ! ignore up-down asymmetry for now hmaxu = 0. rmax = 0. do i=0,nth a = Rpos(rgrid(i), theta(i)) b = Zpos(rgrid(i), theta(i)) if(hmaxu < b ) then hmaxu = b rmax = a endif enddo a = 0.5*(Rpos(rgrid(0), theta(0))-Rpos(rgrid(nth), theta(nth))) e = hmaxu/a t = asin(rcenter(rp)-rmax) d = rcenter(rp) - rcenter(psi_a) end subroutine geofax function rmagaxis(rp) real :: rmagaxis real, intent (in) :: rp if(iflux == 1) then rmagaxis=rmaj else rmagaxis=R_geo endif end function rmagaxis function rcenter(rp) use geq, only: geq_rcenter => rcenter, geq_init_rc => initialize_rcenter use peq, only: peq_rcenter => rcenter, peq_init_rc => initialize_rcenter use veq, only: veq_rcenter => rcenter, veq_init_rc => initialize_rcenter use eeq, only: eeq_rcenter => rcenter, eeq_init_rc => initialize_bound use leq, only: leq_rcenter => rcenter real, intent (in) :: rp real :: rcenter real :: broot, pi integer :: i, init_rc = 1 pi = 2.*acos(0.) if(eqinit == 1) init_rc = 1 if (gen_eq) i = geq_init_rc(init_rc) if (ppl_eq) i = peq_init_rc(init_rc) if (vmom_eq) i = veq_init_rc(init_rc) if (efit_eq) i = eeq_init_rc(init_rc) init_rc = 0 rcenter=rmaj if(vmom_eq) rcenter = veq_rcenter(rp) if(gen_eq) rcenter = geq_rcenter(rp) if(ppl_eq) rcenter = peq_rcenter(rp) if(efit_eq) rcenter = eeq_rcenter(rp) if(local_eq) rcenter = leq_rcenter(rp) end function rcenter function bmodfun(r,thet) use veq, only: veqitem => eqitem, veqB_psi => B_psi use geq, only: eqitem, eqB_psi => B_psi real :: bmodfun real, intent (in) :: r, thet integer :: init = 1 real :: f if(eqinit == 1) init = 1 if(vmom_eq) then call veqitem(r, thet, veqB_psi, f, 'R') bmodfun=f return elseif (gen_eq) then call eqitem(r, thet, eqB_psi, f, 'R') bmodfun=f return else write(*,*) 'Stopping in bmodfun.' write(*,*) 'You must use gen_eq or vmom_eq to call bmodfun.' stop endif 1000 format(1x,11e16.9) end function bmodfun subroutine smoothie(n, a, b, jmax1, jmax2, f1) ! only works for periodic functions s.t. a(1)=a(n) integer :: n, jmax1, jmax2, i, j real, dimension(n) :: a, b, c, d, e real :: f1, dmax ! ! get weighting function for diffusion coefficient d(1)=a(2)-a(n-1) do i=2,n-1 d(i)=a(i+1)-a(i-1) enddo d(n)=a(2)-a(n-1) do j=1,jmax1/2 c(1)=(d(n-1)+d(1)+d(2))/3. do i=2,n c(i)=(d(i-1)+d(i)+d(i+1))/3. enddo c(n)=c(1) do i=1,n d(i)=c(i) enddo enddo do j=1,jmax1/2 c(1)=abs(d(n-1)+d(1)+d(2))/3. do i=2,n c(i)=abs(d(i-1)+d(i)+d(i+1))/3. enddo c(n)=c(1) do i=1,n d(i)=c(i) enddo enddo dmax=0. do i=1,n ! write(*,*) d(i) d(i)=1.0/d(i)**f1 dmax=max(dmax,d(i)) e(i)=d(i) enddo ! write(*,*) 'dmax ',dmax do j=1,jmax1/2 c(1)=e(1)+0.25*d(1)/dmax*(e(2)-2*e(1)+e(n-1)) do i=2,n-1 c(i)=e(i)+0.25*d(i)/dmax*(e(i+1)-2.*e(i)+e(i-1)) enddo c(n)=c(1) do i=1,n e(i)=c(i) enddo enddo do i=1,n d(i)=e(i) enddo ! do i=1,n ! write(*,*) d(i)/dmax ! enddo ! write(*,*) do j=1,jmax2 c(1)=b(1)+0.25*d(1)/dmax*(b(2)-2*b(1)+b(n-1)) do i=2,n-1 c(i)=b(i)+0.25*d(i)/dmax*(b(i+1)-2.*b(i)+b(i-1)) enddo c(n)=c(1) do i=1,n b(i)=c(i) enddo enddo end subroutine smoothie subroutine arclength (ntheta, nperiod, gpar, arcl) integer, intent (in) :: ntheta, nperiod ! real, dimension(-ntgrid:), intent (in) :: theta real, dimension(-ntgrid:), intent (out) :: arcl real, dimension(-ntgrid:), intent (in out) :: gpar integer :: nth, j, k real :: pi, arcfac nth=ntheta/2 if(2*nth /= ntheta) write(*,*) 'ntheta should be even. ',nth, ntheta pi=2.*acos(0.) arcl(-nth)=0. do j=-nth,nth-1 arcl(j+1)=arcl(j)+0.5*(theta(j+1)-theta(j))*(1./gpar(j)+1./gpar(j+1)) enddo arcfac=1./arcl(nth) do j=-nth,nth arcl(j)=arcl(j)*2.*pi*arcfac-pi enddo if(nperiod > 1) call periodic_copy(arcl, 2.*pi) do k=-nperiod+1,nperiod-1 do j=-nth,nth gpar(j+k*ntheta)=2.*pi*arcfac enddo enddo end subroutine arclength subroutine loftheta(rgrid, theta, ltheta) real, dimension (-ntgrid:) :: rgrid, theta, ltheta integer :: i, j !, itot, k real :: rold, told, rpos_old, zpos_old, rpos_new, zpos_new real :: r, thet !, L_tot ltheta(0)=0. ! ! this could be made more efficient ! rold=rgrid(0) told=theta(0) do i=1,nth r=rgrid(i) thet=theta(i) rpos_old = Rpos(rold, told) zpos_old = Zpos(rold, told) rpos_new = Rpos(r, thet) zpos_new = Zpos(r, thet) ltheta(i)=ltheta(i-1) & +sqrt((Rpos_new-Rpos_old)**2 + (Zpos_new-Zpos_old)**2) rold=r told=thet enddo rold=rgrid(0) told=theta(0) do i=-1,-nth,-1 r=rgrid(i) thet=theta(i) rpos_old = Rpos(rold, told) zpos_old = Zpos(rold, told) rpos_new = Rpos(r, thet) zpos_new = Zpos(r, thet) ltheta(i)=ltheta(i+1) & -sqrt((Rpos_new-Rpos_old)**2 + (Zpos_new-Zpos_old)**2) rold=r told=thet enddo end subroutine loftheta subroutine gradl (ltheta, f, dfdl, ext, n) real, dimension (-ntgrid:), intent(in) :: ltheta, f real, dimension (-ntgrid:), intent(out) :: dfdl real, intent(in) :: ext integer :: n real :: pi integer :: i dfdl(-n) = (f(n-1) + ext - f(-n+1)) & /(ltheta(-n) - ltheta(-n+1))/2. do i = -n+1, n-1 dfdl(i) = (f(i+1) - f(i-1))/(ltheta(i+1)-ltheta(i-1)) enddo dfdl(n) = dfdl(-n) end subroutine gradl subroutine th_bishop(rpgrad, th_bish, nth) integer, intent (in) :: nth real, dimension (-ntgrid:, :), intent (in) :: rpgrad real, dimension (-ntgrid:), intent (out) :: th_bish real, dimension (-ntgrid:ntgrid) :: magrp real, dimension (-ntgrid:ntgrid, 2) :: tvec real :: pi integer :: i pi=2*acos(0.) do i=-nth,nth magrp(i)=sqrt(rpgrad(i,1)**2+rpgrad(i,2)**2) enddo ! tvec is the tangent vector in the surface of the poloidal plane do i=-nth,nth tvec(i,1)=rpgrad(i,2)/magrp(i) tvec(i,2)=-rpgrad(i,1)/magrp(i) enddo ! need to find grad(R) and dot it with tvec. ! but grad(R) is simple -- it is just (1,0) ! so simply return acos(tvec(:,1)) ! with corrections for left-hand quadrants ! do i=-nth,nth tvec(i,1) = min(1.,max(-1.,tvec(i,1))) if(tvec(i,2) > 0.) then th_bish(i)=2*pi-acos(max(min(tvec(i,1),1.0),-1.0)) else th_bish(i)=acos(max(min(tvec(i,1),1.0),-1.0)) endif enddo end subroutine th_bishop subroutine B_pol(rgrid, theta, rpgrad, Bpolmag, nth) integer, intent (in) :: nth real, dimension (-ntgrid:), intent (in) :: rgrid, theta real, dimension (-ntgrid:, :), intent (in) :: rpgrad real, dimension (-ntgrid:) , intent (out) :: Bpolmag real, dimension(-ntgrid:ntgrid) :: Rinv integer :: i do i=-nth,nth Rinv(i) = invRfun(rgrid(i),theta(i)) enddo do i=-nth,nth Bpolmag(i)=sqrt(rpgrad(i,1)**2 + rpgrad(i,2)**2)*Rinv(i) enddo end subroutine B_pol subroutine B_mod(rgrid, theta, Bpolmag, Bmod, nth) integer, intent (in) :: nth real, dimension (-ntgrid:), intent (in) :: rgrid, theta, Bpolmag real, dimension(-ntgrid:), intent (out) :: Bmod real :: bi integer :: i bi = btori(rgrid(0),theta(0)) do i=-nth,nth Bmod(i)=sqrt((bi*invRfun(rgrid(i),theta(i)))**2 & +Bpolmag(i)**2) enddo end subroutine B_mod subroutine R_pol(theta, th_bish, ltheta, Rpol, nth) integer, intent (in) :: nth real, dimension(-ntgrid:), intent (in) :: th_bish, theta, ltheta real, dimension(-ntgrid:), intent (out) :: Rpol real, dimension(-ntgrid:ntgrid) :: dthdl real :: pi integer :: i, is pi=2.*acos(0.) ! ! R = 1/(d theta/dl * d th_bish/d theta) ! ! the only tricky part is near th_bish = 0. ! ! find this point: is = nth-1 do i=1,nth-1 if(abs(th_bish(i+1)-th_bish(i)) > pi) then is=i goto 10 endif enddo 10 continue do i=1,is-1 Rpol(i)=(th_bish(i+1)-th_bish(i-1))/(theta(i+1)-theta(i-1)) enddo Rpol(is)=(th_bish(is+1)-2*pi-th_bish(is-1))/(theta(is+1)-theta(is-1)) Rpol(is+1)=(th_bish(is+2)-th_bish(is)-2*pi)/(theta(is+2)-theta(is)) do i=is+2,nth-1 Rpol(i)=(th_bish(i+1)-th_bish(i-1))/(theta(i+1)-theta(i-1)) enddo Rpol(nth)=(th_bish(-nth+1)-th_bish(nth-1)) & /(theta(-nth+1)-theta(nth-1)+2*pi) Rpol(-nth)=Rpol(nth) do i=-nth+1,0 Rpol(i)=(th_bish(i+1)-th_bish(i-1))/(theta(i+1)-theta(i-1)) enddo dthdl(-nth) = (theta(nth-1) - theta(-nth+1) -2*pi) & /(ltheta(nth-1) - ltheta(-nth+1) - 2*ltheta(nth)) do i = -nth+1, nth-1 dthdl(i) = (theta(i+1)-theta(i-1))/(ltheta(i+1)-ltheta(i-1)) enddo dthdl(nth) = dthdl(-nth) do i=-nth,nth Rpol(i) = 1/(Rpol(i)*dthdl(i)) enddo end subroutine R_pol subroutine test(rgrid, theta, Bpolmag, Bmod, Rpol, th_bish, bgrad, nth) integer, intent (in) :: nth real, dimension(-ntgrid:), intent (in) :: rgrid, theta, Bmod, Bpolmag, & Rpol, th_bish real, dimension(-ntgrid:, :), intent(in) :: bgrad real, dimension(-ntgrid:ntgrid) :: rmajor, bbgrad real :: b, biggest, pi, dp, bi, di, bp integer :: i pi = 2.*acos(0.) call rmajortgrid(rgrid, theta, rmajor) if(iflux /= 1) write(*,*) 'should not be using test.' dp = dpfun(rgrid(0),theta(0)) bi = btori(rgrid(0),theta(0)) di = dbtori(rgrid(0),theta(0)) do i = -nth, nth bp=-Bpolmag(i) bbgrad(i) = bp/Bmod(i) & *(bp/Rpol(i) + dp*rmajor(i) & - bi**2*sin(th_bish(i))/rmajor(i)**3/bp) ! bbgrad(i) = -di*bi/rmajor(i)*bp & ! -bi**2*sin(th_bish(i))/rmajor(i)**3 enddo do i=-nth,nth bp=-Bpolmag(i) write(*,100) theta(i), & ! bi/rmajor(i)*bgrad(i,1), & bgrad(i,1), & bbgrad(i), & +Bp**2/Bmod(i)/Rpol(i), & +Bp/Bmod(i)*dp*rmajor(i), & -1./Bmod(i)*bi**2*sin(th_bish(i))/rmajor(i)**3, & ! Rpol(i), & ! th_bish(i), & ! -di*bi/rmajor(i)*Bp, & ! -bi**2*sin(th_bish(i))/rmajor(i)**3, & bbgrad(i) - bgrad(i,1) ! Bmod(i), bgrad(i,2) enddo 100 format(20(1x,g12.6)) end subroutine test subroutine tdef(nthg) real :: pi integer :: nthg, nthsave, i ! logical :: first = .true. ! if(.not.first) return ! first = .false. pi = 2*acos(0.) nthsave=nth nth=nthg-1 ! write(*,*) 'old nth: ',nthsave,' new nth: ',nth ntheta=2*nth ntgrid=(2*nperiod-1)*nth ! redefine theta grid used by the rest of the code. ! note that the gen_eq theta grid has theta(1)=0. ! problems avoided by mtheta grid in eqitem. if(.not.allocated(theta)) allocate(theta(-ntgrid:ntgrid)) theta(0)=0. do i=1,nth theta(i)=i*pi/float(nth) theta(-i)=-theta(i) enddo end subroutine tdef subroutine alloc_module_arrays(n) integer n allocate(grho (-n:n), & bmag (-n:n), & gradpar (-n:n), & cvdrift (-n:n), & cvdrift0 (-n:n), & gbdrift (-n:n), & gbdrift0 (-n:n), & gds2 (-n:n), & gds21 (-n:n), & gds22 (-n:n), & jacob (-n:n), & iob2p (-n:n), & job (-n:n), & jobp (-n:n) ) end subroutine alloc_module_arrays subroutine init_theta(nt) integer, intent(in) :: nt integer i real :: pi pi = 2*acos(0.) ntheta=nt nth = nt / 2 ntgrid = (2*nperiod - 1)*nth if(.not.allocated(theta)) allocate(theta(-ntgrid:ntgrid)) theta = (/ (i*pi/real(nth), i=-ntgrid, ntgrid ) /) end subroutine init_theta subroutine periodic_copy(a, ext) real, dimension(-ntgrid:ntgrid) :: a real ext integer :: i, k, itot if(nperiod == 1) return do k=1-nperiod,-1 do i=-nth,nth itot = i + k*ntheta a(itot) = a(i) + k*ext enddo enddo do k=1,nperiod-1 do i=-nth,nth itot = i + k*ntheta a(itot) = a(i) + k*ext enddo enddo end subroutine periodic_copy subroutine bishop_gradB(rgrid, Bmod, Bpolmag, Rpol, th_bish, ltheta, & gradB) real, dimension(-ntgrid:), intent (in) :: rgrid, Bmod, Bpolmag, & Rpol, th_bish, ltheta real, dimension(-ntgrid:,:), intent(out) :: gradB real, dimension(-ntgrid:ntgrid) :: rmajor, dBdl real :: bi, bp integer :: i call rmajortgrid(rgrid, theta, rmajor) call gradl(ltheta, Bmod, dBdl, 0., nth) bi = btori(rgrid(0),theta(0)) do i=-nth, nth bp=-Bpolmag(i) gradB(i,1)=bp/Bmod(i) *(bp/Rpol(i) + dp_new*rmajor(i) & - bi**2*sin(th_bish(i))/rmajor(i)**3/bp) gradB(i,2)=dBdl(i) enddo end subroutine bishop_gradB subroutine check(veq, geq, eeq, peq, leq) logical, intent(in) :: veq, geq, eeq, peq, leq if(veq .and. geq) then write(*,*) 'Choosing vmom_eq = .true. AND gen_eq = .true. is not permitted.' write(*,*) 'Stopping.' stop endif if(veq .and. eeq) then write(*,*) 'Choosing vmom_eq = .true. AND efit_eq = .true. is not permitted.' write(*,*) 'Stopping.' stop endif if(veq .and. leq) then write(*,*) 'Choosing vmom_eq = .true. AND iflux = 0 is not permitted.' write(*,*) 'Stopping.' stop endif if(veq .and. peq) then write(*,*) 'Choosing vmom_eq = .true. AND ppl_eq = .true. is not permitted.' write(*,*) 'Stopping.' stop endif if(geq .and. eeq) then write(*,*) 'Choosing gen_eq = .true. AND efit_eq = .true. is not permitted.' write(*,*) 'Stopping.' stop endif if(geq .and. peq) then write(*,*) 'Choosing gen_eq = .true. AND ppl_eq = .true. is not permitted.' write(*,*) 'Stopping.' stop endif if(geq .and. leq) then write(*,*) 'Choosing gen_eq = .true. AND iflux = 0 is not permitted.' write(*,*) 'Stopping.' stop endif if(eeq .and. leq) then write(*,*) 'Choosing efit_eq = .true. AND iflux = 0 is not permitted.' write(*,*) 'Stopping.' stop endif if(eeq .and. peq) then write(*,*) 'Choosing efit_eq = .true. AND ppl_eq = .true. is not permitted.' write(*,*) 'Stopping.' stop endif if(peq .and. leq) then write(*,*) 'Choosing ppl_eq = .true. AND iflux = 0 is not permitted.' write(*,*) 'Stopping.' stop endif end subroutine check subroutine grad(rgrid, theta, gradf, char, rp, nth, ntgrid) use geq, only: geq_gradient => gradient use peq, only: peq_gradient => gradient use veq, only: veq_gradient => gradient use eeq, only: eeq_gradient => gradient use leq, only: leq_gradient => gradient integer :: nth, ntgrid real, dimension(-ntgrid:) :: rgrid, theta real, dimension(-ntgrid:,:) :: gradf character*1 :: char real rp if(gen_eq) call geq_gradient(rgrid, theta, gradf, char, rp, nth, ntgrid) if(ppl_eq) call peq_gradient(rgrid, theta, gradf, char, rp, nth, ntgrid) if(vmom_eq) call veq_gradient(rgrid, theta, gradf, char, rp, nth, ntgrid) if(efit_eq) call eeq_gradient(rgrid, theta, gradf, char, rp, nth, ntgrid) if(local_eq) call leq_gradient(rgrid, theta, gradf, char, rp, nth, ntgrid) end subroutine grad subroutine bgrad(rgrid, theta, gradf, char, rp, nth, ntgrid) use geq, only: geq_bgradient => bgradient use peq, only: peq_bgradient => bgradient use veq, only: veq_bgradient => bgradient use eeq, only: eeq_bgradient => bgradient use leq, only: leq_bgradient => bgradient integer :: nth, ntgrid real, dimension(-ntgrid:) :: rgrid, theta real, dimension(-ntgrid:,:) :: gradf character*1 :: char real rp if(gen_eq) call geq_bgradient(rgrid, theta, gradf, char, rp, nth, ntgrid) if(ppl_eq) call peq_bgradient(rgrid, theta, gradf, char, rp, nth, ntgrid) if(vmom_eq) call veq_bgradient(rgrid, theta, gradf, char, rp, nth, ntgrid) if(efit_eq) call eeq_bgradient(rgrid, theta, gradf, char, rp, nth, ntgrid) if(local_eq) call leq_bgradient(rgrid, theta, gradf, char, rp, nth, ntgrid) end subroutine bgrad subroutine rtg(rgrid, rp) use eeq, only: bound, rfun real, dimension(-ntgrid:), intent (out) :: rgrid real, intent (in) :: rp real :: broot integer :: i if(.not. efit_eq) then write(*,*) 'error in rtg. efit_eq = .false.' stop endif do i=-nth,nth rgrid(i)=rfun(rp, theta(i), bound(theta(i))) enddo end subroutine rtg subroutine Hahm_Burrell(i, a, rho_star0) use geq, only: geq_Hahm_Burrell => Hahm_Burrell, phipp, temp use peq, only: peq_Hahm_Burrell => Hahm_Burrell use eeq, only: eeq_Hahm_Burrell => Hahm_Burrell integer i real :: a, fract, squeeze, phi_pp, rho_star0 real :: pbar, temp_local, rp_loc, qpar, Lp real :: rhocm, rp_locm, qparm, Lpm real :: rhocp, rp_locp, qparp, Lpp, g1p real :: temp_localp, temp_localm real :: phi_ppp, phi_ppm, pbarm, pbarp real :: squeezep, squeezem if(gen_eq) call geq_Hahm_Burrell(i, a, rho_star0) if(efit_eq) call eeq_Hahm_Burrell(i, a, rho_star0) if(ppl_eq) call peq_Hahm_Burrell(i, a, rho_star0) if(.not. gen_eq) return ! calculate heat flux stuff fract = f_trap(bmag(-nth:nth)) rhocm = rhoc*(1.-delrho) rhocp = rhoc*(1.+delrho) pbarm = pbarofrho(rhocm) pbar = pbarofrho(rhoc) pbarp = pbarofrho(rhocp) rp_locm = rpofrho(rhocm) rp_loc = rpofrho(rhoc) rp_locp = rpofrho(rhocp) temp_localm = temp(pbarm, a) temp_local = temp(pbar , a) temp_localp = temp(pbarp, a) phi_ppm = phipp(pbarm, a) phi_pp = phipp(pbar , a) phi_ppp = phipp(pbarp, a) squeezem = 1. + rho_star0**2*btori(rp_locm,0.)**2*abs(phi_ppm)*temp_localm squeeze = 1. + rho_star0**2*btori(rp_loc ,0.)**2*abs(phi_pp )*temp_local squeezep = 1. + rho_star0**2*btori(rp_locp,0.)**2*abs(phi_ppp)*temp_localp Lpm = -pfun(rp_locm, 0.)/dpfun(rp_locm, 0.)*drhodpsin Lp = -pfun(rp_loc , 0.)/dpfun(rp_loc , 0.)*drhodpsin Lpp = -pfun(rp_locp, 0.)/dpfun(rp_locp, 0.)*drhodpsin g1p = (1./Lpm-1./Lpp)/(-2.*delrho)*Lp & -1.5*sqrt(2.*squeeze)/1.36*(squeezem-squeezep)/(-2.*delrho) & /(fract+sqrt(2.)/1.36*squeeze**1.5) g1p = g1p*2.5*rho_star0*btori(rp_loc, 0.)*fract & *drhodpsin*sqrt(temp_local) & /(fract+sqrt(2.)/1.36*squeeze**1.5) & *(1.-a)/Lp ! get /(v_t p B_a) qpar = 1./(fract+sqrt(2.)/1.36*squeeze**1.5) qpar = qpar*2.5*rho_star0 qpar = qpar*sqrt(temp_local)*drhodpsin qpar = qpar*btori(rp_loc, 0.)*fract ! This is supposed to be 1/L_TiN factor ! (1.-a)/L_pN = 1/L_TN ! and eventually I need to do further correction to get L_TiN ! when LTi /= LTe qpar = qpar*(1.-a)/Lp ! write(*,*) 'rhoc S T/T_0 qparB 1/L_nN' ! write(*,*) rhoc, squeeze, temp_local, qpar, a/Lp ! write(*,*) drhodpsin, fract, btori(rp_loc, 0.) ! write(*,*) 1. + rho_star0**2*btori(rp_loc ,0.)**2*abs(phi_pp )*temp_local ! write(*,*) 'eta_q ~ ',-Lp/(1.e-15+a)*g1p,' qparN = ',qpar end subroutine Hahm_Burrell function nth_get() integer nth_get nth_get = nth end function nth_get subroutine drho_drp(rp, drhodrp) real :: rp1, rp2, rp, rho1, rho2 real, intent (out) :: drhodrp ! compute d rho / d rp rp1=rp*(1.-delrho) rp2=rp*(1.+delrho) if(rp2 > rpmax) then write(*,*) 'delrho too big: ',rp2, rpmax write(*,*) 'Make delrho smaller.' stop endif if(rp1 < rpmin) then write(*,*) 'delrho too big: ',rp1, rpmin write(*,*) 'Make delrho smaller.' stop endif if(irho == 1) then rho1=sqrt(abs(phi(rp1)/phi(rpmax))) rho2=sqrt(abs(phi(rp2)/phi(rpmax))) elseif(irho == 2) then rho1=0.5*diameter(rp1) rho2=0.5*diameter(rp2) elseif(irho == 3) then rho1=psibarfun(rp1) rho2=psibarfun(rp2) endif drhodrp=(rho2-rho1)/(rp2-rp1) end subroutine drho_drp function fluxavg(f) real, dimension(-nth:nth), intent(in) :: f real, dimension(-nth:nth) :: delth real :: fluxavg delth(-nth+1:nth-1) = 0.5*(theta(-nth+2:nth)-theta(-nth:nth-2)) delth(-nth) = 0.5*(theta(-nth+1)-theta(-nth)) delth(nth) = 0.5*(theta(nth)-theta(nth-1)) fluxavg = sum(f*jacob(-nth:nth)*delth)/sum(jacob(-nth:nth)*delth) end function fluxavg function f_trap(b_mag) real, dimension(-nth:nth), intent (in) :: b_mag real :: f_trap, ftu, ftl, havg, h2avg, h(-nth:nth) ! use expressions from Y. R. Lin-Liu and Bob Miller, coding from Q. P. Liu h = min(b_mag/maxval(b_mag),1.0) havg = fluxavg(h) h2avg = fluxavg(h**2) ftu = 1.0 - h2avg/havg**2*(1.0-sqrt(1.0-havg)*(1.0+0.5*havg)) ftl = 1.0 - h2avg*fluxavg((1.0-sqrt(1.0-h)*(1.0+0.5*h))/h**2) f_trap = 0.75*ftu + 0.25*ftl end function f_trap function hirsh_bs(b_mag, Z_eff, gamma) real, dimension(-nth:nth), intent (in) :: b_mag real :: hirsh_bs, Z_eff, gamma, ft, x, L31, L32, alfa, denom, b2 ft = f_trap(b_mag) x = ft/(1.-ft) denom = 1.414*Z_eff + Z_eff**2 & + x*(0.754 + 2.657*Z_eff + 2.*Z_eff**2) & + x**2*(0.348 + 1.243*Z_eff + Z_eff**2) L31 = x/denom * ((0.754 + 2.21*Z_eff + Z_eff**2) & + x*(0.348 + 1.243*Z_eff + Z_eff**2)) L32 = -x/denom * (0.884 + 2.074*Z_eff) alfa = -1.172/(1. + 0.462*x) hirsh_bs = L31*(1. + alfa*gamma/2.) + L32 * gamma/2. ! We really want J_par / B. We know that J_BS/B = C_0, i.e., ind. of ! theta because div J = 0 only allows that possibility. hirsh_bs is ! a formula for , which is manifestly independent of theta. ! ==> = J_BS/B ! ==> J_BS/B = hirsh_bs / ! b2 = fluxavg(b_mag**2) ! hirsh_bs = hirsh_bs / b2 ! Multiply by remaining factor of I*p' in calling routine end function hirsh_bs end module geometry