From 667f6fd11108b2ab2e6c7940249bab9769980e6d Mon Sep 17 00:00:00 2001 From: Daniele Micheletti Date: Thu, 16 Jul 2015 14:47:56 +0000 Subject: [PATCH] subroutine wall_refl moved in reflections module, rlim,zlim,nlim moved from interp_eqprof to reflections module --- src/gray.f | 106 ++++++------------------------------------ src/interp_eqprof.f90 | 32 +------------ src/reflections.f90 | 102 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 117 insertions(+), 123 deletions(-) diff --git a/src/gray.f b/src/gray.f index e555630..10d025c 100644 --- a/src/gray.f +++ b/src/gray.f @@ -203,10 +203,10 @@ c c subroutine after_onestep(i,istop) use const_and_precisions, only : wp_,pi - use reflections, only : inside + use reflections, only : inside,rlim,zlim,nlim,wall_refl use graydata_flags, only : iwarm,istpr0,istpl0,dst,ipass,igrad use graydata_par, only : psipol0,chipol0 - use interp_eqprof, only : zbmin,zbmax,rlim,zlim,nlim + use interp_eqprof, only : zbmin,zbmax use beamdata, only : nrayr,nrayth,psjki,ppabs,ccci,iiv,tauv,ihcd, . istore,iop,iow,tau1v,yyrfl,ywrk,ext,eyt implicit none @@ -658,8 +658,8 @@ c use graydata_par use graydata_anequil use coreprofiles, only : psdbnd - use interp_eqprof, only : rmxm,rlim,zlim,nlim,zbmin,zbmax, - . btrcen,rcen,alloc_lim + 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 @@ -1213,12 +1213,13 @@ c use dierckx, only : curfit,splev,regrid,bispev,coeff_parder use graydata_flags, only : ipsinorm,sspl,ixp,icocos,neqdsk use graydata_par, only : sgnbphi,sgniphi,factb - use interp_eqprof, only : nsr,nsz,nsf,rlim,zlim,nlim,nr,nz, + use interp_eqprof, only : nsr,nsz,nsf,nr,nz, . psia,psiant,psinop,btrcen,rcen,btaxis,rmaxis,zmaxis,rmnm, . 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 + use reflections, only : rlim,zlim,nlim,alloc_lim implicit none c local constants integer, parameter :: kspl=3 @@ -1427,6 +1428,8 @@ c 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) @@ -2551,8 +2554,8 @@ c computation of flux surface averaged quantities fc=1.0_wp_ psicon(1)=0.0_wp_ - rcon(1,:)=0.0_wp_ - zcon(1,:)=0.0_wp_ + rcon(1,:)=rmaxis + zcon(1,:)=zmaxis pstab(1)=0.0_wp_ rhot_eq(1)=0.0_wp_ rpstab(1)=0.0_wp_ @@ -2965,8 +2968,8 @@ c computation of flux surface averaged quantities fc=1.0_wp_ psicon(1)=0.0_wp_ - rcon(1,:)=0.0_wp_ - zcon(1,:)=0.0_wp_ + rcon(1,:)=rmaxis + zcon(1,:)=zmaxis pstab(1)=0.0_wp_ rhot_eq(1)=0.0_wp_ rpstab(1)=0.0_wp_ @@ -5976,94 +5979,11 @@ c end function inside_plasma c c -c - subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln, - . irfl) - use const_and_precisions, only : wp_,ui=>im,pi - use reflections, only : inters_linewall,inside - use interp_eqprof, only : rlim,zlim,nlim - implicit none -c arguments - integer :: irfl - real(wp_), dimension(3) :: xv,anv,xvrfl,anvrfl,walln - complex(wp_) :: ext,eyt,extr,eytr -c local variables - real(wp_) :: smax,rrm,zzm - real(wp_), dimension(3) :: anv0,vv1,vv2,vv3 - complex(wp_) :: eztr - complex(wp_), dimension(3) :: evin,evrfl -c - anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) - rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2) - zzm=1.0e-2_wp_*xv(3) -c -c computation of reflection coordinates and normal to the wall - call inters_linewall(xv/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), - . nlim,smax,walln) - smax=smax*1.0e2_wp_ - xvrfl=xv+smax*anv0 - irfl=1 - if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then - ! first wall interface is outside-inside - if (dot_product(walln,walln) move in wall/reflections - INTEGER, SAVE :: nlim - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rlim,zlim REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis REAL(wp_), SAVE :: rmnm,rmxm,zmnm,zmxm,dr,dz @@ -90,52 +87,27 @@ contains end subroutine dealloc_equilvec - subroutine alloc_bnd(ier) implicit none integer, intent(out) :: ier - if(nlim.lt.0.or.nbbbs.lt.0) then + if(nbbbs.lt.0) then ier = -1 return end if call dealloc_bnd - allocate(rlim(nlim),zlim(nlim),rbbbs(nbbbs),zbbbs(nbbbs), & + allocate(rbbbs(nbbbs),zbbbs(nbbbs), & stat=ier) if (ier/=0) call dealloc_bnd end subroutine alloc_bnd subroutine dealloc_bnd implicit none - if(allocated(rlim)) deallocate(rlim) - if(allocated(zlim)) deallocate(zlim) if(allocated(rbbbs)) deallocate(rbbbs) if(allocated(zbbbs)) deallocate(zbbbs) end subroutine dealloc_bnd - - subroutine alloc_lim(ier) - implicit none - integer, intent(out) :: ier - - if(nlim.lt.0) then - ier = -1 - return - end if - - call dealloc_lim - allocate(rlim(nlim),zlim(nlim), & - stat=ier) - if (ier/=0) call dealloc_lim - end subroutine alloc_lim - - subroutine dealloc_lim - implicit none - if(allocated(rlim)) deallocate(rlim) - if(allocated(zlim)) deallocate(zlim) - end subroutine dealloc_lim - subroutine equinum_psi(rpsim,zpsim,psinv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz) use dierckx, only : fpbisp diff --git a/src/reflections.f90 b/src/reflections.f90 index 8c7b708..567db67 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -1,9 +1,16 @@ module reflections use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge, zero, one implicit none + + ! === 1D array limiter Rlim_i, Zlim_i + integer, public, save :: nlim + real(wp_), public, dimension(:), allocatable, save :: rlim,zlim + private public :: reflect,inters_linewall,inside public :: linecone_coord,interssegm_coord,interssegm + public :: alloc_lim,wall_refl + contains subroutine reflect(ki,nsurf,ko) @@ -205,5 +212,100 @@ function inside(xc,yc,n,x,y) inside=(mod(locatef(xint,nj,x),2)==1) end function inside +subroutine alloc_lim(ier) + implicit none + integer, intent(out) :: ier + + if(nlim.lt.0) then + ier = -1 + return + end if + + call dealloc_lim + allocate(rlim(nlim),zlim(nlim), & + stat=ier) + if (ier/=0) call dealloc_lim +end subroutine alloc_lim + +subroutine dealloc_lim + implicit none + if(allocated(rlim)) deallocate(rlim) + if(allocated(zlim)) deallocate(zlim) +end subroutine dealloc_lim + +subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl) + implicit none +! arguments + integer :: irfl + real(wp_), dimension(3) :: xv,anv,xvrfl,anvrfl,walln + complex(wp_) :: ext,eyt,extr,eytr +! local variables + real(wp_) :: smax,rrm,zzm + real(wp_), dimension(3) :: anv0,vv1,vv2,vv3 + complex(wp_) :: eztr + complex(wp_), dimension(3) :: evin,evrfl +! + anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) + rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2) + zzm=1.0e-2_wp_*xv(3) +! +! computation of reflection coordinates and normal to the wall + call inters_linewall(xv/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), & + nlim,smax,walln) + smax=smax*1.0e2_wp_ + xvrfl=xv+smax*anv0 + irfl=1 + if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then + ! first wall interface is outside-inside + if (dot_product(walln,walln)