module eeq implicit none private integer :: nw, nh, nbbbs real, allocatable, dimension (:) :: psi_bar, fp, qsf, pressure, beta, spsi_bar real, allocatable, dimension (:) :: dummy, efit_R, efit_Z, sefit_R, sefit_Z real, allocatable, dimension (:,:) :: efit_psi, sefit_psi real, allocatable, dimension (:,:,:) :: dpm, dtm real, allocatable, dimension (:) :: rbbbs, zbbbs, thetab, r_bound !boundary of plasma real :: efit_dR, efit_dZ, psi_N real :: R_mag, Z_mag, B_T0, aminor, beta_0, psi0, psia logical :: init_bound = .true. logical :: init_btori = .true. logical :: init_dbtori = .true. logical :: init_q = .true. logical :: init_pressure = .true. logical :: init_dpressure = .true. logical :: init_beta = .true. logical :: init_R = .true. logical :: init_Z = .true. public :: efit_init, efitin, gradient, eqitem, bgradient, rfun, diameter, rcenter, & Hahm_Burrell public :: invR public :: Rpos public :: Zpos public :: btori, initialize_btori public :: dbtori, initialize_dbtori public :: qfun, initialize_q public :: pfun, initialize_pressure public :: dpfun, initialize_dpressure public :: betafun, initialize_beta public :: bound, initialize_bound public :: psi contains subroutine efitin(eqfile, psi_0, psi_a, rmaj, B_T, amin, initeq, big) ! ! This subroutine reads an EFIT output file containing ! the axisymmetric magnetic field geometry on a rectangular ! domain defined by the coordinates (R,Z). ! ! efit_R R grid ! efit_Z Z grid ! fp F on psibar grid ! efit_psi psi on (R,Z) grid ! R_mag major radius of the magnetic axis ! Z_mag elevation of the magnetic axis ! rwid total width of the domain ! rleft position of leftmost point of domain ! zhei total height of domain ! implicit none real :: xdum, p_0 real :: rwid, rleft, zhei, amin, B_T real :: psi_0, psi_a, rmaj, rcentr, bcentr real, dimension(:), allocatable :: tmp1, tmp2, zp, temp, & zx1, zxm, zy1, zyn real :: zxy11, zxym1, zxy1n, zxymn real :: fitp_surf2 character*80 :: filename, eqfile character char*10 integer :: i, j, init, ndum, initeq, i1, big, nhb, nwb, ierr data init /1/ save init ! Need to generalize initialization condition if equilibrium changes if(initeq == 0) return init=0 i=index(eqfile,' ')-1 filename = eqfile(1:i) open(unit=5,file=filename,status='old',form='formatted') ! Read the data read(5,1000) char, char, char, char, char, i, nw, nh nwb = nw * big nhb = nh * big call alloc_module_arrays(nwb, nwb, nhb, nw, nh) read(5,2020) rwid, zhei, rcentr, rleft, xdum read(5,2020) R_mag, Z_mag, psi_0, psi_a, bcentr ! ! pbar is defined by ! pbar = (psi-psi_0)/(psi_a-psi_0) ! fp and q are functions of pbar ! do i=1,2 read(5,2020)xdum,xdum,xdum,xdum,xdum enddo do i=1,nw spsi_bar(i) = float(i-1) / float(nw-1) sefit_R(i) = rleft +rwid * float(i-1) / float(nw-1) enddo do j=1,nh sefit_Z(j) = ((float(j-1) / float(nh-1) ) - 0.5)*zhei enddo do i=1,nwb psi_bar(i) = float(i-1) / float(nwb-1) efit_R(i) = rleft +rwid * float(i-1) / float(nwb-1) enddo do j=1,nhb efit_Z(j) = ((float(j-1) / float(nhb-1) ) - 0.5)*zhei enddo read(5,2020) (dummy(j), j = 1, nw) call inter_cspl(nw, spsi_bar, dummy, nwb, psi_bar, fp) read(5,2020) (dummy(j), j = 1, nw) call inter_cspl(nw, spsi_bar, dummy, nwb, psi_bar, pressure) read(5,2020) (dummy(j), j = 1, nw) read(5,2020) (dummy(j), j = 1 ,nw) read(5,2020) ((sefit_psi(i,j) , i = 1, nw) , j = 1, nh) allocate(zp(3*nw*nh), temp(nw+2*nh)) allocate(zx1(nh), zxm(nh), zy1(nw), zyn(nw)) call fitp_surf1(nw, nh, sefit_R, sefit_Z, sefit_psi, & nw, zx1, zxm, zy1, zyn, zxy11, zxym1, zxy1n, zxymn, & 255, zp, temp, 1., ierr) do j = 1, nhb do i = 1, nwb efit_psi(i,j) = fitp_surf2(efit_R(i), efit_Z(j), nw, nh, & sefit_R, sefit_Z, sefit_psi, nw, zp, 1.) enddo enddo deallocate(zp, temp) read(5,2020) (dummy(j) , j = 1, nw) call inter_cspl(nw, spsi_bar, dummy, nwb, psi_bar, qsf) nw = nwb nh = nhb read(5,2022) nbbbs, ndum allocate(rbbbs(nbbbs), zbbbs(nbbbs), thetab(nbbbs), r_bound(nbbbs)) allocate(tmp1(nbbbs), tmp2(nbbbs)) read(5,2020) (rbbbs(i), zbbbs(i) , i = 1, nbbbs) ! get r_boundary(theta) thetab = atan2 ((zbbbs-Z_mag), (rbbbs-R_mag)) r_bound = sqrt( (rbbbs - R_mag)**2 + (zbbbs - Z_mag)**2 ) call sort(thetab, r_bound, zbbbs, rbbbs) ! Allow for duplicated points near +- pi: if(thetab(1) == thetab(2)) then thetab(1) = thetab(1) + 4.*acos(0.) call sort(thetab, r_bound, zbbbs, rbbbs) endif if(thetab(nbbbs-1) == thetab(nbbbs)) then thetab(nbbbs) = thetab(nbbbs) - 4.*acos(0.) call sort(thetab, r_bound, zbbbs, rbbbs) endif ! It isn't likely that a duplicate point would exist near theta = 0, ! so I am not allowing this possibility for now. do i=1,nbbbs-1 if(thetab(i) == thetab(i+1)) then ! write(*,*) 'Duplicates near theta = 0 not allowed.' ! write(*,*) i, i+1, ' Stopping.' ! write(*,*) r_bound(i),r_bound(i+1) ! write(*,*) rbbbs(i),rbbbs(i+1) ! write(*,*) zbbbs(i),zbbbs(i+1) ! write(*,*) thetab(i),thetab(i+1) ! stop thetab(i)=thetab(i)*0.9999999999999 endif enddo call a_minor(rbbbs, zbbbs, Z_mag, amin) aminor=amin ! read(5,2020)(xlim(i),ylim(i),i=1,limitr) close(5) deallocate (rbbbs, zbbbs, tmp1, tmp2) r_bound = r_bound / aminor R_mag = R_mag / aminor Z_mag = Z_mag / aminor ! rleft = rleft / aminor ! rwid = rwid / aminor ! zhei = zhei / aminor rcentr = rcentr / aminor efit_R = efit_R / aminor efit_Z = efit_Z / aminor ! should rmaj be R_mag or rcentr? use R_mag for now. rmaj = R_mag B_T0 = abs(bcentr) B_T = B_T0 psi_a = psi_a / (B_T0*aminor**2) psi_0 = psi_0 / (B_T0*aminor**2) psi_N = psi_a - psi_0 psi0 = psi_0 psia = psi_a p_0=pressure(1) ! do i=1,nw ! psi_bar(i) = float(i-1) / float(nw-1) ! efit_R(i) = rleft +rwid * float(i-1) / float(nw-1) ! enddo fp = fp / (B_T0*aminor) ! MKS: beta = 2 mu_0 p / B**2 beta = 8. * (2. * acos(0.)) * pressure * 1.e-7 / B_T0**2 beta_0 = beta(1) pressure = pressure / p_0 efit_psi = efit_psi / (B_T0*aminor**2) ! do j=1,nh ! efit_Z(j) = ((float(j-1) / float(nh-1) ) - 0.5)*zhei ! enddo efit_dR = efit_R(2) - efit_R(1) efit_dZ = efit_Z(2) - efit_Z(1) 1000 format(5(a10),i2,i4,i4) 2020 format (5e16.9) 2022 format (2i5) end subroutine efitin subroutine efit_init real, dimension(nw, nh) :: eqth integer :: i, j do i = 1, nw do j = 1,nh if(efit_Z(j) == Z_mag .and. efit_R(i) == R_mag) then eqth(i,j) = 0. ! value should not matter else eqth(i,j) = atan2( (efit_Z(j)-Z_mag), (efit_R(i)-R_mag)) endif enddo enddo call derm(efit_psi, dpm) call tderm(eqth, dtm) end subroutine efit_init subroutine tderm(f, dfm) implicit none integer i, j real f(:,:), dfm(:,:,:), pi pi = 2.*acos(0.) ! EFIT grid is equally spaced in R, Z -- this routine uses that fact and ! is therefore not completely general. It is fine for EFIT output. i=1 dfm(i,:,1) = -0.5*(3*f(i,:)-4*f(i+1,:)+f(i+2,:))/efit_dR i=nw dfm(i,:,1) = 0.5*(3*f(i,:)-4*f(i-1,:)+f(i-2,:))/efit_dR j=1 dfm(:,j,2) = -0.5*(3*f(:,j)-4*f(:,j+1)+f(:,j+2))/efit_dZ j=nh dfm(:,j,2) = 0.5*(3*f(:,j)-4*f(:,j-1)+f(:,j-2))/efit_dZ do i=2,nw-1 dfm(i,:,1)=0.5*(f(i+1,:)-f(i-1,:))/efit_dR enddo do j=2,nh-1 do i = 1,nw if(f(i,j+1)-f(i,j-1) > pi) then dfm(i,j,2)=0.5*(f(i,j+1)-f(i,j-1)-2.*pi)/efit_dZ else dfm(i,j,2)=0.5*(f(i,j+1)-f(i,j-1))/efit_dZ endif enddo enddo end subroutine tderm subroutine derm(f, dfm) implicit none integer i, j real f(:,:), dfm(:,:,:) ! EFIT grid is equally spaced in R, Z -- this routine uses that fact and ! is therefore not completely general. It is fine for EFIT output. i=1 dfm(i,:,1) = -0.5*(3*f(i,:)-4*f(i+1,:)+f(i+2,:))/efit_dR i=nw dfm(i,:,1) = 0.5*(3*f(i,:)-4*f(i-1,:)+f(i-2,:))/efit_dR j=1 dfm(:,j,2) = -0.5*(3*f(:,j)-4*f(:,j+1)+f(:,j+2))/efit_dZ j=nh dfm(:,j,2) = 0.5*(3*f(:,j)-4*f(:,j-1)+f(:,j-2))/efit_dZ do i=2,nw-1 dfm(i,:,1)=0.5*(f(i+1,:)-f(i-1,:))/efit_dR enddo do j=2,nh-1 dfm(:,j,2)=0.5*(f(:,j+1)-f(:,j-1))/efit_dZ enddo end subroutine derm subroutine gradient(rgrid, theta, grad, char, rp, nth, ntm) integer nth, ntm character*1 char real, dimension(-ntm:), intent(in) :: rgrid, theta real, dimension(-ntm:,:), intent(out) :: grad real tmp1, tmp2, aa, daa, rp, rp_bar integer i, j grad = 0. do i=-nth, nth call eqitem(rgrid(i), theta(i), dpm(:,:,1), grad(i,1)) call eqitem(rgrid(i), theta(i), dpm(:,:,2), grad(i,2)) enddo ! to get grad(pressure), multiply grad(psi) by dpressure/dpsi if(char == 'R') then rp_bar = rp call inter_d_cspl(nw, psi_bar, pressure, 1, rp_bar, aa, daa) grad = grad*daa*0.5* beta_0/psi_N endif end subroutine gradient subroutine bgradient(rgrid, theta, grad, char, rp, nth_used, ntm) implicit none integer nth_used, ntm character*1 char real rgrid(-ntm:), theta(-ntm:), grad(-ntm:,:) real tmp(2), aa, daa, rp, rp_bar real, dimension(nw, nh, 2) :: dbish integer i logical :: first = .true. dbish(:,:,1) = sqrt(dpm(:,:,1)**2 + dpm(:,:,2)**2) dbish(:,:,2) = 0. if(char == 'T') then ! the order of the next two statements matters dbish(:,:,2) = (dtm(:,:,2)*dpm(:,:,1)-dtm(:,:,1)*dpm(:,:,2))/dbish(:,:,1) dbish(:,:,1) = (dtm(:,:,1)*dpm(:,:,1)+dtm(:,:,2)*dpm(:,:,2))/dbish(:,:,1) endif do i=-nth_used,-1 call eqitem(rgrid(i), theta(i), dbish(:,:,1), tmp(1)) call eqitem(rgrid(i), theta(i), dbish(:,:,2), tmp(2)) grad(i,1) = tmp(1) grad(i,2) = tmp(2) enddo do i=0,nth_used call eqitem(rgrid(i), theta(i), dbish(:,:,1), tmp(1)) call eqitem(rgrid(i), theta(i), dbish(:,:,2), tmp(2)) grad(i,1)=tmp(1) grad(i,2)=tmp(2) enddo ! to get grad(pressure), multiply grad(psi) by dpressure/dpsi if(char == 'R') then rp_bar = rp call inter_d_cspl(nw, psi_bar, pressure, 1, rp_bar, aa, daa) do i=-nth_used, nth_used grad(i,1)=grad(i,1)*daa * 0.5*beta_0/psi_N grad(i,2)=grad(i,2)*daa * 0.5*beta_0/psi_N enddo endif end subroutine bgradient subroutine eqitem(r, thetin, f, fstar) integer :: i, j, istar, jstar real, intent (in) :: r, thetin, f(:,:) real, intent (out) :: fstar real st, dt, sr, dr real r_pos, z_pos r_pos = Rpos(r, thetin) z_pos = Zpos(r, thetin) ! find point on R mesh if(r_pos >= efit_R(nw) .or. r_pos <= efit_R(1)) then write(*,*) 'No evaluation of eqitem allowed outside' write(*,*) 'or on edge of R domain' write(*,*) r, thetin, efit_R(nw), r_pos ! call aborter(72,'R problem') stop endif ! ensure point is on Z mesh if(z_pos >= efit_Z(nh) .or. z_pos <= efit_Z(1)) then write(*,*) 'No evaluation of eqitem allowed outside' write(*,*) 'or on edge of Z domain' write(*,*) r, thetin, efit_Z(1), efit_Z(nh), z_pos stop endif istar=0 do i=2,nw if(istar /= 0) cycle if(r_pos < efit_R(i)) then dr = r_pos - efit_R(i-1) sr = efit_R(i) - r_pos istar=i-1 endif enddo ! Now do Z direction jstar=0 do j=1,nh if(jstar /= 0) cycle if(z_pos < efit_Z(j)) then dt = z_pos - efit_Z(j-1) st = efit_Z(j) - z_pos jstar=j-1 endif enddo ! use opposite area stencil to interpolate fstar=f(istar , jstar ) * sr * st & +f(istar + 1, jstar ) * dr * st & +f(istar , jstar + 1) * sr * dt & +f(istar + 1, jstar + 1) * dr * dt fstar = fstar & /abs(efit_R(istar+1)-efit_R(istar)) & /(efit_Z(jstar+1)-efit_Z(jstar)) ! write(*,*) i, dr, dt, sr, st ! write(*,*) f(istar,jstar+1),f(istar+1,jstar+1) ! write(*,*) f(istar,jstar),f(istar+1,jstar) ! write(*,*) efit_R(istar),efit_R(istar+1) ! write(*,*) efit_Z(jstar),efit_Z(jstar+1) end subroutine eqitem function Zpos (r, theta) real, intent (in) :: r, theta real :: Zpos Zpos = Z_mag + r * sin(theta) end function Zpos function Rpos (r, theta) real, intent (in) :: r, theta real :: Rpos Rpos = R_mag + r * cos(theta) end function Rpos function invR (r, theta) real, intent (in) :: r, theta real :: invR invR = 1/(R_mag + r*cos(theta)) end function invR function psi (r, theta) real, intent (in) :: r, theta real :: f, psi if(r == 0.) then psi=psi0 return endif if(r == 1. .and. theta == 0.) then psi=psia return endif call eqitem(r, theta, efit_psi, f) psi = f end function psi function initialize_btori (init) integer :: init, initialize_btori init_btori = .false. if(init == 1) init_btori = .true. initialize_btori = 1 end function initialize_btori function btori (pbar) use splines real :: pbar, btori type (spline), save :: spl if(init_btori) call new_spline(nw, psi_bar, fp, spl) btori = splint(pbar, spl) end function btori function initialize_dbtori (init) integer :: init, initialize_dbtori init_dbtori = .false. if(init == 1) init_dbtori = .true. initialize_dbtori = 1 end function initialize_dbtori function dbtori (pbar) use splines real :: pbar, dbtori type (spline), save :: spl if(init_dbtori) call new_spline(nw, psi_bar, fp, spl) dbtori = dsplint(pbar, spl)/psi_N end function dbtori function initialize_q (init) integer :: init, initialize_q init_q = .false. if(init == 1) init_q = .true. initialize_q = 1 end function initialize_q function qfun (pbar) use splines real :: pbar, qfun type (spline), save :: spl if(init_q) then call new_spline(nw, psi_bar, qsf, spl) endif qfun = splint(pbar, spl) end function qfun function initialize_pressure (init) integer :: init, initialize_pressure init_pressure = .false. if(init == 1) init_pressure = .true. initialize_pressure = 1 end function initialize_pressure function pfun (pbar) use splines real :: pbar, pfun type (spline), save :: spl if(init_pressure) call new_spline(nw, psi_bar, pressure, spl) pfun = splint(pbar, spl) * beta_0/2. end function pfun function initialize_dpressure (init) integer :: init, initialize_dpressure init_dpressure = .false. if(init == 1) init_dpressure = .true. initialize_dpressure = 1 end function initialize_dpressure function dpfun (pbar) use splines real :: pbar, dpfun type (spline), save :: spl if(init_dpressure) call new_spline(nw, psi_bar, pressure, spl) dpfun = dsplint(pbar, spl)/psi_N * beta_0/2. end function dpfun function initialize_beta (init) integer :: init, initialize_beta init_beta = .false. if(init == 1) init_beta = .true. initialize_beta = 1 end function initialize_beta function betafun (pbar) use splines real :: pbar, betafun type (spline), save :: spl if(init_beta) call new_spline(nw, psi_bar, beta, spl) betafun = splint(pbar, spl) end function betafun function initialize_bound (init) integer :: init, initialize_bound init_bound = .false. if(init == 1) init_bound = .true. initialize_bound = 1 end function initialize_bound function bound(theta) use splines real :: theta, bound type (spline), save :: spl integer i if(init_bound) call new_spline(nbbbs, thetab, r_bound, spl) init_bound = .false. bound = splint(theta, spl) end function bound subroutine alloc_module_arrays(np, nw, nh, nws, nhs) integer :: np, nw, nh, nws, nhs allocate(psi_bar(np), fp(np), qsf(np), pressure(np), beta(np)) allocate(dummy(nws), efit_R(nw), efit_Z(nh)) allocate(spsi_bar(nws), sefit_R(nws), sefit_Z(nhs)) allocate(efit_psi(nw, nh), sefit_psi(nws, nhs)) allocate(dpm(nw, nh, 2), dtm(nw, nh, 2)) end subroutine alloc_module_arrays subroutine a_minor(r, z, Z_mag, a) use splines real, dimension(:), intent (in) :: r, z real :: a, Z_mag, r1, r2 integer, parameter :: nz = 5 real, dimension(nz) :: rtmp, ztmp integer i, j, i1, n type (spline) :: spl n = size(r) if(n < nz) then write(*,*) 'nbbbs < nz -- very strange. Stopping.' write(*,*) 'Look in eeq.f90.' stop endif j = 0 do i = nz/2+1,1,-1 j = j + 1 ztmp(j) = z(i) rtmp(j) = r(i) enddo if(r(1) == r(n) .and. z(1) == z(n)) then do i = n-1, n-nz/2, -1 j = j + 1 ztmp(j) = z(i) rtmp(j) = r(i) enddo else do i = n, n-nz/2+1, -1 j = j + 1 ztmp(j) = z(i) rtmp(j) = r(i) enddo endif call new_spline(nz, ztmp, rtmp, spl) r1 = splint(Z_mag, spl) call delete_spline(spl) ! find point near magnetic axis elevation on low field side do i = nz, n if(z(i)-Z_mag > 0.) then i1 = i - 1 exit endif enddo do i = 1, nz rtmp(i) = r(i1 - nz/2 + i - 1) ztmp(i) = z(i1 - nz/2 + i - 1) enddo call new_spline(nz, ztmp, rtmp, spl) r2 = splint(Z_mag, spl) call delete_spline(spl) a = (r2 - r1)/2. end subroutine a_minor function rcenter(rp) real :: rcenter, pi, rp pi = 2.*acos(0.) rcenter = R_mag + 0.5*(rfun(rp, 0., bound(0.))-rfun(rp,pi,bound(pi))) end function rcenter function diameter(rp) real :: diameter, pi, rp pi = 2.*acos(0.) diameter = rfun(rp, 0., bound(0.)) + rfun(rp, pi, bound(pi)) end function diameter function rfun(rp, thet, broot) integer :: i, j, k, imax, jmax, kmax real :: rfun, fa, fb, fbroot, bmult, rootval, thetroot real :: a, b, xerrsec, rp, thet, broot thetroot=thet rootval=rp a=0. b=broot if(broot < 0.) broot = bound(thet) fa =psi(a,thetroot)-rootval fbroot=psi(b,thetroot)-rootval fb=fbroot if(fa*fbroot <= 0.) goto 100 ! need to bracket root. bmult=0.25 jmax=5 kmax=5 do k=1,kmax bmult=bmult*2. imax=10 do j=1,jmax imax=imax*2 do i=1,imax b=broot*(1+bmult*float(i)/float(imax)) fb=psi(b,thetroot)-rootval if(fa*fb <= 0.) goto 100 enddo do i=1,imax b=broot/(1+bmult*float(i)/float(imax)) fb=psi(b,thetroot)-rootval if(fa*fb <= 0.) goto 100 enddo enddo enddo write(*,*) 'yuk' 100 continue xerrsec=1.e-9 ! write(*,*) 'a, b = ',a,b rfun=zbrent(psi, a, b, rootval, thetroot, xerrsec) ! write(*,1000) a, b, fa, fb, rfun, thet 1000 format(1x,11e16.9) end function rfun function zbrent(func, x1, x2, rootval, thetroot, tol) real :: zbrent, func ! real, parameter :: eps = 3.e-8, eps1 = 2.e-5 real, parameter :: eps = 3.e-8, eps1 = 2.e-8 integer, parameter :: itmax = 100 real, intent (in) :: x1, x2, rootval, thetroot, tol real :: a, b, c, fa, fb, fc, d, e, tol1, xm, q, p, r, s integer :: iter a=x1 b=x2 c=0. fa=func(a,thetroot)-rootval fb=func(b,thetroot)-rootval if(fb*fa > 0.) then ! check to see if fa or fb is close to zero already: if(abs(fa) < eps1) then zbrent=a return elseif(abs(fb) < eps1) then zbrent=b return endif write(*,*) a,b,fa,fb write(*,*) 'root must be bracketed for zbrent.' stop endif fc=fb do 11 iter=1,itmax if(fb*fc > 0.) then c=a fc=fa d=b-a e=d endif if(abs(fc) < abs(fb)) then a=b b=c c=a fa=fb fb=fc fc=fa endif tol1=2.*eps*abs(b)+0.5*tol xm=.5*(c-b) if(abs(xm) <= tol1 .or. fb == 0.)then zbrent=b return endif if(abs(e) >= tol1 .and. abs(fa) > abs(fb)) then s=fb/fa if(a == c) then p=2.*xm*s q=1.-s else q=fa/fc r=fb/fc p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.)) q=(q-1.)*(r-1.)*(s-1.) endif if(p > 0.) q=-q p=abs(p) if(2.*p < min(3.*xm*q-abs(tol1*q),abs(e*q))) then e=d d=p/q else d=xm e=d endif else d=xm e=d endif a=b fa=fb if(abs(d) > tol1) then b=b+d else b=b+sign(tol1,xm) endif fb=func(b,thetroot)-rootval 11 continue write(*,*) 'zbrent exceeding maximum iterations.' stop zbrent=b end function zbrent subroutine sort(a, b, c, d) real, dimension(:) :: a, b, c, d real tmp integer :: i, j, jmax jmax = size(a) do j=1,jmax do i=1,jmax-j if(a(i+1) < a(i)) then tmp=a(i); a(i)=a(i+1); a(i+1)=tmp tmp=b(i); b(i)=b(i+1); b(i+1)=tmp tmp=c(i); c(i)=c(i+1); c(i+1)=tmp tmp=d(i); d(i)=d(i+1); d(i+1)=tmp endif enddo enddo end subroutine sort subroutine Hahm_Burrell(irho, a, rho_star0) real, intent(in) :: a character*1 :: char integer :: i, irho real :: gradpsi2, gradpsi_r, gradpsi_z, mag_B, a real :: rho_eq, rho_star0 real :: rp1, rp2, rho1, rho2, drhodpsiq real, dimension(nw) :: gamma, dp, d2p, pres if(irho /= 2) then write(*,*) 'use irho=2 to get correct shearing rate.' endif gamma = 0. do i=2, nw-1 dp(i) = dpfun(psi_bar(i)) enddo do i=3,nw-2 d2p(i) = (dp(i+1)-dp(i-1))/(psi_bar(i+1)-psi_bar(i-1))/psi_N enddo pres=0. do i=3,nw-2 rp1=psi_bar(i+1)*psi_N + psi0 rp2=psi_bar(i-1)*psi_N + psi0 rho1=0.5*diameter(rp1) rho2=0.5*diameter(rp2) drhodpsiq=(rho1-rho2)/(rp1-rp2) pres(i) = pfun(psi_bar(i)) call eqitem(psi_bar(i), 0., dpm(:,:,1), gradpsi_r) call eqitem(psi_bar(i), 0., dpm(:,:,2), gradpsi_z) gradpsi2 = gradpsi_r**2 + gradpsi_z**2 mag_B = sqrt(btori(psi_bar(i))**2 + gradpsi2)/efit_R(i) gamma(i) = (d2p(i)/pres(i)-a*(dp(i)/pres(i))**2) gamma(i) = rho_star0*gradpsi2*gamma(i) & /mag_B*(2.*pres(i)/beta_0)**((1-a)/2.) & *(-pres(i)/(dp(i)/drhodpsiq)) enddo do i=3,nw-2 if(irho == 1) then rho_eq = psi_bar(i) else if(irho == 2) then rho_eq = 0.5 * diameter(psi_bar(i)*psi_N+psi0) else if(irho == 3) then rho_eq = psi_bar(i) endif write(24,1000) i, psi_bar(i), rho_eq, pres(i), gamma(i) ! write(61,*) rho_eq, gamma(i), d2p(i)/pres(i), (dp(i)/pres(i))**2 !, dp(i)/drhodpsiq/pres(i) ! write(63,*) rho_eq, psi_bar(i), pres(i), dp(i), d2p(i) enddo if(irho == 1) write(* ,*) '# irho = 1 produces psi_bar instead of rho_eq' if(irho == 1) write(24,*) '# irho = 1 produces psi_bar instead of rho_eq' 1000 format(i5,11(1x,e16.9)) end subroutine Hahm_Burrell end module eeq