flux_average now in magsurf_data works anal/num , new subr fluxval, some routines moved into equilibrium, minor bugs fixed

Makefile

@@ -44,8 +44,10 @@ dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
 eccd.o: const_and_precisions.o conical.o magsurf_data.o dierckx.o numint.o
 eierf.o: const_and_precisions.o
 gray_params.o: const_and_precisions.o utils.o
-equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o utils.o
+equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o \
-magsurf_data.o: const_and_precisions.o
+               utils.o gray_params.o
+magsurf_data.o: const_and_precisions.o gray_params.o equilibrium.o dierckx.o \
+               simplespline.o utils.o
 math.o: const_and_precisions.o
 minpack.o: const_and_precisions.o
 numint.o: const_and_precisions.o

src/equilibrium.f90

@@ -2,32 +2,35 @@ module equilibrium
   use const_and_precisions, only : wp_
   implicit none
   
-  REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis,sgnbphi
+  real(wp_), save :: btaxis,rmaxis,zmaxis,sgnbphi
-  REAL(wp_), SAVE :: btrcen ! used only for Jcd_ASTRA def.
+  real(wp_), save :: btrcen ! used only for jcd_astra def.
-  REAL(wp_), SAVE :: rcen   ! computed as fpol(a)/btrcen
+  real(wp_), save :: rcen   ! computed as fpol(a)/btrcen
-  REAL(wp_), SAVE :: rmnm,rmxm,zmnm,zmxm
+  real(wp_), save :: rmnm,rmxm,zmnm,zmxm
-  REAL(wp_), SAVE :: zbinf,zbsup
+  real(wp_), save :: zbinf,zbsup
-!  REAL(wp_), SAVE :: rup,zup,rlw,zlw
+!  real(wp_), save :: rup,zup,rlw,zlw
 
-  INTEGER, PARAMETER :: kspl=3,ksplp=kspl+1
+  integer, parameter :: kspl=3,ksplp=kspl+1
-  ! === 2D spline psi(R,z), normalization and derivatives ==========
+
-  INTEGER, SAVE :: nsr, nsz
+! === 2d spline psi(r,z), normalization and derivatives ==========
-  REAL(wp_), SAVE :: psia, psiant, psinop
+  integer, save :: nsr, nsz
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tr,tz
+  real(wp_), save :: psia, psiant, psinop
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq,   cceq01, cceq10, &
+  real(wp_), dimension(:), allocatable, save :: tr,tz
+  real(wp_), dimension(:), allocatable, save :: cceq,   cceq01, cceq10, &
                                                 cceq20, cceq02, cceq11
-  ! === 1D spline Fpol(psi) ========================================
+
-  ! INTEGER, SAVE :: npsiest
+! === 1d spline fpol(psi) ========================================
-  INTEGER, SAVE :: nsf
+! integer, save :: npsiest
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfp, cfp
+  integer, save :: nsf
-  REAL(wp_), SAVE :: fpolas
+  real(wp_), dimension(:), allocatable, save :: tfp, cfp
-  ! === 1D spline rhot(rhop), rhop(rhot), q(psi) ===================
+  real(wp_), save :: fpolas
+
+! === 1d spline rhot(rhop), rhop(rhot), q(psi) ===================
 ! computed on psinr,rhopnr [,rhotnr] arrays
-  INTEGER, SAVE :: nq,nrho
+  integer, save :: nq,nrho
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,rhopr,rhotr
+  real(wp_), dimension(:), allocatable, save :: psinr,rhopr,rhotr
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cq,crhop,crhot
+  real(wp_), dimension(:,:), allocatable, save :: cq,crhop,crhot
-  REAL(wp_), SAVE :: phitedge,aminor
+  real(wp_), save :: phitedge,aminor
-  REAL(wp_), SAVE :: q0,qa,alq
+  real(wp_), save :: q0,qa,alq
 
 contains
 
@@ -146,15 +149,14 @@ contains
     if(in==0) psin=(psin-psiaxis)/psia
   end subroutine read_eqdsk
 
+
+
   subroutine read_equil_an(filenm,rv,zv,fpol,q,unit)
     use utils, only : get_free_unit
     implicit none
 ! arguments
     character(len=*), intent(in) :: filenm
     integer, optional, intent(in) :: unit
-!    integer, intent(in) :: isgnbphi
-!    real(wp_), intent(in) :: factb
-!    real(wp_), intent(out) :: rvac,rax,zax
     real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q
 ! local variables
     integer :: u
@@ -199,7 +201,6 @@ contains
     real(wp_), intent(inout) :: psia
     real(wp_), dimension(:), intent(inout) :: fpol,q
     integer, intent(in) :: cocosin, cocosout
-!    real(wp_), intent(out) :: isign,bsign
     integer, intent(out), optional :: ierr
 ! local variables
     real(wp_) :: isign,bsign
@@ -222,12 +223,9 @@ contains
     end if
 
 ! convert cocosin to cocosout
-    if (phiccw.neqv.phiccwout) then
+
 ! opposite direction of toroidal angle phi in cocosin and cocosout
-!      bsign=-bsign
+    if (phiccw.neqv.phiccwout) fpol=-fpol
-!      isign=-isign
-      fpol=-fpol
-    end if
 ! q has opposite sign for given sign of Bphi*Ip
     if (qpos .neqv. qposout) q=-q
 ! psi and Ip signs don't change accordingly
@@ -301,7 +299,7 @@ contains
     real(wp_), dimension(:), allocatable :: fvpsi,wf,wrk
     integer, dimension(:), allocatable :: iwrk
     integer :: ier,ixploc,info
-!
+
 ! compute array sizes and prepare working space arrays
     nr=size(rv)
     nz=size(zv)
@@ -310,17 +308,17 @@ contains
     nzest=nz+ksplp
     lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest)
     liwrk=nz+nr+nrest+nzest+kspl
-!
+
     npsi=size(psinr)
     npsest=npsi+ksplp
     lwrkf=npsi*ksplp+npsest*(7+3*kspl)
-!
+
     allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest)))
-!
+
 ! spline fitting/interpolation of psin(i,j) and derivatives
-!
+
 ! length in m !!!
-!
+
     rmnm=rv(1)
     rmxm=rv(nr)
     zmnm=zv(1)
@@ -362,9 +360,9 @@ contains
     call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,2,0,cceq20,lw20,ier)
     call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier)
     call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier)
-!
+
 ! spline interpolation of fpol(i)
-!
+
 ! allocate knots and spline coefficients arrays
     if (allocated(tfp)) deallocate(tfp)
     if (allocated(cfp)) deallocate(cfp)
@@ -380,17 +378,17 @@ contains
     sgnbphi=sign(one,fpolas)
 ! free temporary arrays
     deallocate(wrk,iwrk,wf)
-!
+
 ! re-normalize psi after spline computation
-!
+
 ! start with un-corrected psi
-!
+
     psia=psiwbrad
     psinop=0.0_wp_
     psiant=1.0_wp_
-!
+
 ! use provided boundary to set an initial guess for the search of O/X points
-!
+
     nbnd=0
     if(present(rbnd).and.present(zbnd)) then
       nbnd=min(size(rbnd),size(zbnd))
@@ -408,9 +406,9 @@ contains
       rbmin=rv(2)
       rbmax=rv(nr-1)
     end if
-!
+
 ! search for exact location of the magnetic axis
-!
+
     if(present(rax)) then
       rax0=rax
     else
@@ -423,9 +421,9 @@ contains
     end if
     call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info)
     print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp
-!
+
 ! search for X-point   if ixp not = 0
-!
+
     if(present(ixp)) then
       ixploc=ixp
     else
@@ -473,10 +471,10 @@ contains
       print'(a,4f8.4)','no X-point  ',rbinf,zbinf,rbsup,zbsup
     end if
     print*,' '
-!
+
 ! save Bt value on axis (required in flux_average and used in Jcd def)
 ! and vacuum value B0 at ref. radius R0 (used in Jcd_astra def)
-!
+
     call equinum_fpol(0.0_wp_,btaxis)
     btaxis=btaxis/rmaxis
     if(present(r0)) then
@@ -522,9 +520,9 @@ contains
     integer, dimension(liwrk) :: iwrk
     real(wp_), dimension(1) :: rrs,zzs,ffspl
     real(wp_), dimension(lwrk) :: wrk
-!
+
 !   here lengths are measured in meters
-!
+
     if (rpsim.le.rmxm .and. rpsim.ge.rmnm .and. &
         zpsim.le.zmxm .and. zpsim.ge.zmnm) then
 
@@ -608,6 +606,7 @@ contains
   end subroutine equinum_fpol
 
   subroutine bfield(rpsim,zpsim,bphi,br,bz)
+    use gray_params, only : iequil
     implicit none
 ! arguments
     real(wp_), intent(in) :: rpsim,zpsim
@@ -615,12 +614,17 @@ contains
 ! local variables
     real(wp_) :: psin,fpol
 
+    if (iequil < 2) then
+      call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br)
+      if (present(bphi)) bphi=bphi/rpsim
+    else
       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
+    end if
     if (present(br)) br=-br/rpsim
     if (present(bz)) bz= bz/rpsim
   end subroutine bfield
@@ -646,9 +650,8 @@ contains
 
     psinr=psinq
     call difcs(psinr,q,nq,iopt,cq,ier)
-!
+
 ! Toroidal flux phi = 2*pi*Integral q dpsi
-!
     phit(1)=0.0_wp_
     do k=1,nq-1
       dx=psinr(k+1)-psinr(k)
@@ -818,7 +821,7 @@ contains
     real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
 ! local variables
     real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
-!
+
     select case(iflag)
     case(1)
       call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz)
@@ -851,7 +854,6 @@ contains
     real(wp_), dimension(n) :: xvec,fvec,f0
     real(wp_), dimension(lwa) :: wa
     real(wp_), dimension(ldfjac,n) :: fjac
-! common/external functions/variables
 
     xvec(1)=rz
     xvec(2)=zz
@@ -860,8 +862,8 @@ contains
     tol = sqrt(comp_eps)
     call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
     if(info.gt.1)  then
-      print'(a,i2,a,2f8.4)','    info subr points_tgo =',info, &
+      print'(a,i2,a,5f8.4)','    info subr points_tgo =',info, &
-        '  R,z coord.',xvec
+        '  R,z coord.',xvec,rz,zz,psin0
     end if
     rf=xvec(1)
     zf=xvec(2)
@@ -927,12 +929,13 @@ contains
     end if
     if (allocated(psinr)) deallocate(psinr)
     allocate(psinr(nq),rhotn(nq),rhopn(nq))
+    
     dr=one/(nq-1)
     rhotn(1)=zero
     psinr(1)=zero
     res=zero
     fq0=zero
-    do i=2,n
+    do i=2,nq
       rn=(i-1)*dr
       qq=q0+(q1-q0)*rn**qexp
       fq1=rn/qq
@@ -941,11 +944,13 @@ contains
       rhotn(i)=rn
       psinr(i)=res
     end do
+    
     phitedge=btaxis*aminor**2 ! temporary
     psia=res*phitedge
     phitedge=pi*phitedge    ! final
     psinr=psinr/res
     rhopn=sqrt(psinr)
+    
     call set_rhospl(rhopn,rhotn)
   end subroutine set_equian
 
@@ -958,8 +963,7 @@ contains
     real(wp_), intent(out), optional :: psinv,fpolv,dfpv,dpsidr,dpsidz, &
       ddpsidrr,ddpsidzz,ddpsidrz
 ! local variables
-    real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,snt,rhop,rhot
+    real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,snt,rhop,rhot,btaxqq
-!    real(wp_) :: frhopol
 
 !  simple model for equilibrium: large aspect ratio
 !  outside plasma: analytical continuation, not solution Maxwell eqs
@@ -987,12 +991,14 @@ contains
 
     if(rn <= 1.0_wp_) then
       qq=q0+(qa-q0)*rn**alq
-      dpsidrp=btaxis*aminor*rn/qq
+      btaxqq=btaxis/qq
+      dpsidrp=btaxqq*rpm
       dqq=alq*(qa-q0)*rn**(alq-1.0_wp_)
-      d2psidrp=btaxis/qq*(1.0_wp_-rn*dqq/qq)
+      d2psidrp=btaxqq*(1.0_wp_-rn*dqq/qq)
     else
-      dpsidrp=btaxis*rpm/qa
+      btaxqq=btaxis/qa
-      d2psidrp=btaxis/qa
+      dpsidrp=btaxqq*rpm
+      d2psidrp=btaxqq
     end if
 
     if(present(fpolv)) fpolv=btaxis*rmaxis
@@ -1000,10 +1006,80 @@ contains
 
     if(present(dpsidr)) dpsidr=dpsidrp*cst
     if(present(dpsidz)) dpsidz=dpsidrp*snt
-    if(present(ddpsidrr)) ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp
+    if(present(ddpsidrr)) ddpsidrr=btaxqq*snt**2+cst**2*d2psidrp
-    if(present(ddpsidrz)) ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm)
+    if(present(ddpsidrz)) ddpsidrz=cst*snt*(d2psidrp-btaxqq)
-    if(present(ddpsidzz)) ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp
+    if(present(ddpsidzz)) ddpsidzz=btaxqq*cst**2+snt**2*d2psidrp
-
   end subroutine equian
 
+
+
+  subroutine tor_curr(rpsim,zpsim,ajphi)
+    use const_and_precisions, only : wp_,ccj=>mu0inv
+    use gray_params, only : iequil
+    implicit none
+! arguments
+    real(wp_) :: rpsim,zpsim,ajphi
+! local variables
+    real(wp_) :: bzz,dbvcdc13,dbvcdc31
+    real(wp_) :: dpsidr,ddpsidrr,ddpsidzz
+    
+    if(iequil < 2) then
+      call equian(rpsim,zpsim,dpsidr=dpsidr,                     &
+                  ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
+    else 
+      call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,                &
+                       ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
+    end if 
+    bzz= dpsidr/rpsim
+    dbvcdc13=-ddpsidzz/rpsim
+    dbvcdc31= ddpsidrr/rpsim-bzz/rpsim
+    ajphi=ccj*(dbvcdc13-dbvcdc31)
+  end subroutine tor_curr
+
+
+
+  subroutine psi_raxis(psin,r1,r2)
+    use const_and_precisions, only : wp_
+    use gray_params, only : iequil
+    use dierckx, only : profil,sproota
+    implicit none
+! local constants
+    integer, parameter :: mest=4
+! arguments
+    real(wp_) :: psin,r1,r2
+! local variables
+    integer :: iopt,ier,m
+    real(wp_) :: zc,val
+    real(wp_), dimension(mest) :: zeroc
+    real(wp_), dimension(nsr) :: czc
+
+    if (iequil < 2) then
+      val=frhotor(sqrt(psin))
+      r1=rmaxis-val*aminor
+      r2=rmaxis+val*aminor
+    else
+      iopt=1
+      zc=zmaxis
+      call profil(iopt,tr,nsr,tz,nsz,cceq,kspl,kspl,zc,nsr,czc,ier)
+      if(ier.gt.0) print*,'    profil =',ier
+      val=psin*psiant+psinop
+      call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
+      r1=zeroc(1)
+      r2=zeroc(2)
+    end if  
+ end subroutine psi_raxis
+
+
+
+  subroutine tor_curr_psi(psin,ajphi)
+    use const_and_precisions, only : wp_
+    implicit none
+! arguments
+    real(wp_) :: psin,ajphi
+! local variables
+    real(wp_) :: r1,r2
+    call psi_raxis(psin,r1,r2)
+    call tor_curr(r2,zmaxis,ajphi)
+  end subroutine tor_curr_psi
+
 end module equilibrium

src/graycore.f90

@@ -13,6 +13,7 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp,   &
      rtrparam_type, hcdparam_type, set_codepar, igrad, iequil, iprof
   use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff
   use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl
+  use magsurf_data, only : flux_average
   use beamdata, only : init_rtr, dealloc_beam, nstep, dst
   use reflections, only : set_lim
   implicit none
@@ -49,7 +50,7 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp,   &
 
   if(iequil<2) then
     call set_equian(rv(1),zv(1),rv(2), fpol(1)/rv(1), qpsi(1),qpsi(2),qpsi(3))
-    call flux_average_an
+    call flux_average
   else
     call set_eqspl(rv,zv,psin, psia, psinr,fpol, eqp%ssplps,eqp%ssplf, rvac,  &
        rax,zax, rbnd,zbnd, eqp%ixp)
@@ -97,6 +98,7 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp,   &
   call prfile
   call paraminit
   call vectinit
+
   if(igrad==0) then
     call ic_rt(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, &
        w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0)
@@ -104,8 +106,9 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp,   &
     call ic_gb(x0,y0,z0,anx0,any0,anz0,ak0,xgcn,bres, &
        w1,w2,ri1,ri2,phiw,phir,sox,psipol0,chipol0)
   end if
-!  if(ierr==0) return
+
 ! beam/ray propagation
+
 ! st0=0.0_wp_
 ! if(index_rt.gt.1) st0=strfl11
   do i=1,nstep
@@ -117,6 +120,7 @@ subroutine gray(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp,   &
     call after_onestep(p0,sox,ak0,xgcn,bres,i,istop,ierr)
     if(istop.eq.1) exit
   end do
+
 ! postprocessing
   call after_gray_integration(pec,icd,dpdv,jcd)
 

src/magsurf_data.f90

@@ -2,23 +2,23 @@ module magsurf_data
   use const_and_precisions, only : wp_
   implicit none
   
-  INTEGER, SAVE :: npsi, npoints !# sup mag, # punti per sup
+  integer, save :: npsi, npoints !# sup mag, # punti per sup
-  INTEGER, SAVE :: njpt, nlmt
+  integer, save :: njpt, nlmt
 
-  REAL(wp_), SAVE :: rarea
+  real(wp_), save :: rarea
 
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psicon,pstab,rhot_eq, &
+  real(wp_), dimension(:), allocatable, save :: psicon,pstab,rhot_eq, &
    rhotqv,bav,varea,vcurrp,vajphiav,qqv,ffc,vratja,vratjb
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rpstab
+  real(wp_), dimension(:), allocatable, save :: rpstab
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: vvol,rri,rbav,bmxpsi,bmnpsi
+  real(wp_), dimension(:), allocatable, save :: vvol,rri,rbav,bmxpsi,bmnpsi
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tjp,tlm,ch,ch01
+  real(wp_), dimension(:), allocatable, save :: tjp,tlm,ch,ch01
 
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: rcon,zcon
+  real(wp_), dimension(:,:), allocatable, save :: rcon,zcon
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cdadrhot,cdvdrhot
+  real(wp_), dimension(:,:), allocatable, save :: cdadrhot,cdvdrhot
   
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cvol,crri,crbav,cbmx,cbmn,carea,cfc
+  real(wp_), dimension(:,:), allocatable, save :: cvol,crri,crbav,cbmx,cbmn,carea,cfc
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: crhotq
+  real(wp_), dimension(:,:), allocatable, save :: crhotq
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cratja,cratjb,cratjpl
+  real(wp_), dimension(:,:), allocatable, save :: cratja,cratjb,cratjpl
 
 
 contains
@@ -103,4 +103,383 @@ contains
     if(allocated(cratjb))  deallocate(cratjb)
   end subroutine dealloc_surfvec
 
+
+
+  subroutine flux_average
+    use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv
+    use gray_params, only : iequil
+    use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, &
+                            equian,equinum_psi,bfield,frhotor,fq,tor_curr
+    use simplespline, only : difcs
+    use dierckx, only : regrid,coeff_parder
+    use utils, only : get_free_unit
+    implicit none
+! local constants
+    integer, parameter :: nnintp=101,ncnt=100,nlam=101,ksp=3,     &
+             njest=nnintp+ksp+1,nlest=nlam+ksp+1,     &
+             lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, &
+             kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam
+! local variables
+    integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr,u56,unit
+    integer, dimension(kwrk) :: iwrk
+    real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator,     &
+                 ratio_cdbtor,ratio_pltor,fc,height,height2,r2iav,currp,     &
+                 area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, &
+                 dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0,     &
+                 shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi,     &
+                 bphi,brr,bzz,riav,fp,psinjp,rhopjp,rhotjp,qq,rup,rlw,zup,zlw
+    real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl
+    real(wp_), dimension(2*ncnt) :: dlpv
+    real(wp_), dimension(2*ncnt+1) :: bv,bpv
+    real(wp_), dimension(nlam) :: alam,weights
+    real(wp_), dimension(nnintp,nlam) :: fhlam
+    real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam
+    real(wp_), dimension(lwrk) :: wrk
+    real(wp_), dimension(:), allocatable :: rctemp,zctemp
+! common/external functions/variables
+    real(wp_) :: fpolv,ddpsidrr,ddpsidzz
+
+    u56=get_free_unit()
+    open(file='f56.txt',unit=u56)
+
+    npsi=nnintp
+    ninpr=(npsi-1)/10
+    npoints = 2*ncnt+1
+
+    call alloc_surfvec(ierr)
+    if(allocated(tjp)) deallocate(tjp)
+    if(allocated(tlm)) deallocate(tlm)
+    if(allocated(ch)) deallocate(ch)
+    allocate(tjp(njest),tlm(nlest),ch((njest-4)*(nlest-4)),   &
+             rctemp(npoints),zctemp(npoints),stat=ierr)
+    if (ierr.ne.0) return
+
+! computation of flux surface averaged quantities
+
+    write(71,*)' #i psin R  z'
+
+    dlam=1.0_wp_/dble(nlam-1)
+    do l=1,nlam-1
+      alam(l)=dble(l-1)*dlam
+      fhlam(1,l)=sqrt(1.0_wp_-alam(l))
+      ffhlam(l)=fhlam(1,l)
+      dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l))
+      weights(l)=1.0_wp_
+    end do 
+    weights(1)=0.5_wp_
+    weights(nlam)=0.5_wp_
+    alam(nlam)=1.0_wp_
+    fhlam(1,nlam)=0.0_wp_
+    ffhlam(nlam)=0.0_wp_
+    dffhlam(nlam)=-99999.0_wp_
+
+    jp=1
+    anorm=2.0_wp_*pi*rmaxis/abs(btaxis)
+    dvdpsi=2.0_wp_*pi*anorm
+    dadpsi=2.0_wp_*pi/abs(btaxis)
+    b2av=btaxis**2
+    ratio_cdator=abs(btaxis/btrcen)
+    ratio_cdbtor=1.0_wp_
+    ratio_pltor=1.0_wp_
+    fc=1.0_wp_
+    if(iequil < 2) then
+      call equian(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
+    else
+      call equinum_psi(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
+    end if       
+    qq=btaxis/sqrt(ddpsidrr*ddpsidzz)
+    ajphiav=-ccj*(ddpsidrr+ddpsidzz)/rmaxis
+    
+    psicon(1)=0.0_wp_
+    rcon(1,:)=rmaxis
+    zcon(1,:)=zmaxis
+    pstab(1)=0.0_wp_
+    rpstab(1)=0.0_wp_
+    vcurrp(1)=0.0_wp_
+    vajphiav(1)=ajphiav
+    bmxpsi(1)=abs(btaxis)
+    bmnpsi(1)=abs(btaxis)
+    bav(1)=abs(btaxis)
+    rbav(1)=1.0_wp_
+    rri(1)=rmaxis
+    varea(1)=0.0_wp_
+    vvol(1)=0.0_wp_
+    vratjpl(1)=ratio_pltor
+    vratja(1)=ratio_cdator
+    vratjb(1)=ratio_cdbtor
+    ffc(1)=fc
+    qqv(1)=qq
+    dadrhotv(1)=0.0_wp_
+    dvdrhotv(1)=0.0_wp_
+
+    rup=rmaxis
+    rlw=rmaxis
+    zup=zmaxis+(zbsup-zmaxis)/10.0_wp_
+    zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_
+
+    do jp=2,npsi
+      height=dble(jp-1)/dble(npsi-1)
+      if(jp.eq.npsi) height=0.9999_wp_
+      ipr=0
+      jpr=mod(jp,ninpr)
+      if(jpr.eq.1) ipr=1
+      rhopjp=height
+      psinjp=height*height
+      rhotjp=frhotor(rhopjp)
+      psicon(jp)=height
+
+      if(iequil<2) then
+        call contours_psi_an(psinjp,rctemp,zctemp,ipr) 
+      else
+        call contours_psi(psinjp,rup,zup,rlw,zlw,rctemp,zctemp,ipr)
+      end if
+      rcon(jp,:) = rctemp
+      zcon(jp,:) = zctemp
+
+      r2iav=0.0_wp_
+      anorm=0.0_wp_
+      dadpsi=0.0_wp_
+      currp=0.0_wp_
+      b2av=0.0_wp_
+      area=0.0_wp_
+      volume=0.0_wp_
+      ajphiav=0.0_wp_
+      bbav=0.0_wp_
+      bmmx=-1.0e+30_wp_
+      bmmn=1.0e+30_wp_
+
+      call tor_curr(rctemp(1),zctemp(1),ajphi0)
+      call bfield(rctemp(1),zctemp(1),bphi,br=brr,bz=bzz)
+      fpolv=bphi*rctemp(1)
+      btot0=sqrt(bphi**2+brr**2+bzz**2)
+      bpoloid0=sqrt(brr**2+bzz**2)
+      bv(1)=btot0
+      bpv(1)=bpoloid0
+      rpsim0=rctemp(1)
+
+      do inc=1,npoints-1
+        inc1=inc+1
+        dla=sqrt((rctemp(inc)-rmaxis)**2+(zctemp(inc)-zmaxis)**2)
+        dlb=sqrt((rctemp(inc1)-rmaxis)**2+(zctemp(inc1)-zmaxis)**2)
+        dlp=sqrt((rctemp(inc1)-rctemp(inc))**2+(zctemp(inc1)-zctemp(inc))**2)
+        drc=(rctemp(inc1)-rctemp(inc))
+
+! compute length, area and volume defined by psi=psinjp=height^2
+        ph=0.5_wp_*(dla+dlb+dlp)
+        area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp)
+        area=area+sqrt(area2)
+        rzp=rctemp(inc1)*zctemp(inc1)
+        rz=rctemp(inc)*zctemp(inc)
+        volume=pi*(rzp+rz)*drc+volume
+
+! compute line integrals on the contour psi=psinjp=height^2
+        rpsim=rctemp(inc1)
+        zpsim=zctemp(inc1)
+        call bfield(rpsim,zpsim,br=brr,bz=bzz)
+        call tor_curr(rpsim,zpsim,ajphi)
+        bphi=fpolv/rpsim
+        btot=sqrt(bphi**2+brr**2+bzz**2)
+        bpoloid=sqrt(brr**2+bzz**2)
+        dlpv(inc)=dlp
+        bv(inc1)=btot
+        bpv(inc1)=bpoloid
+        
+        dlph=0.5_wp_*dlp
+        anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0)
+        dadpsi=dadpsi+dlph*(1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0))
+        currp=currp+dlph*(bpoloid+bpoloid0)
+        b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid)
+        bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0)
+        r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2)+1.0_wp_/(bpoloid0*rpsim0**2))
+        ajphiav=ajphiav+dlph*(ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim))
+        
+        ajphi0=ajphi
+        rpsim0=rpsim
+        bpoloid0=bpoloid
+        btot0=btot      
+
+! computation maximum/minimum B values on given flux surface
+        if(btot.le.bmmn) bmmn=btot
+        if(btot.ge.bmmx) bmmx=btot
+      end do     
+
+! bav=<B> [T] ,  b2av=<B^2> [T^2] , rbav=<B>/b_min
+! anorm = int d l_p/B_p = dV/dpsi/(2pi)
+! r2iav=<1/R^2> [m^-2] ,
+! riav=<1/R> [m^-1] = dA/dpsi/(dV/dpsi/(2pi)),
+! rri = <B>/(|R B_tor|<1/R^2>) , used to compute I_tor  [m^-1]
+! currp = plasma current within psi=const
+
+      bbav=bbav/anorm
+      r2iav=r2iav/anorm
+      dvdpsi=2.0_wp_*pi*anorm
+      riav=dadpsi/anorm
+      b2av=b2av/anorm
+      vcurrp(jp)=ccj*currp
+      vajphiav(jp)=ajphiav/dadpsi
+
+! area == varea, volume == vvol
+! flux surface minor radius == (area/pi)^1/2
+! ratio_cdator = Jcd_astra/J_phi    Jcd_astra = <Jcd.B>/B0
+! ratio_cdbtor = Jcd_jintrac/J_phi    Jcd_jintrac  = <Jcd.B>/<B>
+! ratio_pltor  = Jcd_||/J_phi       Jcd_|| = <Jcd.b>
+      pstab(jp)=psinjp
+      rpstab(jp)=rhopjp
+      vvol(jp)=abs(volume)
+      varea(jp)=area
+      bav(jp)=bbav
+      rbav(jp)=bbav/bmmn
+      bmxpsi(jp)=bmmx
+      bmnpsi(jp)=bmmn
+      rri(jp)=bav(jp)/abs(fpolv*r2iav)
+      ratio_cdator=abs(b2av*riav/(fpolv*r2iav*btrcen))
+      ratio_cdbtor=abs(b2av*riav/(fpolv*r2iav*bbav))
+      ratio_pltor=abs(bbav*riav/(fpolv*r2iav))
+      vratjpl(jp)=ratio_pltor
+      vratja(jp)=ratio_cdator
+      vratjb(jp)=ratio_cdbtor
+      qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi))
+      qqv(jp)=qq
+      dadrhotv(jp)=phitedge*frhotor(rhopjp)/fq(psinjp)*dadpsi/pi
+      dvdrhotv(jp)=phitedge*frhotor(rhopjp)/fq(psinjp)*dvdpsi/pi
+      
+! computation of fraction of circulating/trapped fraction fc, ft 
+! and of function H(lambda,rhop)
+! ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
+      fc=0.0_wp_
+      shlam=0.0_wp_
+      do l=nlam,1,-1
+        lam=alam(l)
+        srl=0.0_wp_
+        rl2=1.0_wp_-lam*bv(1)/bmmx
+        rl0=0.0_wp_
+        if(rl2.gt.0) rl0=sqrt(rl2)
+        do inc=1,npoints-1
+          rl2=1.0_wp_-lam*bv(inc+1)/bmmx
+          rl=0.0_wp_
+        if(rl2.gt.0) rl=sqrt(rl2)
+        srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc))
+        rl0=rl
+        end do
+        srl=srl/anorm
+        dhlam=0.5_wp_/srl
+        fc=fc+lam/srl*weights(l)
+        if(l.eq.nlam) then
+        fhlam(jp,l)=0.0_wp_
+        ffhlam(nlam*(jp-1)+l)=0.0_wp_
+        dffhlam(nlam*(jp-1)+l)=-dhlam
+        dhlam0=dhlam
+        else
+        shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam
+        fhlam(jp,l)=shlam
+        dffhlam(nlam*(jp-1)+l)=-dhlam
+        dhlam0=dhlam
+        end if
+      end do
+      fc=0.75_wp_*b2av/bmmx**2*fc*dlam
+      ffc(jp)=fc
+      
+      ccfh=bmmn/bmmx/fc
+      do l=1,nlam
+        ffhlam(nlam*(jp-1)+l)=ccfh*fhlam(jp,l)
+        dffhlam(nlam*(jp-1)+l)=ccfh*dffhlam(nlam*(jp-1)+l)
+      end do
+    end do
+
+    write(u56,*)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol |I_pl| <J_phi> fc ratJa ratJb'
+
+    do jp=1,npsi
+      if(jp.eq.npsi) then
+        rpstab(jp)=1.0_wp_
+        pstab(jp)=1.0_wp_
+      end if
+      rhotjp=frhotor(rpstab(jp))
+      write(u56,99) rpstab(jp),rhotjp,bav(jp),bmxpsi(jp),bmnpsi(jp),  &
+                   varea(jp),vvol(jp),vcurrp(jp),vajphiav(jp),ffc(jp), &
+                   vratja(jp),vratjb(jp)
+    end do 
+
+    close(u56)
+    
+! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs
+! used for computations of dP/dV and J_cd
+    iopt=0
+    call difcs(rpstab,vvol,npsi,iopt,cvol,ier)
+    iopt=0
+    call difcs(rpstab,rbav,npsi,iopt,crbav,ier)
+    iopt=0
+    call difcs(rpstab,rri,npsi,iopt,crri,ier)
+    iopt=0
+    call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier)
+    iopt=0
+    call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier)
+    iopt=0
+    call difcs(rpstab,vratja,npsi,iopt,cratja,ier)
+    iopt=0
+    call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier)
+    iopt=0
+    call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier)
+    iopt=0
+    call difcs(rpstab,varea,npsi,iopt,carea,ier)
+    iopt=0
+    call difcs(rpstab,ffc,npsi,iopt,cfc,ier)
+    iopt=0
+    call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier)
+    iopt=0
+    call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier)
+!    iopt=0
+!    call difcs(rpstab,qqv,npsi,iopt,cqq,ier)
+
+! spline interpolation of H(lambda,rhop) and dH/dlambda
+    iopt=0
+    s=0.0_wp_
+    call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam,zero,one,zero,one,  &
+                ksp,ksp,s,njest,nlest,njp,tjp,nlm,tlm,ch,fp,  &
+                wrk,lwrk,iwrk,kwrk,ier) 
+    njpt=njp
+    nlmt=nlm
+      
+ 99 format(20(1x,e12.5))
+ 
+  end subroutine flux_average
+
+
+
+  subroutine fluxval(rhop,area,vol,dervol,dadrhot,dvdrhot,    &
+                     rri,rbav,bmn,bmx,fc,ratja,ratjb,ratjpl)
+    use const_and_precisions, only : wp_
+    use utils, only : locate
+    use simplespline, only :spli,splid
+    implicit none
+! arguments
+    real(wp_), intent(in) :: rhop
+    real(wp_), intent(out), optional :: vol,area,rri,rbav,dervol,bmn,bmx,fc,   &
+               ratja,ratjb,ratjpl,dadrhot,dvdrhot
+! local variables
+    integer :: ip
+    real(wp_) :: drh
+
+    call locate(rpstab,npsi,rhop,ip)
+    ip=min(max(1,ip),npsi-1)
+    drh=rhop-rpstab(ip)
+
+    if (present(area)) area=spli(carea,npsi,ip,drh)
+    if (present(vol)) vol=spli(cvol,npsi,ip,drh)
+
+    if (present(dervol)) dervol=splid(cvol,npsi,ip,drh)
+    if (present(dadrhot)) dadrhot=spli(cdadrhot,npsi,ip,drh)
+    if (present(dvdrhot)) dvdrhot=spli(cdvdrhot,npsi,ip,drh)
+
+    if (present(rri)) rri=spli(crri,npsi,ip,drh)
+    if (present(rbav)) rbav=spli(crbav,npsi,ip,drh)
+    if (present(bmn)) bmn=spli(cbmn,npsi,ip,drh)
+    if (present(bmx)) bmx=spli(cbmx,npsi,ip,drh)
+    if (present(fc)) fc=spli(cfc,npsi,ip,drh)
+
+    if (present(ratja)) ratja=spli(cratja,npsi,ip,drh)
+    if (present(ratjb)) ratjb=spli(cratjb,npsi,ip,drh)
+    if (present(ratjpl)) ratjpl=spli(cratjpl,npsi,ip,drh)
+
+  end subroutine fluxval
+
 end module magsurf_data