From c06fbf3d4f5115add279076a58a8a0f389e36c00 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 17 Sep 2015 15:46:21 +0000 Subject: [PATCH] cleaned arguments list in diel_tens_fr --- src/dispersion.f90 | 7 +- src/gray.f | 139 +++++------------------------ src/interp_eqprof.f90 | 199 ++++++++++++++++++++++++++++++++++++++++++ src/reflections.f90 | 22 +++-- 4 files changed, 241 insertions(+), 126 deletions(-) diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 3898575..f199454 100644 --- a/src/dispersion.f90 +++ b/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 diff --git a/src/gray.f b/src/gray.f index 4f79a44..c9e5e5e 100644 --- a/src/gray.f +++ b/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) - open(file=filenmeqq(1:leqq)//'.eqdsk', - . status= 'unknown', unit=neqdsk) + 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) + 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,24 +1245,13 @@ 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 -c (see http://fusion.gat.com/efit/g_eqdsk.html) -c -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 - + nr=size(rv) + nz=size(zv) + nbbbs=size(rbbbs) + nrest=nr+4 nzest=nz+4 lwrk=4*(nr+nz)+11*(nrest+nzest)+nrest*nz+nzest+54 @@ -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 diff --git a/src/interp_eqprof.f90 b/src/interp_eqprof.f90 index b8dbd65..055e26a 100644 --- a/src/interp_eqprof.f90 +++ b/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 diff --git a/src/reflections.f90 b/src/reflections.f90 index 567db67..42f5fa6 100644 --- a/src/reflections.f90 +++ b/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) - !remove solutions with s<=0 - ni = ni-1 - si(1) = si(2) - ti(1) = ti(2) - end do + !discard solutions with s<=0 + nneg=0 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)=zero .and. ti(j)<=one) then !check intersection is in r,z range and keep the closest sint = si(j)