From 4bb3841049844716990ce430094f444cd83b39c0 Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Fri, 6 Mar 2020 10:48:12 +0000
Subject: [PATCH] Modified test in density extrapolation at boundary

---
 src/coreprofiles.f90 | 16 +++++++++++++---
 src/magsurf_data.f90 | 12 ++++++------
 2 files changed, 19 insertions(+), 9 deletions(-)

diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90
index c23aeb3..d882b92 100644
--- a/src/coreprofiles.f90
+++ b/src/coreprofiles.f90
@@ -225,7 +225,7 @@ contains
 ! local variables
     integer, parameter :: iopt=0, kspl=3
     integer :: n, npest, lwrkf, ier
-    real(wp_) :: xb, xe, fp, xnv, xxp,xxm,delta2
+    real(wp_) :: xb, xe, fp, xnv, xxp,xxm,delta2,ssplne_loc
     real(wp_), dimension(:), allocatable :: wf, wrkf
     integer, dimension(:), allocatable :: iwrkf
     real(wp_), dimension(1) :: dedge,ddedge,d2dedge
@@ -235,6 +235,8 @@ contains
     lwrkf=n*4+npest*16
     allocate(wrkf(lwrkf),iwrkf(npest),wf(n))
 
+    ssplne_loc=ssplne
+
 !   if necessary, reallocate spline arrays
     if(.not.allocated(psrad)) then
       allocate(psrad(n),ct(n,4),cz(n,4))
@@ -263,8 +265,15 @@ contains
     xb=zero
     xe=psin(n)
     wf(:)=one
-    call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne,npest, &
+    call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne_loc,npest, &
        nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
+! if ier=-1 data are re-fitted using sspl=0
+    if(ier==-1) then
+        write(*,*) 'density curfit: ier=-1. Re-fitting with interpolating spline'
+        ssplne_loc=0.0_wp_
+        call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne_loc,npest, &
+           nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
+    end if
 
 !   compute polinomial extrapolation matching the spline boundary up to the
 !   2nd order derivative, extending the profile up to psi=psdbnd where
@@ -284,7 +293,7 @@ contains
     delta2=(ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp
     xnv=psnpp-ddenpp/d2denpp
     if(delta2 < zero) then
-      if(xnv > psnpp) psdbnd=min(psdbnd,xnv)
+!       if(xnv > psnpp) psdbnd=min(psdbnd,xnv)
     else
       xxm=xnv-sqrt(delta2)
       xxp=xnv+sqrt(delta2)
@@ -293,6 +302,7 @@ contains
       else if (xxp > psnpp) then
         psdbnd=min(psdbnd,xxp)
       end if    
+      write(*,*) 'density psdbnd=',psdbnd
     end if
 
     deallocate(iwrkf,wrkf,wf)
diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90
index 81fbeb0..510e9c9 100644
--- a/src/magsurf_data.f90
+++ b/src/magsurf_data.f90
@@ -176,7 +176,7 @@ contains
     else
       call equinum_psi(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
     end if       
-    qq=btaxis/sqrt(ddpsidrr*ddpsidzz)
+    qq=abs(btaxis)/sqrt(ddpsidrr*ddpsidzz)
     ajphiav=-ccj*(ddpsidrr+ddpsidzz)/rmaxis
     
     psicon(1)=0.0_wp_
@@ -412,7 +412,7 @@ contains
     do jp=1,npsi
       call print_fluxav(rpstab(jp),frhotor(rpstab(jp)),bav(jp),bmxpsi(jp),  &
                      bmnpsi(jp),varea(jp),vvol(jp),vcurrp(jp),vajphiav(jp), &
-                     ffc(jp),vratja(jp),vratjb(jp))
+                     ffc(jp),vratja(jp),vratjb(jp),qqv(jp))
     end do 
 
     ninpr=(npsi-1)/10
@@ -558,18 +558,18 @@ contains
 
 
   subroutine print_fluxav(psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, &
-                          ffc,ratja,ratjb)
+                          ffc,ratja,ratjb,qq)
     use const_and_precisions, only : wp_, comp_tiny
     use units, only : uflx
     implicit none
 ! arguments
     real(wp_), intent(in) :: psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, &
-                          ffc,ratja,ratjb
+                          ffc,ratja,ratjb,qq
     
     if (psin < comp_tiny) &
-      write(uflx,*)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol |I_pl| <J_phi> fc ratJa ratJb'
+      write(uflx,*)' #rhop rhot |<B>| |Bmx| |Bmn| Area Vol |I_pl| <J_phi> fc ratJa ratJb qq'
     write(uflx,'(20(1x,e12.5))') psin,rhot,bav,bmx,bmn,area,vol,currp,ajphiav, &
-                          ffc,ratja,ratjb
+                          ffc,ratja,ratjb,qq
   end subroutine print_fluxav
 
 end module magsurf_data