From 75b77b01c1fe4a90ec804c9e4811d80106d06cb1 Mon Sep 17 00:00:00 2001
From: Daniela Farina <daniela.farina@istp.cnr.it>
Date: Tue, 9 Jun 2015 17:09:53 +0000
Subject: [PATCH]  frhotor_eq now is a function of rhop

---
 src/gray.f | 28 ++++++++++++++++------------
 1 file changed, 16 insertions(+), 12 deletions(-)

diff --git a/src/gray.f b/src/gray.f
index cd3f591..6a89916 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -1983,8 +1983,10 @@ c
       subroutine rhotor(nr)
       implicit real*8(a-h,o-z)
       parameter(nnw=501)
-      dimension psinr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),crhot(nnw,4)
+      dimension psinr(nnw),rhopnr(nnw),qpsi(nnw),rhotnr(nnw),cq(nnw,4),
+     *          crhot(nnw,4)
       common/psinr/psinr
+      common/rhopnr/rhopnr
       common/qpsi/qpsi
       common/rhotsx/rhotsx
       common/crhot/crhot
@@ -2010,7 +2012,8 @@ c
 c     spline interpolation of rhotor
 c
       iopt=0
-      call difcs(psinr,rhotnr,nr,iopt,crhot,ier)
+      rhopnr=sqrt(psinr)
+      call difcs(rhopnr,rhotnr,nr,iopt,crhot,ier)
       return
       end
 
@@ -2029,27 +2032,28 @@ c
       return
       end
 
-      function frhotor_eq(psi)
+      function frhotor_eq(rhop)
       implicit real*8(a-h,o-z)
       parameter(nnw=501)
-      dimension psinr(nnw),crhot(nnw,4)
-      common/psinr/psinr
+      dimension rhopnr(nnw),crhot(nnw,4)
+      common/rhopnr/rhopnr
       common/eqnn/nr,nz,npp,nintp
       common/crhot/crhot
 c
-      irt=int((nr-1)*psi+1)
+c      irt=int((nr-1)*psi+1)
 c      if(irt.eq.0) irt=1
 c      if(irt.eq.nr) irt=nr-1
+      call locate(rhopnr,nr,rhop,irt)
       irt=min(max(1,irt),nr-1)
-      dps=psi-psinr(irt)
-      frhotor_eq=spli(crhot,nr,irt,dps)
+      drh=rhop-rhopnr(irt)
+      frhotor_eq=spli(crhot,nr,irt,drh)
       return
       end
 
       function frhotor(psi)
       implicit real*8(a-h,o-z)
       common/iieq/iequil
-      if(iequil.eq.2) frhotor=frhotor_eq(psi)
+      if(iequil.eq.2) frhotor=frhotor_eq(sqrt(psi))
       if(iequil.eq.1) frhotor=frhotor_an(sqrt(psi))
       return
       end
@@ -2609,9 +2613,9 @@ c     ratio_pltor  = Jcd_||/J_phi         Jcd_|| = <Jcd.b>
         qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi))
         qqv(jp)=qq
         
-        dadrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2)
+        dadrhotv(jp)=phitedge*frhotor_eq(height)/fq_eq(height2)
      .   *dadpsi/pi
-        dvdrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2)
+        dvdrhotv(jp)=phitedge*frhotor_eq(height)/fq_eq(height2)
      .   *dvdpsi/pi
 c        
 c      write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp)
@@ -2678,7 +2682,7 @@ c  ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
               rpstab(jp)=1.0d0
               pstab(jp)=1.0d0
         end if
-        rhot_eq=frhotor_eq(pstab(jp))
+        rhot_eq=frhotor_eq(sqrt(pstab(jp)))
         write(56,99) pstab(jp),rhot_eq,rhotqv(jp),
      .    bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), 
      .    vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)