From 2856725b4959ad5d6eb7b50f063a3aa5549d2a01 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Fri, 18 Sep 2015 10:24:24 +0000 Subject: [PATCH] partial changes uploaded by mistake in the previous revision. gray, interp_eqprof and reflections were in an inconsistent state --- src/gray.f | 139 ++++++++++++++++++++++++----- src/interp_eqprof.f90 | 199 ------------------------------------------ src/reflections.f90 | 22 ++--- 3 files changed, 123 insertions(+), 237 deletions(-) diff --git a/src/gray.f b/src/gray.f index c9e5e5e..4f79a44 100644 --- a/src/gray.f +++ b/src/gray.f @@ -658,14 +658,12 @@ 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,psia, - . rmaxis,zmaxis,fpolas,psinr,qpsi,rv,zv,psin, - . fpol,rbbbs,zbbbs,read_eqdsk,change_cocos,eq_scal + use interp_eqprof, only : rmxm,zbmin,zbmax,btrcen,rcen use reflections, only : rlim,zlim,nlim,alloc_lim use beamdata, only : nrayr,nrayth,nstep implicit none c local variables - integer :: nfil,iox,ierr,leqq,isgniphi,isgnbphi + integer :: nfil,iox,ierr,leqq real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,zeff, . w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc @@ -792,7 +790,7 @@ c c c signum of toroidal B and I c icocos index for equilibrium from COCOS - O. Sauter Feb 2012 - read(2,*) isgnbphi,isgniphi,icocos + read(2,*) sgnbphi,sgniphi,icocos c c read parameters for analytical density and temperature c profiles if iprof=0 @@ -922,24 +920,16 @@ c read equilibrium data from file if iequil=2 c if (iequil.eq.2) then 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) - call change_cocos(psia,fpol,qpsi,icocos,3) - call eq_scal(psia,fpol,isgniphi,isgnbphi,factb) - nlim=size(rlim) + leqq=len_trim(filenmeqq) + open(file=filenmeqq(1:leqq)//'.eqdsk', + . status= 'unknown', unit=neqdsk) call equidata -c -c locate btot=bres -c - call bfield_res + close(neqdsk) 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 @@ -1229,7 +1219,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,read_eqdsk + . fpol,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd use reflections, only : rlim,zlim,nlim,alloc_lim implicit none c local constants @@ -1245,13 +1235,24 @@ 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 - nr=size(rv) - nz=size(zv) - nbbbs=size(rbbbs) - +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 + nrest=nr+4 nzest=nz+4 lwrk=4*(nr+nz)+11*(nrest+nzest)+nrest*nz+nzest+54 @@ -1274,14 +1275,83 @@ 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 @@ -1352,9 +1422,28 @@ 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 @@ -1505,6 +1594,10 @@ 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 055e26a..b8dbd65 100644 --- a/src/interp_eqprof.f90 +++ b/src/interp_eqprof.f90 @@ -255,203 +255,4 @@ 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 42f5fa6..567db67 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,nneg + integer :: i,j,ni,iint real(wp_), dimension(2) :: si,ti real(wp_) :: drw,dzw,xint,yint,rint,l,kxy real(wp_) :: tol @@ -46,21 +46,13 @@ 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) - !discard solutions with s<=0 - nneg=0 - do j=1,ni - if (si(j)<=tol) then - nneg=j - else - exit - end if + 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 -! 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 + do j=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)