updated bfield routine, replacing btor/bpol; btotal passed explicitly to cniteq; equian partially modified (psinv passed as argument)

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
+      a=reshape(btotal,(/nreq*nzeq/))
-         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
 c
       do ico=1,icmx
         rcon(ico)=0.0_wp_
         zcon(ico)=0.0_wp_
       enddo
 c
-      nrqmax=nr
+      nrqmax=nreq
-      nreq=nr
+      drgrd=rqgrid(2)-rqgrid(1)
-      nzeq=nz
+      dzgrd=zqgrid(2)-zqgrid(1)
-      drgrd=dr
-      dzgrd=dz
 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

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)
@@ -244,4 +239,48 @@ contains
     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