subroutine geo(dpsidrho, bi) ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! ! Initialize qinp, shat, b_mag, rkperp2, omegad, bgrad for general geometry ! from a run of eiktest. use itg_data use mp use gryffin_layouts use fitpack_routines use gryffin_grid, only: ld, r, mrr, md, y0, r0n, r0, meq implicit none integer :: ntgrid real, dimension(:), allocatable :: theta, bmag, gbdrift, gbdrift0, & gds2, gds21, gds22, cvdrift, cvdrift0, gradpar, logb, tgradrho, dum real :: dpsidrho, drhodpsin, rmaj, bi character*40 char integer i, ntg, l, m, n ! ! read in data ! if(proc0) then open(unit=21,file=runname(1:lrunname)//'.geo',status='old') read(21,*) char read(21,*) ntgrid, i, i, drhodpsin, rmaj, shr, bi endif call broadcast(ntgrid) call broadcast(drhodpsin) call broadcast(rmaj) call broadcast(shr) call broadcast(bi) dpsidrho = 1./drhodpsin ntg = 2*ntgrid+1 allocate(theta(ntg), bmag(ntg), gbdrift(ntg), gbdrift0(ntg), & gds2(ntg), gds21(ntg), gds22(ntg), cvdrift(ntg), & cvdrift0(ntg), gradpar(ntg), logb(ntg), tgradrho(ntg)) if(proc0) then read(21,*) char do i=1,ntg read(21,*) gbdrift(i), gradpar(i), tgradrho(i), theta(i) gbdrift(i) = gbdrift(i)*rmaj/4. gradpar(i) = gradpar(i)*rmaj*epsn tgradrho(i) = tgradrho(i)/rmaj/epsn enddo endif call broadcast(gbdrift) call broadcast(gradpar) call broadcast(tgradrho) call broadcast(theta) if(proc0) then read(21,*) char do i=1,ntg read(21,*) cvdrift(i), gds2(i), bmag(i) cvdrift(i) = cvdrift(i)*rmaj/4. logb(i) = log(bmag(i)) enddo endif call broadcast(cvdrift) call broadcast(gds2) call broadcast(bmag) call broadcast(logb) if(proc0) then read(21,*) char do i=1,ntg read(21,*) gds21(i), gds22(i) enddo endif call broadcast(gds21) call broadcast(gds22) if(proc0) then read(21,*) char do i=1,ntg read(21,*) cvdrift0(i), gbdrift0(i) cvdrift0(i) = cvdrift0(i)*rmaj/4. gbdrift0(i) = gbdrift0(i)*rmaj/4. enddo endif call broadcast(cvdrift0) call broadcast(gbdrift0) if(proc0) close(21) allocate (dum(ld)) call inter_d_cspl(ntg,theta,logb,ld,r,dum,bgrad) deallocate (dum) call inter_cspl(ntg,theta,bmag,ld,r,b_mag) call inter_cspl(ntg,theta,gbdrift, ld,r,gbd ) call inter_cspl(ntg,theta,gbdrift0,ld,r,gbd0) call inter_cspl(ntg,theta,cvdrift, ld,r,cvd ) call inter_cspl(ntg,theta,cvdrift0,ld,r,cvd0) call inter_cspl(ntg,theta,gds2, ld,r,gd2 ) call inter_cspl(ntg,theta,gds21,ld,r,gd21) call inter_cspl(ntg,theta,gds22,ld,r,gd22) call inter_cspl(ntg,theta,tgradrho,ld,r,gradrho) gpar = gradpar(1) do l=1,ld bgrad(l)=bgrad(l)*gpar*epsn enddo ! make omegad and rkperp2: !cfpp$ skip R do i=1,nspecies do n = n_low, n_high do m = m_low, m_high do l=1,ld rkx(l,m,n)=(mrr(m,n)/y0)*sqrt(gd22(l))/b_mag(l) ! ! The 1/B above needs to be handled consistently throughout ! the code. ! omega_gb(l,m,n,i)=epsn*rky(m,n)*rho(i)*vt(i) & *(gbd(l)+r0(m,n)*gbd0(l)) omega_kap(l,m,n,i)=epsn*rky(m,n)*rho(i)*vt(i) & *(cvd(l)+r0(m,n)*cvd0(l)) rkperp2(l,m,n)=rky2(m,n) & *(gd2(l)+2.*r0(m,n)*gd21(l)+r0(m,n)**2*gd22(l)) ! ! special cases for k_theta=0, k_r=const. ! if ((m == meq).and.(md /= 1) .and. idx_local(m_lo, meq, n)) then rkx(l,m,n) = -1./y0*r0n(n)*shr/b_mag(l) omega_gb(l,m,n,i)=epsn*rho(i)*vt(i)/y0*r0n(n)*gbd0(l) omega_kap(l,m,n,i)=epsn*rho(i)*vt(i)/y0*r0n(n)*cvd0(l) rkperp2(l,m,n)=(r0n(n)/y0)**2*gd22(l) endif if (md == 1 .and. mrr(m,n) == 0) then rkx(l,m,n)=1./y0 omega_gb(l,m,n,i)=epsn*rho(i)*vt(i)/y0*gbd0(l) omega_kap(l,m,n,i)=epsn*rho(i)*vt(i)/y0*cvd0(l) rkperp2(l,m,n)=1./y0**2 endif ! ! should be krho ! rkx2(l,m,n)=rkx(l,m,n)**2 rkperp2(l,m,n)=rkperp2(l,m,n)/b_mag(l)**2 enddo enddo enddo enddo !cfpp$ noskip R return end subroutine geo