From baf53b932b0e48a43f647bf9ab5bb24a5f835f6e Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Tue, 23 Apr 2024 17:00:40 +0200 Subject: [PATCH] simplify memory management This change replaces pointers with automatic arrays to greatly simplify the memory management in the main subroutine: - All arrays are defined in a single location and with their final dimension explicitely shown. - The allocation/deallocation is performed automatically when entering/leaving the gray_main routine. --- src/beamdata.f90 | 101 +--------------- src/eccd.f90 | 8 +- src/gray_core.f90 | 91 ++++++++------- src/gray_params.f90 | 5 + src/main.f90 | 23 +--- src/multipass.f90 | 272 +++++++++++++++++--------------------------- 6 files changed, 169 insertions(+), 331 deletions(-) diff --git a/src/beamdata.f90 b/src/beamdata.f90 index 652d750..03fe358 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -8,18 +8,10 @@ module beamdata contains - subroutine init_btr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) + subroutine init_btr(rtrparam) use gray_params, only : raytracing_parameters use const_and_precisions, only : half,two - type(raytracing_parameters), intent(inout) :: rtrparam - real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, & - gri,psjki,ppabs,ccci - real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri - real(wp_), dimension(:), intent(out), pointer :: tau0,alphaabs0, & - dids0,ccci0,p0jk - complex(wp_), dimension(:), intent(out), pointer :: ext, eyt - integer, dimension(:), intent(out), pointer :: iiv + type(raytracing_parameters), intent(in) :: rtrparam integer :: jray1 @@ -41,7 +33,6 @@ contains end if nray = (nrayr-1)*nrayth + 1 jkray1 = (jray1-2)*nrayth + 2 - rtrparam%nray = nray if(nrayr>1) then twodr2 = two*(rwmax/(nrayr-1))**2 @@ -50,13 +41,6 @@ contains end if nstep=rtrparam%nstep - - ! Allocate for each ray: - ! y = (x, k), - ! yp = dy/dt, - ! etc. - call alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) end subroutine init_btr @@ -134,85 +118,4 @@ contains end if end function rayi2j - - - function rayi2k(i) result(kt) - integer, intent(in) :: i - integer :: kt - -! kt = max(1, mod(i-2,nrayth) + 1) - if (i>1) then - kt = mod(i-2,nrayth) + 1 - else - kt = 1 - end if - end function rayi2k - - - - function rayjk2i(jr,kt) result(i) - integer, intent(in) :: jr,kt - integer :: i - -! i = max(1, (jr-2)*nrayth + kt + 1) - if (jr>1) then - i = (jr-2)*nrayth + kt + 1 - else - i = 1 - end if - end function rayjk2i - - - - subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, & - gri,psjki,ppabs,ccci - real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri - real(wp_), dimension(:), intent(out), pointer :: tau0,alphaabs0, & - dids0,ccci0,p0jk - complex(wp_), dimension(:), intent(out), pointer :: ext, eyt - integer, dimension(:), intent(out), pointer :: iiv - - call dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - - allocate(ywork(6,nray), ypwork(6,nray), gri(3,nray), ggri(3,3,nray), & - xc(3,nrayth,nrayr), du1(3,nrayth,nrayr), & - psjki(nray,nstep), ppabs(nray,nstep), ccci(nray,nstep), & - tau0(nray), alphaabs0(nray), dids0(nray), ccci0(nray), & - p0jk(nray), ext(nray), eyt(nray), iiv(nray)) - end subroutine alloc_beam - - - - subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & - tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) - real(wp_), dimension(:,:), intent(out), pointer :: ywork,ypwork, & - gri,psjki,ppabs,ccci - real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri - real(wp_), dimension(:), intent(out), pointer :: tau0,alphaabs0, & - dids0,ccci0,p0jk - complex(wp_), dimension(:), intent(out), pointer :: ext, eyt - integer, dimension(:), intent(out), pointer :: iiv - - if (associated(ywork)) deallocate(ywork) - if (associated(ypwork)) deallocate(ypwork) - if (associated(xc)) deallocate(xc) - if (associated(du1)) deallocate(du1) - if (associated(gri)) deallocate(gri) - if (associated(ggri)) deallocate(ggri) - if (associated(psjki)) deallocate(psjki) - if (associated(ppabs)) deallocate(ppabs) - if (associated(ccci)) deallocate(ccci) - if (associated(tau0)) deallocate(tau0) - if (associated(alphaabs0)) deallocate(alphaabs0) - if (associated(dids0)) deallocate(dids0) - if (associated(ccci0)) deallocate(ccci0) - if (associated(p0jk)) deallocate(p0jk) - if (associated(ext)) deallocate(ext) - if (associated(eyt)) deallocate(eyt) - if (associated(iiv)) deallocate(iiv) - end subroutine dealloc_beam - end module beamdata diff --git a/src/eccd.f90 b/src/eccd.f90 index 9d335b7..de585ad 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -104,7 +104,7 @@ contains subroutine setcdcoeff_notrap(zeff,cst2,eccdpar) real(wp_), intent(in) :: zeff real(wp_), intent(out) :: cst2 - real(wp_), dimension(:), pointer, intent(out) :: eccdpar + real(wp_), dimension(:), allocatable, intent(out) :: eccdpar cst2=0.0_wp_ allocate(eccdpar(1)) @@ -122,7 +122,7 @@ contains use conical, only : fconic real(wp_), intent(in) :: zeff,rbn,rbx real(wp_), intent(out) :: cst2 - real(wp_), dimension(:), pointer, intent(out) :: eccdpar + real(wp_), dimension(:), allocatable, intent(out) :: eccdpar real(wp_) :: alams,pa,fp0s cst2=1.0_wp_-rbx @@ -145,7 +145,7 @@ contains integer, parameter :: ksp=3 real(wp_), intent(in) :: zeff,rbx,fc,amu,rhop real(wp_), intent(out) :: cst2 - real(wp_), dimension(:), pointer, intent(out) :: eccdpar + real(wp_), dimension(:), allocatable, intent(out) :: eccdpar real(wp_), dimension(nlmt) :: chlm integer :: nlm,ierr,npar @@ -195,7 +195,7 @@ contains rdu,rdu2t,duu,uu1,uu2,xx1,xx2,resj,resp,eji,epp,anum,denom, & cstrdut,anucc real(wp_), dimension(lw) :: w - real(wp_), dimension(:), pointer :: apar=>null() + real(wp_), dimension(:), allocatable :: apar real(wp_), dimension(0:1) :: uleft,uright ! common/external functions/variables real(wp_), external :: fcur diff --git a/src/gray_core.f90 b/src/gray_core.f90 index ead7e1d..7989a32 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -18,13 +18,12 @@ contains use beamdata, only : pweight, rayi2jk use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd use magsurf_data, only : flux_average, dealloc_surfvec - use beamdata, only : init_btr, dealloc_beam + use beamdata, only : init_btr use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & rhop_tab, rhot_tab use utils, only : vmaxmin, inside - use multipass, only : alloc_multipass, dealloc_multipass, initbeam, & - initmultipass, turnoffray, plasma_in, plasma_out, & - wall_out + use multipass, only : initbeam, initmultipass, turnoffray, & + plasma_in, plasma_out, wall_out use logger, only : log_info, log_debug ! subroutine arguments @@ -50,39 +49,56 @@ contains real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx real(wp_) :: pabs_beam,icd_beam,cpl_beam1,cpl_beam2,cpl_cbeam1,cpl_cbeam2 - real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl - real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg - - ! Ray variables - real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null() - real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null() - - ! i: integration step, jk: global ray index - integer :: i, jk - integer :: iox,nharm,nhf,nnd,iokhawa,ierrn,ierrhcd, index_rt, parent_index_rt integer :: ip,ib,iopmin,ipar,child_index_rt integer :: igrad_b,istop_pass,nbeam_pass,nlim logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff - real(wp_), dimension(:,:,:), pointer :: yynext=>null(),yypnext=>null() - real(wp_), dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null() - real(wp_), dimension(:,:), pointer :: taus=>null(),stnext=>null(), & - yw0=>null(),ypw0=>null(),cpls=>null() - real(wp_), dimension(:), pointer :: p0ray=>null(),tau0=>null(),alphaabs0=>null(), & - dids0=>null(),ccci0=>null(),tau1=>null(),etau1=>null(),cpl1=>null(),lgcpl1=>null() - real(wp_), dimension(:), pointer :: p0jk=>null() - real(wp_), dimension(:), pointer :: jphi_beam=>null(),pins_beam=>null(), & - currins_beam=>null(), dpdv_beam=>null(),jcd_beam=>null(),stv=>null(), & - psipv=>null(),chipv=>null() - complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null() - integer, dimension(:), pointer :: iiv=>null(),iop=>null(),iow=>null() - logical, dimension(:), pointer :: iwait=>null() - logical, dimension(:,:), pointer :: iroff=>null() + ! i: integration step, jk: global ray index + integer :: i, jk ! buffer for formatting log messages character(256) :: msg + real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl + real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg + + associate ( nray => params%raytracing%nray, & + nrayr => params%raytracing%nrayr, & + nrayth => params%raytracing%nrayth, & + nstep => params%raytracing%nstep, & + nbeam_tot => 2**(params%raytracing%ipass+1)-2, & + nbeam_max => 2**(params%raytracing%ipass)) + block + + ! ray variables + real(wp_), dimension(6, nray) :: yw, ypw + real(wp_), dimension(3, nray) :: gri + real(wp_), dimension(3, 3, nray) :: ggri + real(wp_), dimension(3, nrayth, nrayr) :: xc, du1 + real(wp_), dimension(nray) :: ccci0, p0jk + real(wp_), dimension(nray) :: tau0, alphaabs0 + real(wp_), dimension(nray) :: dids0 + integer, dimension(nray) :: iiv + real(wp_), dimension(nray, nstep) :: psjki, ppabs, ccci + complex(wp_), dimension(nray) :: ext, eyt + + ! multipass variables + logical, dimension(nray) :: iwait + integer, dimension(nray) :: iow, iop + logical, dimension(nray, nbeam_tot) :: iroff + real(wp_), dimension(6, nray, nbeam_max-2) :: yynext, yypnext + real(wp_), dimension(6, nray) :: yw0, ypw0 + real(wp_), dimension(nray, nbeam_tot) :: stnext + real(wp_), dimension(nray) :: stv, p0ray + real(wp_), dimension(nray, nbeam_tot) :: taus, cpls + real(wp_), dimension(nray) :: tau1, etau1, cpl1, lgcpl1 + real(wp_), dimension(0:nbeam_tot) :: psipv, chipv + + ! beam variables + real(wp_), dimension(params%output%nrho) :: jphi_beam, pins_beam + real(wp_), dimension(params%output%nrho) :: currins_beam, dpdv_beam, jcd_beam + ! ======== set environment BEGIN ======== ! Number of limiter contourn points nlim = size(data%equilibrium%zlim) @@ -95,8 +111,7 @@ contains call launchangles2n(params%antenna, anv0) ! Initialise the ray variables (beamtracing) - call init_btr(params%raytracing, yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, & - tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv) + call init_btr(params%raytracing) ! Initialise the dispersion module if(params%ecrh_cd%iwarm > 1) call expinit @@ -110,10 +125,6 @@ contains nnd = size(rhop_tab) ! number of radial profile points end if - call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & - stv, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, & - pins_beam, currins_beam, dpdv_beam, jcd_beam, psipv, chipv) - ! Allocate memory for the results... allocate(results%dpdv(params%output%nrho)) allocate(results%jcd(params%output%nrho)) @@ -612,13 +623,11 @@ contains ! ========== free memory BEGIN ========== call dealloc_surfvec - call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, & - alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) call dealloc_pec - call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0, & - stnext,stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & - pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) ! =========== free memory END =========== + + end block + end associate end subroutine gray_main @@ -1736,7 +1745,7 @@ contains integer :: nlarmor, ithn, ierrcd real(wp_) :: mu, rbn, rbx real(wp_) :: zeff, cst2, bmxi, bmni, fci - real(wp_), dimension(:), pointer :: eccdpar=>null() + real(wp_), dimension(:), allocatable :: eccdpar real(wp_) :: effjcd, effjcdav, Btot complex(wp_) :: e(3) @@ -1831,7 +1840,7 @@ contains ithn, cst2, fjncl, eccdpar, effjcd, iokhawa, ierrcd) end select error = error + ierrcd - if (associated(eccdpar)) deallocate(eccdpar) + if (allocated(eccdpar)) deallocate(eccdpar) ! current drive efficiency [A⋅m/W] effjcdav = rbavi*effjcd diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 6c06a84..ea76ff6 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -400,6 +400,8 @@ contains err = 1 end do + ! set computed values + params%raytracing%nray = 1 + (params%raytracing%nrayr-1) * params%raytracing%nrayth contains function ini_handler(section, name, value) result(err) @@ -493,6 +495,9 @@ contains params%output%ipec = mod(params%output%ipec, 2) close(u) + + ! set computed values + params%raytracing%nray = 1 + (params%raytracing%nrayr-1) * params%raytracing%nrayth end subroutine read_gray_params diff --git a/src/main.f90 b/src/main.f90 index 3160d85..5a03e2b 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -472,7 +472,7 @@ contains use gray_params, only : gray_parameters, print_parameters use beams, only : launchangles2n, xgygcoeff use magsurf_data, only : flux_average, dealloc_surfvec - use beamdata, only : init_btr, dealloc_beam + use beamdata, only : init_btr use pec, only : pec_init, postproc_profiles, dealloc_pec, & rhop_tab, rhot_tab use gray_core, only : print_headers, print_finals, print_pec, & @@ -492,29 +492,12 @@ contains real(wp_) :: chipol, psipol, st real(wp_) :: drhotp, drhotj, dpdvmx, jphimx - real(wp_), dimension(3) :: anv0 - real(wp_), dimension(:, :), pointer :: yw=>null(), ypw=>null(), gri=>null() - real(wp_), dimension(:, :, :), pointer :: xc=>null(), du1=>null(), ggri=>null() - - real(wp_), dimension(:, :), pointer :: psjki=>null(), ppabs=>null(), ccci=>null() - real(wp_), dimension(:), pointer :: tau0=>null(), alphaabs0=>null(), & - dids0=>null(), ccci0=>null() - real(wp_), dimension(:), pointer :: p0jk=>null() - complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null() - integer, dimension(:), pointer :: iiv=>null() - ! ======== set environment BEGIN ======== ! Compute X=(ω_pe/ω)² and Y=ω_ce/ω (with B=1) call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn) - ! Compute the initial cartesian wavevector (anv0) - call launchangles2n(params%antenna, anv0) - ! Initialise the ray variables (beamtracing) - call init_btr(params%raytracing, yw, ypw, xc, du1, & - gri, ggri, psjki, ppabs, ccci, & - tau0, alphaabs0, dids0, ccci0, & - p0jk, ext, eyt, iiv) + call init_btr(params%raytracing) ! Initialise the dispersion module if (params%ecrh_cd%iwarm > 1) call expinit @@ -556,8 +539,6 @@ contains ! Free memory call dealloc_surfvec ! for fluxval - call dealloc_beam(yw, ypw, xc, du1, gri, ggri, psjki, ppabs, ccci, & - tau0, alphaabs0, dids0, ccci0, p0jk, ext, eyt, iiv) call dealloc_pec end subroutine sum_profiles diff --git a/src/multipass.f90 b/src/multipass.f90 index 1bec753..ea0b6c7 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -8,9 +8,6 @@ module multipass implicit none - integer, save :: nbeam_max ! max n of beams active at a time - integer, save :: nbeam_tot ! total n of beams - contains subroutine plasma_in(i, x, N, Bres, sox, cpl, psi, chi, iop, ext, eyt, perfect) @@ -18,15 +15,15 @@ contains use const_and_precisions, only: cm ! subroutine arguments - integer, intent(in) :: i ! ray index - real(wp_), intent(in) :: x(3), N(3) ! position, refactive index - real(wp_), intent(in) :: Bres ! resonant B field - integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X - real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X) - real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles - integer, intent(inout), pointer :: iop(:) ! inside/outside plasma flag - complex(wp_), intent(inout), pointer :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y) - logical, intent(in) :: perfect ! whether to assume perfect coupling + integer, intent(in) :: i ! ray index + real(wp_), intent(in) :: x(3), N(3) ! position, refactive index + real(wp_), intent(in) :: Bres ! resonant B field + integer, intent(in) :: sox ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X + real(wp_), intent(out) :: cpl(2) ! power coupling vector (O, X) + real(wp_), intent(out) :: psi, chi ! polarisation ellipse angles + integer, intent(inout) :: iop(:) ! inside/outside plasma flag + complex(wp_), intent(inout) :: ext(:), eyt(:) ! ray polarisation vector (e_x, e_y) + logical, intent(in) :: perfect ! whether to assume perfect coupling ! local variables real(wp_) :: R, z, cosphi, sinphi, B_phi, B_R, B_z @@ -75,20 +72,23 @@ contains end subroutine plasma_in - subroutine plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! ray exits plasma -! arguments - integer, intent(in) :: i ! ray index - real(wp_), dimension(3), intent(in) :: xv,anv - real(wp_), intent(in) :: bres - integer, intent(in) :: sox - integer, dimension(:), intent(inout), pointer :: iop ! in/out plasma flag - complex(wp_), dimension(:), intent(out), pointer :: ext,eyt -! local variables - real(wp_) :: rm,csphi,snphi,bphi,br,bz - real(wp_), dimension(3) :: bv,xmv -! - iop(i)=iop(i)+1 ! in->out - + subroutine plasma_out(i, xv, anv, bres, sox, iop, ext, eyt) + ! Ray exits plasma + + ! subroutine arguments + integer, intent(in) :: i ! ray index + real(wp_), intent(in) :: xv(3), anv(3) + real(wp_), intent(in) :: bres + integer, intent(in) :: sox + integer, intent(inout) :: iop(:) ! in/out plasma flag + complex(wp_), intent(out) :: ext(:), eyt(:) + + ! local variables + real(wp_) :: rm, csphi, snphi, bphi, br, bz + real(wp_), dimension(3) :: bv, xmv + + iop(i)=iop(i)+1 ! in->out + xmv=xv*0.01_wp_ ! convert from cm to m rm=sqrt(xmv(1)**2+xmv(2)**2) csphi=xmv(1)/rm @@ -97,31 +97,29 @@ contains bv(1)=br*csphi-bphi*snphi bv(2)=br*snphi+bphi*csphi bv(3)=bz - call pol_limit(anv,bv,bres,sox,ext(i),eyt(i)) ! polarization at plasma exit - + call pol_limit(anv,bv,bres,sox,ext(i),eyt(i)) ! polarization at plasma exit end subroutine plasma_out -! ------------------------------ -! subroutine wall_in(i) ! ray enters vessel -! integer, intent(in) :: i ! ray index -! -! iow(i)=iow(i)+1 -! end subroutine wall_in -! ------------------------------ - subroutine wall_out(i,ins,xv,anv,bres,sox,psipol1,chipol1,iow,iop,ext,eyt) ! ray exits vessel -! arguments - integer, intent(in) :: i ! ray index - logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) - real(wp_), dimension(3), intent(inout) :: xv,anv - real(wp_), intent(in) :: bres - integer, intent(in) :: sox - real(wp_), intent(out) :: psipol1,chipol1 - integer, dimension(:), intent(inout), pointer :: iow,iop ! in/out vessel and plasma flags - complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt -! local variables + + + subroutine wall_out(i, ins, xv, anv, bres, sox, psipol1, chipol1, iow, iop, ext, eyt) + ! Ray exits vessel + + ! subroutine arguments + integer, intent(in) :: i ! ray index + logical, intent(in) :: ins ! inside plasma? (ins=1 plasma/wall overlap) + real(wp_), intent(inout) :: xv(3), anv(3) + real(wp_), intent(in) :: bres + integer, intent(in) :: sox + real(wp_), intent(out) :: psipol1,chipol1 + integer, intent(inout) :: iow(:) ! in/out vessel and plasma flags + integer, intent(inout) :: iop(:) ! in/out vessel and plasma flags + complex(wp_), intent(inout) :: ext(:), eyt(:) + + ! local variables integer :: irfl real(wp_), dimension(3) :: xvrfl,anvrfl,walln complex(wp_) :: ext1,eyt1 -! + iow(i)=iow(i)+1 ! out->in if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! plasma-wall overlapping call wall_refl(xv-dst*anv,anv,ext(i),eyt(i),xvrfl,anvrfl,ext1,eyt1,walln,irfl) ! ray reflects at wall @@ -138,23 +136,26 @@ contains end if end subroutine wall_out -! ------------------------------ - subroutine initbeam(i,iroff,iboff,iwait,stv,jphi_beam,pins_beam,currins_beam, & - dpdv_beam,jcd_beam) ! initialization at beam propagation start + + + subroutine initbeam(i, iroff, iboff, iwait, stv, jphi_beam, pins_beam, & + currins_beam, dpdv_beam, jcd_beam) + ! Initialization at beam propagation start + use logger, only : log_info, log_warning -! arguments - integer, intent(in) :: i ! beam index - logical, dimension(:,:), intent(in), pointer :: iroff ! global ray status (F = active, T = inactive) - logical, intent(out) :: iboff - logical, dimension(:), intent(out), pointer :: iwait - real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam, & - currins_beam,dpdv_beam,jcd_beam,stv + ! subroutine arguments + integer, intent(in) :: i ! beam index + logical, intent(in) :: iroff(:,:) ! global ray status (F=active, T=inactive) + logical, intent(out) :: iboff + logical, intent(out) :: iwait(:) + real(wp_), dimension(:), intent(out) :: jphi_beam, pins_beam, currins_beam + real(wp_), dimension(:), intent(out) :: dpdv_beam, jcd_beam, stv character(256) :: msg ! buffer for formatting log messages - iboff = .false. ! beam status (F = active, T = inactive) - iwait = iroff(:,i) ! copy ray status for current beam from global ray status - if(all(iwait)) then ! no rays active => stop beam + iboff = .false. ! beam status (F = active, T = inactive) + iwait = iroff(:,i) ! copy ray status for current beam from global ray status + if(all(iwait)) then ! no rays active => stop beam iboff = .true. else if (any(iwait)) then ! only some rays active @@ -162,131 +163,70 @@ contains call log_warning(msg, mod='multipass', proc='initbeam') end if - stv = zero ! starting step - - jphi_beam = zero ! 1D beam profiles + stv = zero ! starting step + jphi_beam = zero ! 1D beam profiles pins_beam = zero currins_beam = zero dpdv_beam = zero jcd_beam = zero end subroutine initbeam -! ------------------------------ - subroutine initmultipass(i,iox,iroff,yynext,yypnext,yw0,ypw0,stnext,p0ray, & - taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) ! initialization before pass loop -! arguments - logical, intent(in) :: i ! ipol - integer, intent(in) :: iox ! mode active on 1st pass - logical, dimension(:,:), intent(out), pointer :: iroff ! global ray status (F = active, T = inactive) - real(wp_), dimension(:), intent(out), pointer :: p0ray,tau1,etau1,cpl1,lgcpl1, & - psipv,chipv - real(wp_), dimension(:,:), intent(out), pointer :: yw0,ypw0,stnext,taus,cpls - real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext -! - iroff = .false. ! global ray status (F = active, T = inactive) - if(.not. i) call turnoffray(0,1,3-iox,iroff) ! !ipol => stop other mode (iox=1/2 -> stop ib=2/1 at first pass) - yynext = zero ! starting beam coordinates (1) - yypnext = zero ! starting beam coordinates (2) - yw0 = zero ! temporary beam coordinates (1) - ypw0 = zero ! temporary beam coordinates (2) - stnext = zero ! starting beam step - p0ray = zero ! starting beam power - taus = zero ! beam tau from previous passes + + + subroutine initmultipass(i, iox, iroff, yynext, yypnext, yw0, ypw0, stnext, p0ray, & + taus, tau1, etau1, cpls, cpl1, lgcpl1, psipv, chipv) + ! Initialization before pass loop + + ! subroutine arguments + logical, intent(in) :: i ! ipol + integer, intent(in) :: iox ! mode active on 1st pass + logical, dimension(:,:), intent(out) :: iroff ! global ray status (F = active, T = inactive) + real(wp_), dimension(:), intent(out) :: p0ray, tau1, etau1, cpl1 + real(wp_), dimension(:), intent(out) :: lgcpl1, psipv, chipv + real(wp_), dimension(:,:), intent(out) :: yw0, ypw0, stnext, taus, cpls + real(wp_), dimension(:,:,:), intent(out) :: yynext, yypnext + + iroff = .false. ! global ray status (F = active, T = inactive) + if(.not. i) call turnoffray(0,1,3-iox,iroff) ! !ipol => stop other mode (iox=1/2 -> stop ib=2/1 at first pass) + yynext = zero ! starting beam coordinates (1) + yypnext = zero ! starting beam coordinates (2) + yw0 = zero ! temporary beam coordinates (1) + ypw0 = zero ! temporary beam coordinates (2) + stnext = zero ! starting beam step + p0ray = zero ! starting beam power + taus = zero ! beam tau from previous passes tau1 = zero etau1 = one - cpls = one ! beam coupling from previous passes + cpls = one ! beam coupling from previous passes cpl1 = one lgcpl1 = zero - psipv = zero ! psi polarization angle at vacuum-plasma boundary - chipv = zero ! chi polarization angle at vacuum-plasma boundary + psipv = zero ! psi polarization angle at vacuum-plasma boundary + chipv = zero ! chi polarization angle at vacuum-plasma boundary end subroutine initmultipass -! ------------------------------ - subroutine turnoffray(jk,ip,ib,iroff) ! turn off ray propagation -! arguments - integer, intent(in) :: jk, ip, ib ! ray (0=all rays), pass, beam indexes - logical, dimension(:,:), intent(out), pointer :: iroff ! global ray status (F = active, T = inactive) -! local variables + + + subroutine turnoffray(jk, ip, ib, iroff) + ! Turn off ray propagation + + ! subroutine arguments + integer, intent(in) :: jk, ip, ib ! ray (0=all rays), pass, beam indexes + logical, dimension(:,:), intent(out) :: iroff ! global ray status (F = active, T = inactive) + + ! local variables integer :: ipx, i1, i2 -! - if(jk==0) then ! stop all rays - do ipx=ip,ipass ! from pass ip to last pass - i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1 ! first derived beam at pass ipx - i2 = 2**ipx-2+2**(ipx-ip)*ib ! last derived beam at pass ipx (i1=i2 for ipx=ip) + + if(jk==0) then ! stop all rays + do ipx=ip,ipass ! from pass ip to last pass + i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1 ! first derived beam at pass ipx + i2 = 2**ipx-2+2**(ipx-ip)*ib ! last derived beam at pass ipx (i1=i2 for ipx=ip) iroff(:,i1:i2) = .true. end do - else ! only stop ray jk + else ! only stop ray jk do ipx=ip,ipass i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1 i2 = 2**ipx-2+2**(ipx-ip)*ib iroff(jk,i1:i2) = .true. end do end if - end subroutine turnoffray -! ------------------------------ - subroutine alloc_multipass(dim,iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & - stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & - pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) - integer :: dim - logical, dimension(:), intent(out), pointer :: iwait - logical, dimension(:,:), intent(out), pointer :: iroff - integer, dimension(:), intent(out), pointer :: iop,iow - real(wp_), dimension(:), intent(out), pointer :: jphi_beam,pins_beam,currins_beam, & - dpdv_beam,jcd_beam,stv,tau1,etau1,cpl1,lgcpl1,p0ray,psipv,chipv - real(wp_), dimension(:,:), intent(out), pointer :: taus,cpls,stnext,yw0,ypw0 - real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext - - call dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext,stv, & - p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam,pins_beam,currins_beam, & - dpdv_beam,jcd_beam,psipv,chipv) - - nbeam_max = 2**ipass ! max n of beams active at a time - nbeam_tot = 2**(ipass+1)-2 ! total n of beams - - allocate(iwait(nray),iroff(nray,nbeam_tot),iop(nray),iow(nray), & - yynext(6,nray,nbeam_max-2),yypnext(6,nray,nbeam_max-2), & - yw0(6,nray),ypw0(6,nray),stnext(nray,nbeam_tot),stv(nray), & - p0ray(nray),taus(nray,nbeam_tot),tau1(nray),etau1(nray), & - cpls(nray,nbeam_tot),cpl1(nray),lgcpl1(nray),jphi_beam(dim), & - pins_beam(dim),currins_beam(dim),dpdv_beam(dim),jcd_beam(dim), & - psipv(0:nbeam_tot),chipv(0:nbeam_tot)) - - end subroutine alloc_multipass -! ------------------------------ - subroutine dealloc_multipass(iwait,iroff,iop,iow,yynext,yypnext,yw0,ypw0,stnext, & - stv,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,jphi_beam, & - pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv) - logical, dimension(:), intent(out), pointer :: iwait - logical, dimension(:,:), intent(out), pointer :: iroff - integer, dimension(:), intent(out), pointer :: iop,iow - real(wp_), dimension(:), intent(out), pointer :: stv,p0ray,tau1,etau1,cpl1,lgcpl1, & - jphi_beam,pins_beam,currins_beam,dpdv_beam,jcd_beam,psipv,chipv - real(wp_), dimension(:,:), intent(out), pointer :: yw0,ypw0,stnext,taus,cpls - real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext - - if (associated(iwait)) deallocate(iwait) - if (associated(iroff)) deallocate(iroff) - if (associated(iop)) deallocate(iop) - if (associated(iow)) deallocate(iow) - if (associated(yynext)) deallocate(yynext) - if (associated(yypnext)) deallocate(yypnext) - if (associated(yw0)) deallocate(yw0) - if (associated(ypw0)) deallocate(ypw0) - if (associated(stnext)) deallocate(stnext) - if (associated(stv)) deallocate(stv) - if (associated(p0ray)) deallocate(p0ray) - if (associated(taus)) deallocate(taus) - if (associated(tau1)) deallocate(tau1) - if (associated(etau1)) deallocate(etau1) - if (associated(cpls)) deallocate(cpls) - if (associated(cpl1)) deallocate(cpl1) - if (associated(lgcpl1)) deallocate(lgcpl1) - if (associated(jphi_beam)) deallocate(jphi_beam) - if (associated(pins_beam)) deallocate(pins_beam) - if (associated(currins_beam)) deallocate(currins_beam) - if (associated(dpdv_beam)) deallocate(dpdv_beam) - if (associated(jcd_beam)) deallocate(jcd_beam) - if (associated(psipv)) deallocate(psipv) - if (associated(chipv)) deallocate(chipv) - end subroutine dealloc_multipass end module multipass