diff --git a/src/gray.f b/src/gray.f index 14bb662..e555630 100644 --- a/src/gray.f +++ b/src/gray.f @@ -1851,7 +1851,7 @@ c c subroutine bfield_res use const_and_precisions, only : wp_ - use interp_eqprof, only : rv,zv,nr,nz,btotal + use interp_eqprof, only : rv,zv,nr,nz,bfield implicit none c local constants integer, parameter :: icmx=2002 @@ -1860,6 +1860,7 @@ c local variables integer, dimension(10) :: ncpts real(wp_) :: btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb real(wp_), dimension(icmx) :: rrcb,zzcb + real(wp_), dimension(nr,nz) :: btotal c common/external functions/variables real(wp_) :: bres c @@ -1888,7 +1889,8 @@ c do n=1,5 bbb=bres/dble(n) if (bbb.ge.btmn.and.bbb.le.btmx) then - call cniteq(bbb,nconts,ncpts,nctot,rrcb,zzcb,1) + call cniteq(rv,zv,btotal,nr,nz,bbb, + . nconts,ncpts,nctot,rrcb,zzcb) do inc=1,nctot write(70,113) inc,bbb,rrcb(inc),zzcb(inc) end do @@ -2221,57 +2223,42 @@ c c c c - subroutine cniteq(h,ncon,npts,icount,rcon,zcon,ichoi) + subroutine cniteq(rqgrid,zqgrid,btotal,nreq,nzeq,h, + . ncon,npts,icount,rcon,zcon) c v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. c (based on an older code) c use const_and_precisions, only : wp_ - use interp_eqprof, only : dr,dz,nr,nz,rqgrid=>rv,zqgrid=>zv,psin, - . btotal implicit none c local constants integer, parameter :: icmx=2002 c arguments - integer :: ncon,icount,ichoi + integer :: nreq,nzeq + real(wp_) :: rqgrid(nreq),zqgrid(nzeq),btotal(nreq,nzeq) + integer :: ncon,icount integer, dimension(10) :: npts real(wp_) :: h real(wp_), dimension(icmx) :: rcon,zcon c local variables - integer :: i,j,k,l,ico,nrqmax,nreq,nzeq,iclast,mpl,ix,jx, + integer :: i,j,k,l,ico,nrqmax,iclast,mpl,ix,jx, . mxr,n1,jm,jfor,lda,ldb,jabs,jnb,kx,ikx,itm,inext,in integer, dimension(3,2) :: ja integer, dimension(1000) :: lx real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y - real(wp_), dimension(nr*nz) :: a + real(wp_), dimension(nreq*nzeq) :: a c data px/0.5d0/ c - if(ichoi.eq.0) then - do j=1,nz - do i=1,nr - a(nr*(j-1)+i)=psin(i,j) - enddo - enddo - endif -c - if(ichoi.eq.1) then - do j=1,nz - do i=1,nr - a(nr*(j-1)+i)=btotal(i,j) - enddo - enddo - endif + a=reshape(btotal,(/nreq*nzeq/)) c do ico=1,icmx rcon(ico)=0.0_wp_ zcon(ico)=0.0_wp_ enddo c - nrqmax=nr - nreq=nr - nzeq=nz - drgrd=dr - dzgrd=dz + nrqmax=nreq + drgrd=rqgrid(2)-rqgrid(1) + dzgrd=zqgrid(2)-zqgrid(1) c ncon = 0 do i=1,10 @@ -2492,7 +2479,7 @@ c use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use magsurf_data use interp_eqprof, only : btrcen,btaxis,rmaxis,zmaxis,phitedge, - . equinum_fpol + . equinum_fpol,bfield use simplespline, only : difcs use dierckx, only : regrid,coeff_parder implicit none @@ -2619,7 +2606,7 @@ c bmmn=1.0e+30_wp_ call tor_curr(rctemp(1),zctemp(1),ajphi0) - call bpol(rctemp(1),zctemp(1),brr,bzz) + call bfield(rctemp(1),zctemp(1),br=brr,bz=bzz) call equinum_fpol(height2,fpolv) bphi=fpolv/rctemp(1) btot0=sqrt(bphi**2+brr**2+bzz**2) @@ -2649,7 +2636,7 @@ c compute line integral on the contour psi=height^2 rpsim=rctemp(inc1) zpsim=zctemp(inc1) - call bpol(rpsim,zpsim,brr,bzz) + call bfield(rpsim,zpsim,br=brr,bz=bzz) call tor_curr(rpsim,zpsim,ajphi) ! call equinum_fpol(height2,fpolv) bphi=fpolv/rpsim @@ -2918,6 +2905,8 @@ c local variables real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam real(wp_), dimension(lwrk) :: wrk real(wp_), dimension(:), allocatable :: rctemp,zctemp + real(wp_) :: psinv !unused. can be removed when equian is modified + !to accept optional arguments c common/external functions/variables real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv real(wp_) :: frhotor_an @@ -3029,7 +3018,7 @@ c bmmx=-1.0e+30_wp_ bmmn=1.0e+30_wp_ - call equian(rctemp(1),zctemp(1)) + call equian(rctemp(1),zctemp(1),psinv) dbvcdc13=-ddpsidzz/rctemp(1) dbvcdc31= ddpsidrr/rctemp(1)-bzz/rctemp(1) ajphi=ccj*(dbvcdc13-dbvcdc31) @@ -3063,7 +3052,7 @@ c compute line integral on the contour psi=height^2 rpsim=rctemp(inc1) zpsim=zctemp(inc1) - call equian(rpsim,zpsim) + call equian(rpsim,zpsim,psinv) brr=-dpsidz/rpsim bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim @@ -4158,7 +4147,7 @@ c rrm=1.0e-2_wp_*rr c if(iequil.eq.1) then - call equian(rrm,zzm) + call equian(rrm,zzm,psinv) end if c if(iequil.eq.2) then @@ -4266,7 +4255,7 @@ c c c c - subroutine equian(rrm,zzm) + subroutine equian(rrm,zzm,psinv) use const_and_precisions, only : wp_ use graydata_par, only : sgnbphi,sgniphi use graydata_anequil, only : b0,rr0m,zr0m,rpam,q0,qa,alq @@ -4274,14 +4263,14 @@ c implicit none c arguments real(wp_) :: rrm,zzm + real(wp_) :: psinv c local variables real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,sgn,snt,rhop,rhot c common/external functions/variables - real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv, + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv, . dfpv,xg,xgcn,dxgdpsi,dens,ddens real(wp_) :: frhopol c - common/psival/psinv common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/fpol/fpolv,dfpv @@ -4338,56 +4327,6 @@ c end c c -c - subroutine bfield(rpsim,zpsim,bphi,brr,bzz) - use const_and_precisions, only : wp_ - implicit none -c arguments - real(wp_) :: rpsim,zpsim,bphi,brr,bzz -c - call btor(rpsim,zpsim,bphi) - call bpol(rpsim,zpsim,brr,bzz) -c - return - end -c -c -c - subroutine btor(rpsim,zpsim,bphi) - use const_and_precisions, only : wp_ - use interp_eqprof, only : equinum_fpol, equinum_psi - implicit none -c arguments - real(wp_) :: rpsim,zpsim,bphi -c local variables - real(wp_) :: psinv,fpolv -c - call equinum_psi(rpsim,zpsim,psinv) - call equinum_fpol(psinv,fpolv) - bphi=fpolv/rpsim -c - return - end -c -c -c - subroutine bpol(rpsim,zpsim,brr,bzz) - use const_and_precisions, only : wp_ - use interp_eqprof, only : equinum_psi - implicit none -c arguments - real(wp_) :: rpsim,zpsim,brr,bzz -c local variables - real(wp_) :: dpsidr,dpsidz -c - call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,dpsidz=dpsidz) - brr=-dpsidz/rpsim - bzz= dpsidr/rpsim -c - return - end -c -c c subroutine tor_curr_psi(h,ajphi) use const_and_precisions, only : wp_ @@ -6015,13 +5954,11 @@ c implicit none c arguments real(wp_) :: rrm,zzm -c common/external functions/variables +c local variables real(wp_) :: psinv -c - common/psival/psinv c if(iequil.eq.1) then - call equian(rrm,zzm) + call equian(rrm,zzm,psinv) else call equinum_psi(rrm,zzm,psinv) end if diff --git a/src/interp_eqprof.f90 b/src/interp_eqprof.f90 index 1da17e1..86d302f 100644 --- a/src/interp_eqprof.f90 +++ b/src/interp_eqprof.f90 @@ -41,10 +41,6 @@ module interp_eqprof ! computed on psinr,rhopnr [,rhotnr] arrays REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: crhot,cq - ! !!! 2D B(R,z) array. Computed in bfield_res, - ! !!! used by cniteq to plot resonant field contour lines. - REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: btotal - contains subroutine alloc_equilvec(ier) @@ -66,7 +62,7 @@ contains call dealloc_equilvec allocate(rv(nr),zv(nz),tr(nrest),tz(nzest),tfp(nrest),cfp(nrest), & - btotal(nr,nz),cceq(nr*nz),cceq01(lw01),cceq10(lw10),cceq02(lw02), & + cceq(nr*nz),cceq01(lw01),cceq10(lw10),cceq02(lw02), & cceq20(lw20),cceq11(lw11),psin(nr,nz),psinr(nr),rhopnr(nr),fpol(nr), & qpsi(nr), stat=ier) if (ier/=0) call dealloc_equilvec @@ -80,7 +76,6 @@ contains if(allocated(tz)) deallocate(tz) if(allocated(tfp)) deallocate(tfp) if(allocated(cfp)) deallocate(cfp) - if(allocated(btotal)) deallocate(btotal) if(allocated(cceq)) deallocate(cceq) if(allocated(cceq01)) deallocate(cceq01) if(allocated(cceq10)) deallocate(cceq10) @@ -243,5 +238,49 @@ contains if (present(dfpv)) dfpv=0._wp_ end if end subroutine equinum_fpol + +! subroutine btor(rpsim,zpsim,bphi) +! implicit none +!! arguments +! real(wp_), intent(in) :: rpsim,zpsim +! real(wp_), intent(out) :: bphi +!! local variables +! real(wp_) :: psinv,fpolv +! +! call equinum_psi(rpsim,zpsim,psinv) +! call equinum_fpol(psinv,fpolv) +! bphi=fpolv/rpsim +! end subroutine btor + +! subroutine bpol(rpsim,zpsim,brr,bzz) +! implicit none +!! arguments +! real(wp_), intent(in) :: rpsim,zpsim +! real(wp_), intent(out) :: brr,bzz +!! local variables +! real(wp_) :: dpsidr,dpsidz +! +! call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,dpsidz=dpsidz) +! brr=-dpsidz/rpsim +! bzz= dpsidr/rpsim +! end subroutine bpol + + subroutine bfield(rpsim,zpsim,bphi,br,bz) + implicit none +! arguments + real(wp_), intent(in) :: rpsim,zpsim + real(wp_), intent(out), optional :: bphi,br,bz +! local variables + real(wp_) :: psin,fpol + + call equinum_psi(rpsim,zpsim,psinv=bphi,dpsidr=bz,dpsidz=br) + if (present(bphi)) then + psin=bphi + call equinum_fpol(psin,fpol) + bphi=fpol/rpsim + end if + if (present(br)) br=-br/rpsim + if (present(bz)) bz= bz/rpsim + end subroutine bfield end module interp_eqprof