cleaned arguments list in diel_tens_fr

src/dispersion.f90

@@ -124,7 +124,7 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
   if (fast.eq.1) then
     call diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
   else
-    call diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
+    call diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
   end if
 !
   do
@@ -250,7 +250,7 @@ end subroutine warmdisp
 !
 !
 !
-subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
+subroutine diel_tens_fr(xg,yg,mu,npl,e330,epsl,lrm,fast)
 ! Fully relativistic case computation of dielectric tensor elements
 ! up to third order in Larmor radius for hermitian part
 !
@@ -258,11 +258,12 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
   implicit none
 ! arguments
   integer :: lrm,fast
-  real(wp_) :: xg,yg,mu,npl,cr,ci
+  real(wp_) :: xg,yg,mu,npl
   complex(wp_) :: e330
   complex(wp_), dimension(3,3,lrm) :: epsl
 ! local variables
   integer :: i,j,l,is,k,lm
+  real(wp_) :: cr,ci
   real(wp_) :: npl2,dnl,w,asl,bsl,cmxw,fal
   real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr
   real(wp_), dimension(lrm,0:2,lrm) :: ri

src/gray.f

@@ -658,12 +658,14 @@ c
       use graydata_par
       use graydata_anequil
       use coreprofiles, only : psdbnd,read_profiles,tene_scal,set_prfspl
-      use interp_eqprof, only : rmxm,zbmin,zbmax,btrcen,rcen
+      use interp_eqprof, only : rmxm,zbmin,zbmax,btrcen,rcen,psia,
+     . rmaxis,zmaxis,fpolas,psinr,qpsi,rv,zv,psin,
+     . fpol,rbbbs,zbbbs,read_eqdsk,change_cocos,eq_scal
       use reflections, only : rlim,zlim,nlim,alloc_lim
       use beamdata, only : nrayr,nrayth,nstep
       implicit none
 c local variables
-      integer :: nfil,iox,ierr,leqq
+      integer :: nfil,iox,ierr,leqq,isgniphi,isgnbphi
       real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,zeff,
      . w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta
       real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc
@@ -790,7 +792,7 @@ c
 c
 c     signum of toroidal B and I
 c     icocos index for equilibrium from COCOS - O. Sauter Feb 2012
-      read(2,*) sgnbphi,sgniphi,icocos
+      read(2,*) isgnbphi,isgniphi,icocos
 c
 c     read parameters for analytical density and temperature
 c     profiles if iprof=0
@@ -920,16 +922,24 @@ c     read equilibrium data  from file if iequil=2
 c
       if (iequil.eq.2) then
         neqdsk=99
-        leqq=len_trim(filenmeqq)
+        call read_eqdsk(trim(filenmeqq)//'.eqdsk',rv,zv,psin,psia,
-         open(file=filenmeqq(1:leqq)//'.eqdsk',
+     .     psinr,fpol,qpsi,rcen,rmaxis,zmaxis,rbbbs,zbbbs,rlim,zlim,
-     .        status= 'unknown', unit=neqdsk)
+     .     ipsinorm,0,ipsinorm,neqdsk) !,idesc,ifreefmt,neqdsk)
+        call change_cocos(psia,fpol,qpsi,icocos,3)
+        call eq_scal(psia,fpol,isgniphi,isgnbphi,factb)
+        nlim=size(rlim)
         call equidata
-        close(neqdsk)
+c
+c     locate  btot=bres
+c
+        call bfield_res
         
 c     print density, temperature, safecty factor, toroidal current dens
 c     versus psi, rhop, rhot
         call print_prof
       end if
+      sgnbphi=isgnbphi
+      sgnbphi=isgniphi
 
       if (iequil.eq.1)   call bres_anal
 
@@ -1219,7 +1229,7 @@ c
      . rmxm,zmnm,zmxm,dr,dz,zbmin,zbmax,fpolas,phitedge,rup,zup,
      . rlw,zlw,tr,tz,tfp,cc=>cceq,cfp,cc01=>cceq01,cc10=>cceq10,
      . cc20=>cceq20,cc02=>cceq02,cc11=>cceq11,psinr,qpsi,rv,zv,psin,
-     . fpol,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd
+     . fpol,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd,read_eqdsk
       use reflections, only : rlim,zlim,nlim,alloc_lim
       implicit none
 c local constants
@@ -1235,23 +1245,12 @@ c local variables
      . rhot15,psi15,rhot2,psi2,rhotsx
       real(wp_), dimension(:), allocatable :: fpoli,fvpsi,
      . wrkf,wf,wrk
-      character(len=48) :: stringa
 c common/external functions/variables
       real(wp_) :: frhotor
 c
-c     read from file eqdsk
+      nr=size(rv)
-c     (see http://fusion.gat.com/efit/g_eqdsk.html)
+      nz=size(zv)
-c
+      nbbbs=size(rbbbs)
-c     spline interpolation of psi and derivatives
-c
-      if(icocos.gt.0) then 
-        read (neqdsk,'(a48,3i4)') stringa,idum,nr,nz
-      else
-        read (neqdsk,*) nr,nz
-      end if
-c
-      call alloc_equilvec(ierr)
-      if (ierr.ne.0) stop
       
       nrest=nr+4
       nzest=nz+4
@@ -1275,83 +1274,14 @@ c
       allocate(fvpsi(nr*nz),wf(nr),wrk(lwrk),iwrk(liwrk),
      . fpoli(nr),wrkf(lwrkf),iwrkf(nrest))
 c
-      if(ipsinorm.eq.0) then
-        read (neqdsk,2020) drnr1,dznz1,rcen,rleft,zmid
-        read (neqdsk,2020) rmaxis,zmaxis,psiaxis,psiedge,xdum
-        read (neqdsk,2020) xdum,xdum,xdum,xdum,xdum
-        read (neqdsk,2020) xdum,xdum,xdum,xdum,xdum
-        read (neqdsk,2020) (fpol(i),i=1,nr)
-        read (neqdsk,2020) (xdum,i=1,nr)
-        read (neqdsk,2020) (xdum,i=1,nr)
-        read (neqdsk,2020) (xdum,i=1,nr)
-        read (neqdsk,2020) ((psin(i,j),i=1,nr),j=1,nz)
-        read (neqdsk,2020) (qpsi(i),i=1,nr)
-      else
-        read (neqdsk,*) drnr1,dznz1,rcen,rleft,zmid
-        read (neqdsk,*) rmaxis,zmaxis,psiaxis,psiedge,xdum
-        read (neqdsk,*) xdum,xdum,xdum,xdum,xdum
-        read (neqdsk,*) xdum,xdum,xdum,xdum,xdum
-        read (neqdsk,*) (fpol(i),i=1,nr)
-        read (neqdsk,*) (xdum,i=1,nr)
-        read (neqdsk,*) (xdum,i=1,nr)
-        read (neqdsk,*) (xdum,i=1,nr)
-        read (neqdsk,*) ((psin(i,j),i=1,nr),j=1,nz)
-        read (neqdsk,*) (qpsi(i),i=1,nr)
-      end if
-2020  format (5e16.9)
-
-      if(ipsinorm.eq.0) psin=(psin-psiaxis)/(psiedge-psiaxis)
       psia0=psiedge-psiaxis
 c
-c  compensate for different reference systems
-c
-      icocosmod=mod(icocos,10)
-c
-      if (mod(icocos,2).eq.0.and.icocosmod.gt.0) then
-c  icocos mod 2 = 0: toroidal angle phi CW (opposite to gray convention)
-        do i=1,nr
-          fpol(i)=-fpol(i)
-        end do
-      end if
-c
-      if (icocosmod.eq.1 .or. icocosmod.eq.4 .or. 
-     &    icocosmod.eq.5 .or. icocosmod.eq.8) then
-c  icocos mod 10 = 1,4,5,8: psi increasing with CCW Ip
-c  icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip
-        psia0=-psia0
-      end if
-c
 c     length in m !!!
-c
-      dr=drnr1/dble(nr-1)
-      dz=dznz1/dble(nz-1)
-      rv(1)=rleft
-      zv(1)=zmid-dznz1/2.0_wp_
-      dpsinr=1.0_wp_/dble(nr-1)
-c
-      do i=1,nr
-        psinr(i)=(i-1)*dpsinr
-        rv(i)=rv(1)+(i-1)*dr
-      end do
-c
-      do j=1,nz
-        zv(j)=zv(1)+(j-1)*dz
-      end do
 c
       rmnm=rv(1)
       rmxm=rv(nr)
       zmnm=zv(1)
       zmxm=zv(nz)
-c
-c     psi function
-c
-c  icocos=0: adapt psi to force Ip sign, otherwise maintain psi
-      if (icocosmod.ne.0) sgniphi=sign(1.0_wp_,-psia0)
-c  icocos>10: input psi is in Wb
-c  icocos<10: input psi is in Wb/rad (gray convention)
-      if (icocos.ge.10) psia0=psia0/(2.0_wp_*pi)
-c
-      psia=sign(psia0,-sgniphi)*factb
 c
       do j=1,nz
         do i=1,nr
@@ -1422,28 +1352,9 @@ c
       call splev(tfp,nsf,cfp,3,psinr,fpoli,nr,ier)
       fpolas=fpoli(nr)
       btrcen=fpolas/rcen
-c
+
-c     read plasma boundaries from eqdsk file
-c
-      read (neqdsk,*) nbbbs,nlim
-c
       call alloc_bnd(ierr)
-      if (ierr.ne.0) stop
       call alloc_lim(ierr)
-      if (ierr.ne.0) stop
-c
-      if(nbbbs.gt.0) then
-        if(ipsinorm.eq.1)
-     .   read(neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs)
-        if(ipsinorm.eq.0)
-     .   read(neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs)
-      end if
-      if(nlim.gt.0) then
-        if(ipsinorm.eq.1)
-     .   read(neqdsk,*) (rlim(i),zlim(i),i=1,nlim)
-        if(ipsinorm.eq.0) 
-     .   read(neqdsk,2020) (rlim(i),zlim(i),i=1,nlim)
-      end if
 c
 c     compute max and min z of last closed surface
 c
@@ -1594,10 +1505,6 @@ c     locate  psi surface for q=1.5 and q=2
         print'(3(a,f8.5))','psi_2   = ',psi2,'  rhop_2   = '
      .   ,sqrt(psi2),'  rhot_2   = ',rhot2
       end if
-c
-c     locate  btot=bres
-c
-      call bfield_res
 c
       deallocate(iwrk,iwrkf,fpoli,fvpsi,wrkf,wf,wrk)
       return

src/interp_eqprof.f90

@@ -255,4 +255,203 @@ contains
     if (present(bz)) bz= bz/rpsim
   end subroutine bfield
   
+  subroutine eq_scal(psia,fpol,isign,bsign,factb)
+    use const_and_precisions, only : one
+    implicit none
+    real(wp_), intent(inout) :: psia
+    integer, intent(inout) :: isign,bsign
+    real(wp_), dimension(:), intent(inout) :: fpol
+    real(wp_), intent(in) :: factb
+    
+    ! isign and bsign ignored on input if equal to 0
+    ! used to assign the direction of Bphi and Ipla BEFORE scaling
+    ! cocos=3 assumed: CCW direction is >0
+    ! Bphi and Ipla scaled by the same factor factb to keep q unchanged
+    ! factb<0 reverses the directions of Bphi and Ipla
+    if(isign/=0) psia=sign(psia,real(-isign,wp_))
+    if(bsign/=0 .and. bsign*fpol(size(fpol))<0) fpol=-fpol
+    psia=psia*factb
+    fpol=fpol*factb
+    isign=int(sign(one,-psia))
+    bsign=int(sign(one,fpol(size(fpol))))
+  end subroutine eq_scal
+
+  subroutine read_eqdsk(filenm,rv,zv,psin,psia,psinr,fpol,q,rvac,rax,zax, &
+                        rbnd,zbnd,rlim,zlim,ipsinorm,idesc,ifreefmt,unit)
+    use const_and_precisions, only : one
+    use utils, only : get_free_unit
+    implicit none
+! arguments
+    character(len=*), intent(in) :: filenm
+    real(wp_), intent(out) :: psia,rvac,rax,zax
+    real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,psinr,fpol,q
+    real(wp_), dimension(:), allocatable, intent(out) :: rbnd,zbnd,rlim,zlim
+    real(wp_), dimension(:,:), allocatable, intent(out) :: psin
+    integer, optional, intent(in) :: ipsinorm,idesc,ifreefmt,unit
+! local variables
+    integer, parameter :: indef=0,iddef=1,iffdef=0
+    integer :: in,id,iffmt,u,idum,i,j,nr,nz,nbnd,nlim
+    character(len=48) :: string
+    real(wp_) :: dr,dz,dps,rleft,zmid,zleft,xdum,psiedge,psiaxis
+    
+! set default values if optional arguments are absent
+    in=indef; id=iddef; iffmt=iffdef
+    if(present(ipsinorm)) in=ipsinorm
+    if(present(idesc)) id=idesc
+    if(present(ifreefmt)) iffmt=ifreefmt
+    if (present(unit)) then
+      u=unit
+    else
+      u=get_free_unit()
+    end if
+
+! open G EQDSK file (see http://fusion.gat.com/efit/g_eqdsk.html)
+    open(file=trim(filenm),status='old',unit=u)
+
+! get size of main arrays and allocate them
+    if (id==1) then
+      read (u,'(a48,3i4)') string,idum,nr,nz
+    else
+      read (u,*) nr,nz
+    end if
+    if (allocated(rv)) deallocate(rv)
+    if (allocated(zv)) deallocate(zv)
+    if (allocated(psin)) deallocate(psin)
+    if (allocated(psinr)) deallocate(psinr)
+    if (allocated(fpol)) deallocate(fpol)
+    if (allocated(q)) deallocate(q)
+    allocate(rv(nr),zv(nz),psin(nr,nz),psinr(nr),fpol(nr),q(nr))
+
+! store 0D data and main arrays
+    if (iffmt==1) then
+      read (u,*) dr,dz,rvac,rleft,zmid
+      read (u,*) rax,zax,psiaxis,psiedge,xdum
+      read (u,*) xdum,xdum,xdum,xdum,xdum
+      read (u,*) xdum,xdum,xdum,xdum,xdum
+      read (u,*) (fpol(i),i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) ((psin(i,j),i=1,nr),j=1,nz)
+      read (u,*) (q(i),i=1,nr)
+    else
+      read (u,'(5e16.9)') dr,dz,rvac,rleft,zmid
+      read (u,'(5e16.9)') rax,zax,psiaxis,psiedge,xdum
+      read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum
+      read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum
+      read (u,'(5e16.9)') (fpol(i),i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') ((psin(i,j),i=1,nr),j=1,nz)
+      read (u,'(5e16.9)') (q(i),i=1,nr)
+    end if
+
+! get size of boundary and limiter arrays and allocate them
+    read (u,*) nbnd,nlim
+    if (allocated(rbnd)) deallocate(rbnd)
+    if (allocated(zbnd)) deallocate(zbnd)
+    if (allocated(rlim)) deallocate(rlim)
+    if (allocated(zlim)) deallocate(zlim)
+    allocate(rbnd(nbnd),zbnd(nbnd),rlim(nlim),zlim(nlim))
+
+! store boundary and limiter data
+    if(nbnd>0) then
+      if (iffmt==1) then
+        read(u,*) (rbnd(i),zbnd(i),i=1,nbnd)
+      else
+        read(u,'(5e16.9)') (rbnd(i),zbnd(i),i=1,nbnd)
+      end if
+    end if
+    if(nlim>0) then
+      if (iffmt==1) then
+        read(u,*) (rlim(i),zlim(i),i=1,nlim)
+      else
+        read(u,'(5e16.9)') (rlim(i),zlim(i),i=1,nlim)
+      end if
+    end if
+
+! reading of G EQDSK file completed    
+    close(u)
+    
+! build rv,zv,psinr arrays and normalize psin
+    zleft=zmid-0.5_wp_*dz
+    dr=dr/(nr-1)
+    dz=dz/(nz-1)
+    dps=one/(nr-1)
+    do i=1,nr
+      psinr(i)=(i-1)*dps
+      rv(i)=rleft+(i-1)*dr
+    end do
+    do i=1,nz
+      zv(i)=zleft+(i-1)*dz
+    end do
+    psia=psiedge-psiaxis
+    if(in==0) psin=(psin-psiaxis)/psia
+
+  end subroutine read_eqdsk
+  
+  subroutine change_cocos(psia,fpol,q,cocosin,cocosout,ierr)
+    use const_and_precisions, only : zero,one,pi
+    implicit none
+! arguments
+    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
+    integer :: exp2pi,exp2piout,npsi
+    logical :: phiccw,psiincr,qpos,phiccwout,psiincrout,qposout
+
+    call decode_cocos(cocosin,exp2pi,phiccw,psiincr,qpos)
+    call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout)
+
+!   check sign consistency
+    isign=sign(one,psia)
+    if (.not.psiincr) isign=-isign
+    bsign=sign(one,fpol(size(fpol)))
+    if (qpos.neqv.isign*bsign*q(size(q))>zero) then
+!     warning: sign inconsistency found among q, Ipla and Bref
+      q=-q
+      if(present(ierr)) ierr=1
+    else
+      if(present(ierr)) ierr=0
+    end if
+
+! convert cocosin to cocosout
+    if (phiccw.neqv.phiccwout) then
+!     opposite direction of toroidal angle phi in cocosin and cocosout
+!      bsign=-bsign
+!      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
+    if ((phiccw.eqv.phiccwout) .neqv. (psiincr.eqv.psiincrout)) psia=-psia
+!   convert Wb to Wb/rad or viceversa
+    psia=psia*(2.0_wp_*pi)**(exp2piout-exp2pi)
+  end subroutine change_cocos
+
+  subroutine decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos)
+    implicit none
+    integer, intent(in) :: cocos
+    integer, intent(out) :: exp2pi
+    logical, intent(out) :: phiccw,psiincr,qpos
+    integer :: cmod10,cmod4
+
+    cmod10=mod(cocos,10)
+    cmod4=mod(cmod10,4)
+!  cocos>10 psi in Wb, cocos<10 psi in Wb/rad
+    exp2pi=(cocos-cmod10)/10
+!  cocos mod 10 = 1,3,5,7: toroidal angle phi increasing CCW
+    phiccw=(mod(cmod10,2)==1)
+!  cocos mod 10 = 1,2,5,6: psi increasing with positive Ip
+    psiincr=(cmod4==1 .or. cmod4==2)
+!  cocos mod 10 = 1,2,7,8: q positive for positive Bphi*Ip
+    qpos=(cmod10<3 .or. cmod10>6)
+  end subroutine decode_cocos
+
 end module interp_eqprof

src/reflections.f90

@@ -35,7 +35,7 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
   real(wp_), dimension(nw), intent(in) :: rw,zw
   real(wp_), intent(out) :: sint
   real(wp_), dimension(3), intent(out) :: normw
-  integer :: i,j,ni,iint
+  integer :: i,j,ni,iint,nneg
   real(wp_), dimension(2) :: si,ti
   real(wp_) :: drw,dzw,xint,yint,rint,l,kxy
   real(wp_) :: tol
@@ -46,13 +46,21 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
   do i=1,nw-1
     !search intersections with i-th wall segment
     call linecone_coord(xv,kv,rw(i:i+1),zw(i:i+1),si,ti,ni)
-    do while (ni>0 .and. si(1)<=tol)
+    !discard solutions with s<=0
-      !remove solutions with s<=0
+    nneg=0
-      ni = ni-1
-      si(1) = si(2)
-      ti(1) = ti(2)
-    end do
     do j=1,ni
+      if (si(j)<=tol) then
+        nneg=j
+      else
+        exit
+      end if
+    end do
+!    do while (ni>0 .and. si(1)<=tol)
+!      ni = ni-1
+!      si(1) = si(2)   ???
+!      ti(1) = ti(2)   ???
+!    end do
+    do j=nneg+1,ni
       if ((si(j)<sint .or. iint==0) .and. ti(j)>=zero .and. ti(j)<=one) then
         !check intersection is in r,z range and keep the closest
         sint = si(j)