bugs fixed in eccdeff, density and zeffan

src/coreprofiles.f90

@@ -67,7 +67,11 @@ contains
         if(ier.gt.0) print*,ier
         if(abs(dens).lt.1.0e-10_wp_) dens=zero
       end if
-      if(dens.lt.zero) print*,'    DENSITY NEGATIVE',dens
+!      if(dens.lt.zero) print*,'    DENSITY NEGATIVE',dens
+      if(dens.lt.zero) then
+        dens=zero
+        ddens=zero
+      end if  
     end if
   end subroutine density
 
@@ -267,4 +271,4 @@ contains
     if(allocated(cfn)) deallocate(cfn)
   end subroutine unset_prfspl
 
-end module coreprofiles
+end module coreprofiles

src/eccd.f90

@@ -245,6 +245,7 @@ contains
 
       call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp, &
                    epp,neval,ier,liw,lw,last,iw,w)
+             
       if (ier.gt.0) then
         ierr=90
         return
@@ -275,6 +276,7 @@ contains
 
       resj=0.0_wp_
       do i=0,1
+        resji=0.0_wp_
         xx1=amu*(anpl*uleft(i)+ygn-1.0_wp_)
         xx2=amu*(anpl*uright(i)+ygn-1.0_wp_)
         if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
@@ -879,4 +881,4 @@ contains
     f = t**(3+r2)/g**3 * (gp1/(g+1))**r2
   END FUNCTION HSL_f
 
-end module eccd
+end module eccd

src/gray.f

@@ -665,8 +665,8 @@ c
       use beamdata, only : nrayr,nrayth,nstep
       implicit none
 c local variables
-      integer :: nfil,iox,ierr,isgniphi,isgnbphi
-      real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,zeff,
+      integer :: nfil,iox,ierr,isgniphi,isgnbphi,ifreefmt
+      real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,
      . w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta,ssplf
       real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc,
      . rv,zv,fpol
@@ -776,7 +776,7 @@ c     ipsinorm  0 standard EQDSK format, 1 format Portone summer 2004
 c     sspl      spline parameter for psi interpolation
 c
       read(2,*) filenmeqq
-      read(2,*) ipsinorm,sspl
+      read(2,*) ipsinorm,sspl,ifreefmt
 c
 c     factb     factor for magnetic field (only for numerical equil)
 c     scaling adopted: beta=const, qpsi=const, nustar=const
@@ -802,11 +802,12 @@ c
       if(iprof.eq.0) then
         read(2,*) dens0,aln1,aln2
         read(2,*) te0,dte0,alt1,alt2
+        read(2,*) zeffan
       else
         read(2,*) dummy,dummy,dummy
         read(2,*) dummy,dummy,dummy,dummy
+        read(2,*) dummy
       end if
-      read(2,*) zeff
 c
 c     read parameters for analytical simple cylindical equilibrium
 c     if iequil=0,1
@@ -926,7 +927,7 @@ c
         neqdsk=99
         call read_eqdsk(trim(filenmeqq)//'.eqdsk',rv,zv,psin,psia,
      .     psinr,fpol,qpsi,rcen,rmaxis,zmaxis,rbbbs,zbbbs,rlim,zlim,
-     .     ipsinorm,0,ipsinorm,neqdsk) !,idesc,ifreefmt,neqdsk)
+     .     ipsinorm,1,ifreefmt,neqdsk) !,idesc,ifreefmt,neqdsk)
         call change_cocos(psia,fpol,qpsi,icocos,3)
         call eq_scal(psia,fpol,isgniphi,isgnbphi,factb)
         ssplf=0.01_wp_
@@ -982,7 +983,7 @@ c     set simple limiter as two cylindrical walls at rwallm and r00
       if (iprof.eq.1) then
         write(nfil,907) trim(filenmprf)
       else
-        write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeff
+        write(nfil,913) dens0,aln1,aln2,te0,dte0,alt1,alt2,zeffan
       end if    
       write(nfil,911) fghz,p0mw,iox
       write(nfil,902) x00,y00,z00