subroutine wall_refl moved in reflections module, rlim,zlim,nlim moved from interp_eqprof to reflections module

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,
+      use interp_eqprof, only : rmxm,zbmin,zbmax,btrcen,rcen
-     . btrcen,rcen,alloc_lim
+      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_
+      rcon(1,:)=rmaxis
-      zcon(1,:)=0.0_wp_
+      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_
+      rcon(1,:)=rmaxis
-      zcon(1,:)=0.0_wp_
+      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)<tiny(walln)) then
-          ! wall never hit
-          xvrfl=xv
-          anvrfl=anv0
-          extr=ext
-          eytr=eyt
-          irfl=0
-          return
-        end if
-        ! search second wall interface (inside-outside)
-        call inters_linewall(xvrfl/1.0e2_wp_,anv0,rlim(1:nlim),
-     .                       zlim(1:nlim),nlim,smax,walln)
-        smax=smax*1.0e2_wp_
-        xvrfl=xvrfl+smax*anv0
-        irfl=2
-      end if
-c
-c     rotation matrix from local to lab frame
-      vv1(1)=anv0(2)
-      vv1(2)=-anv0(1)
-      vv1(3)=0.0_wp_
-      vv2(1)=anv0(1)*anv0(3)
-      vv2(2)=anv0(2)*anv0(3)
-      vv2(3)=-anv0(1)*anv0(1)-anv0(2)*anv0(2)
-      vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
-      vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
-      vv3=anv0
-c
-      evin=ext*vv1+eyt*vv2
-c     wave vector and electric field after reflection in lab frame
-      anvrfl=anv0-2.0_wp_*
-     .       (anv0(1)*walln(1)+anv0(2)*walln(2)+anv0(3)*walln(3))*walln
-      evrfl=-evin+2.0_wp_*
-     .       (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln
-c
-      vv1(1)=anvrfl(2)
-      vv1(2)=-anvrfl(1)
-      vv1(3)=0.0_wp_
-      vv2(1)=anvrfl(1)*anvrfl(3)
-      vv2(2)=anvrfl(2)*anvrfl(3)
-      vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2)
-      vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
-      vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
-      vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2)
-c
-      extr=dot_product(vv1,evrfl)
-      eytr=dot_product(vv2,evrfl)
-      eztr=dot_product(vv3,evrfl)
-c
-      return
-      end
-c
-c
 c
       subroutine vacuum_rt(xvstart,anv,xvend,ivac)
       use const_and_precisions, only : wp_
-      use reflections, only : inters_linewall,inside
+      use reflections, only : inters_linewall,inside,rlim,zlim,nlim
       use graydata_flags, only : dst
-      use interp_eqprof, only : rlim,zlim,nlim
       implicit none
 c arguments
       integer :: ivac

src/interp_eqprof.f90

@@ -14,9 +14,6 @@ module interp_eqprof
   ! === 1D array plasma boundary Rbnd_i, Zbnd_i ====================
   INTEGER, SAVE :: nbbbs
   REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rbbbs,zbbbs
-  ! === 1D array limiter Rlim_i, Zlim_i   ==> 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

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)<tiny(walln)) then
+      ! wall never hit
+      xvrfl=xv
+      anvrfl=anv0
+      extr=ext
+      eytr=eyt
+      irfl=0
+      return
+    end if
+    ! search second wall interface (inside-outside)
+    call inters_linewall(xvrfl/1.0e2_wp_,anv0,rlim(1:nlim), &
+                         zlim(1:nlim),nlim,smax,walln)
+    smax=smax*1.0e2_wp_
+    xvrfl=xvrfl+smax*anv0
+    irfl=2
+  end if
+!
+! rotation matrix from local to lab frame
+  vv1(1)=anv0(2)
+  vv1(2)=-anv0(1)
+  vv1(3)=0.0_wp_
+  vv2(1)=anv0(1)*anv0(3)
+  vv2(2)=anv0(2)*anv0(3)
+  vv2(3)=-anv0(1)*anv0(1)-anv0(2)*anv0(2)
+  vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
+  vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
+  vv3=anv0
+!
+  evin=ext*vv1+eyt*vv2
+! wave vector and electric field after reflection in lab frame
+  anvrfl=anv0-2.0_wp_* &
+         (anv0(1)*walln(1)+anv0(2)*walln(2)+anv0(3)*walln(3))*walln
+  evrfl=-evin+2.0_wp_* &
+         (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln
+!
+  vv1(1)=anvrfl(2)
+  vv1(2)=-anvrfl(1)
+  vv1(3)=0.0_wp_
+  vv2(1)=anvrfl(1)*anvrfl(3)
+  vv2(2)=anvrfl(2)*anvrfl(3)
+  vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2)
+  vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
+  vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
+  vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2)
+!
+  extr=dot_product(vv1,evrfl)
+  eytr=dot_product(vv2,evrfl)
+  eztr=dot_product(vv3,evrfl)
+end
+
 end module reflections