From f600108884a21e2f2993cda82db6517410d8dc6a Mon Sep 17 00:00:00 2001 From: Daniele Micheletti Date: Thu, 28 Mar 2019 09:50:28 +0000 Subject: [PATCH] trunk/graycore: minor fixes --- src/graycore.f90 | 65 +++++++++++++++++++++++------------------------ src/multipass.f90 | 64 +++++++++++++++++++++++++--------------------- 2 files changed, 67 insertions(+), 62 deletions(-) diff --git a/src/graycore.f90 b/src/graycore.f90 index 7884d95..7a1e0b1 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -185,7 +185,7 @@ contains nbeam_pass=1 ! max n of beam per pass index_rt=0 ! global beam index: 1,O 2,X 1st pass - ! / \ / \ + ! | | | | do ip=1,ipass ! 3,O 4,X 5,O 6,X 2nd pass pabs_pass = zero @@ -194,7 +194,7 @@ contains nbeam_pass = 2*nbeam_pass ! max n of beams in current pass if(ip.gt.1) then - du1 = zero ! igrad=0 from 2nd pass + du1 = zero gri = zero ggri = zero if(ip.eq.ipass) cpl = (/zero,zero/) ! no successive passes @@ -220,12 +220,12 @@ contains call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) call print_headers((/' '/),index_rt) - if(ip.eq.1) then ! first pass - igrad_b = igrad + if(ip.eq.1) then ! 1st pass + igrad_b = igrad ! * input value, igrad_b=0 from 2nd pass - tau1 = zero + tau1 = zero ! * tau from previous passes etau1 = one - cpl1 = one + cpl1 = one ! * coupling from previous passes lgcpl1 = zero p0ray = p0jk ! * initial beam power @@ -233,7 +233,6 @@ contains call set_pol(yw,bres,sox,psipol,chipol,ext,eyt) ! * initial polarization do jk=1,nray - if(iwait(jk)) cycle zzm = yw(3,jk)*0.01_wp_ rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ @@ -244,7 +243,7 @@ contains end if end do - else ! successive passes + else ! 2nd+ passes ipar = (index_rt+1)/2-1 ! * parent beam index yw = yynext(:,:,ipar) ! * starting coordinates from ypw = yypnext(:,:,ipar) ! parent beam last step @@ -307,7 +306,7 @@ contains if(ent_pl) then ! ray enters plasma call plasma_in(jk,xv,anv,bres,sox,cpl,psipol,chipol,iop,ext,eyt) - if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray never entered in plasma) => continue current pass + if(iop(jk).eq.1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass if(ipol.eq.0) then ! + IF single mode propagation cpl = cpl0 @@ -333,29 +332,29 @@ contains igrad_b = 0 ! + switch to ray-tracing iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray - if(ip.lt.ipass) then ! not last pass - yynext(:,jk,index_rt) = yw0(:,jk) ! + copy starting coordinates - yypnext(:,jk,index_rt) = ypw0(:,jk) ! for next pass from last step - stnext(jk,index_rt) = stv(jk) - dst ! + starting step for next pass = last step + if(ip.lt.ipass) then ! + not last pass + yynext(:,jk,index_rt) = yw0(:,jk) ! . copy starting coordinates + yypnext(:,jk,index_rt) = ypw0(:,jk) ! for next pass from last step + stnext(jk,index_rt) = stv(jk) - dst ! . starting step for next pass = last step - if(cpl(1).lt.etaucr) then ! + low coupled power for O-mode => de-activate derived rays + if(cpl(1).lt.etaucr) then ! . low coupled power for O-mode => de-activate derived rays call turnoffray(jk,ip+1,2*ib-1,iroff) if(cpl(1).le.comp_tiny) cpl(1)=etaucr end if - if(cpl(2).lt.etaucr) then ! + low coupled power for X-mode => de-activate derived rays + if(cpl(2).lt.etaucr) then ! . low coupled power for X-mode => de-activate derived rays call turnoffray(jk,ip+1,2*ib,iroff) if(cpl(2).le.comp_tiny) cpl(2)=etaucr end if - taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! + starting tau for next O-mode pass - cpls(jk,iO) = cpl1(jk) * cpl(1) ! + cumulative coupling for next O-mode pass - cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! + cumulative coupling for next X-mode pass + taus(jk,iO:iO+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass + cpls(jk,iO) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass + cpls(jk,iO+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass - if(jk.eq.1) then ! + polarization angles at plasma boundary for central ray + if(jk.eq.1) then ! . polarization angles at plasma boundary for central ray psipv(iO:iO+1) = psipol chipv(iO:iO+1) = chipol end if - else + else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass cpl = (/zero,zero/) end if end if @@ -365,7 +364,7 @@ contains end if if(ent_wl) then ! ray enters vessel - iow(jk)=iow(jk)+1 + iow(jk)=iow(jk)+1 ! * out->in else if(ext_wl) then ! ray exit vessel call wall_out(jk,ins_pl,xv,anv,bres,sox,psipol,chipol,iow,iop,ext,eyt) @@ -385,7 +384,7 @@ contains chipv(index_rt) = chipol end if - if(ins_pl) then ! * plasma-wall overlapping => wall crossing=plasma crossing => end of current pass + if(ins_pl) then ! * plasma-wall overlapping => wall+plasma crossing => end of current pass iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray if(ip.lt.ipass) then ! + not last pass @@ -418,11 +417,11 @@ contains iopmin = min(iopmin,iop(jk)) if(ip.lt.ipass) then ! not last pass - yw0(:,jk) = yw(:,jk) ! * store current coordinates - ypw0(:,jk) = ypw(:,jk) ! in case pass ends on next step + yw0(:,jk) = yw(:,jk) ! * store current coordinates in case + ypw0(:,jk) = ypw(:,jk) ! current pass ends on next step end if - ! compute ECRH&CD + ! compute ECRH&CD if (inside plasma & power available>0 & ray still active) if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. & (tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check tekev=temp(psinv) @@ -456,14 +455,14 @@ contains dids0(jk)=dids ccci0(jk)=ccci(jk,i) - if(iwait(jk)) then + if(iwait(jk)) then ! copy values from last pass for inactive ray ppabs(jk,i:nstep) = ppabs(jk,i-1) ccci(jk,i:nstep) = ccci(jk,i-1) psjki(jk,i:nstep) = psjki(jk,i-1) else call print_output(i,jk,stv(jk),p0ray(jk)/etau1(jk),xv,psinv, & btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, & - nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) !****** p0ray(jk)/etau1(jk) fattore normalizzazione per dids da controllare + nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) end if end do @@ -490,7 +489,7 @@ contains istop_pass = istop_pass +1 ! * +1 non propagating beam if(ip.lt.ipass) call turnoffray(0,ip,ib,iroff) ! * de-activate derived beams exit - else if(all(iwait)) then ! all rays in current beam are waiting for next pass + else if(all(iwait)) then ! all rays in current beam are waiting for next pass => do not increase istop_pass exit end if end do @@ -514,10 +513,10 @@ contains icd_pass(iox) = icd_pass(iox) + icd_beam if(ip.lt.ipass) then - cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1) / sum(p0ray * exp(-tau0)) ! average O-mode coupling - cpl_beam2 = one - cpl_beam1 - cpl_cbeam1 = cpls(1,iO)/cpl1(1) - cpl_cbeam2 = one - cpl_cbeam1 + cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1) / sum(p0ray * exp(-tau0)) ! average O-mode coupling for next beam + cpl_beam2 = one - cpl_beam1 ! average X-mode coupling for next beam + cpl_cbeam1 = cpls(1,iO)/cpl1(1) ! central ray O-mode coupling for next beam + cpl_cbeam2 = one - cpl_cbeam1 ! central ray X-mode coupling for next beam else cpl_beam1 = zero cpl_beam2 = zero @@ -533,7 +532,7 @@ contains write(*,'(a,f9.4)') 'I_tot (kA) = ',icd_beam*1.0e3_wp_ if(ip.lt.ipass) then write(*,'(a,2(f9.4,1x))') 'Coupling (average, O/X):',cpl_beam1,cpl_beam2 ! average coupling for next O/X beams - write(*,'(a,2(f9.4,1x))') 'Coupling (c. beam, O/X):',cpl_cbeam1,cpl_cbeam2 ! central beam coupling for next O/X beams + write(*,'(a,2(f9.4,1x))') 'Coupling (ctr ray, O/X):',cpl_cbeam1,cpl_cbeam2 ! central ray coupling for next O/X beams end if call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, & diff --git a/src/multipass.f90 b/src/multipass.f90 index e0393a9..e0bd4f7 100755 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -8,28 +8,29 @@ module multipass implicit none - integer, save :: nbeam_max - integer, save :: nbeam_tot + 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,xv,anv,bres,sox,cpl,psipol1,chipol1,iop,ext,eyt) ! ray enters plasma implicit none +! arguments integer, intent(in) :: i ! ray index real(wp_), dimension(3), intent(in) :: xv,anv real(wp_), intent(in) :: bres,sox real(wp_), dimension(2), intent(out) :: cpl ! coupling (O/X) real(wp_), intent(out) :: psipol1,chipol1 - integer, dimension(:), intent(inout), pointer :: iop + integer, dimension(:), intent(inout), pointer :: iop ! in/out plasma flag complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt -! +! local variables real(wp_) :: rm,csphi,snphi,bphi,br,bz real(wp_) :: qq1,uu1,vv1,qq,uu,vv,powcpl real(wp_), dimension(3) :: bv,xmv complex(wp_) :: ext1,eyt1 ! - iop(i)=iop(i)+1 + iop(i)=iop(i)+1 ! out->in xmv=xv*0.01_wp_ ! convert from cm to m rm=sqrt(xmv(1)**2+xmv(2)**2) @@ -43,7 +44,7 @@ contains call pol_limit(anv,bv,bres,sox,ext1,eyt1) call stokes_ce(ext1,eyt1,qq1,uu1,vv1) ! stokes parameter at plasma entrance call stokes_ce(ext(i),eyt(i),qq,uu,vv) ! stokes parameter at plasma exit/wall reflection - powcpl = half*(one + vv*vv1+uu*uu1+qq*qq1) + powcpl = half*(one + vv*vv1+uu*uu1+qq*qq1) ! coupling for incoming mode if(sox.eq.-one) then ! incoming mode = O cpl=(/powcpl,one-powcpl/) @@ -51,7 +52,7 @@ contains cpl=(/one-powcpl,powcpl/) end if - if(i.eq.1) then + if(i.eq.1) then ! polarization angles at plasma entrace for central ray call polellipse(qq1,uu1,vv1,psipol1,chipol1) psipol1=psipol1/degree ! convert from rad to degree chipol1=chipol1/degree @@ -63,16 +64,17 @@ contains ! ------------------------------ subroutine plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) ! ray exits plasma implicit none +! arguments integer, intent(in) :: i ! ray index real(wp_), dimension(3), intent(in) :: xv,anv real(wp_), intent(in) :: bres,sox - integer, dimension(:), intent(inout), pointer :: iop + 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 + 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) @@ -95,23 +97,24 @@ contains ! ------------------------------ subroutine wall_out(i,ins,xv,anv,bres,sox,psipol1,chipol1,iow,iop,ext,eyt) ! ray exits vessel implicit none +! 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,sox real(wp_), intent(out) :: psipol1,chipol1 - integer, dimension(:), intent(inout), pointer :: iow,iop + integer, dimension(:), intent(inout), pointer :: iow,iop ! in/out vessel and plasma flags complex(wp_), dimension(:), intent(inout), pointer :: ext,eyt - ! +! local variables integer :: irfl real(wp_) :: qq1,uu1,vv1 real(wp_), dimension(3) :: xvrfl,anvrfl,walln complex(wp_) :: ext1,eyt1 - ! - iow(i)=iow(i)+1 - if(ins) call plasma_out(i,xv,anv,bres,sox,iop,ext,eyt) - call wall_refl(xv-dst*anv,anv,ext(i),eyt(i),xvrfl,anvrfl,ext1,eyt1,walln,irfl) - ext(i) = ext1 +! + 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 + ext(i) = ext1 ! save parameters at wall reflection eyt(i) = eyt1 xv = xvrfl anv = anvrfl @@ -131,8 +134,9 @@ contains subroutine initbeam(i,iroff,iboff,iwait,stv,jphi_beam,pins_beam,currins_beam, & dpdv_beam,jcd_beam) ! initialization at beam propagation start implicit none +! arguments integer, intent(in) :: i ! beam index - logical, dimension(:,:), intent(in), pointer :: iroff + 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, & @@ -144,7 +148,7 @@ contains iboff = .true. write(*,*) write(*,'("Beam ",i5," inactive")'),i - else if(.not.all(.not.iwait)) then ! no all rays active + else if(.not.all(.not.iwait)) then ! only some rays active write(*,*) write(*,'("WARNING: not all rays in beam ",i5," are active")'),i end if @@ -161,25 +165,26 @@ contains subroutine initmultipass(i,iox,iroff,yynext,yypnext,yw0,ypw0,stnext,p0ray, & taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) ! initialization before pass loop implicit none +! arguments integer, intent(in) :: i, iox ! ipol, mode active on 1st pass - logical, dimension(:,:), intent(out), pointer :: iroff + 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(i.eq.0) call turnoffray(0,1,3-iox,iroff) ! ipol = 0 => stop other mode (iox=1/2 -> stop ib=2/1) + if(i.eq.0) call turnoffray(0,1,3-iox,iroff) ! ipol = 0 => 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 ! starting beam tau - tau1 = zero ! total beam tau + taus = zero ! beam tau from previous passes + tau1 = zero etau1 = one - cpls = one + cpls = one ! beam coupling from previous passes cpl1 = one lgcpl1 = zero psipv = zero ! psi polarization angle at vacuum-plasma boundary @@ -188,15 +193,16 @@ contains ! ------------------------------ subroutine turnoffray(jk,ip,ib,iroff) ! turn off ray propagation implicit none +! arguments integer, intent(in) :: jk, ip, ib ! ray (0=all rays), pass, beam indexes - logical, dimension(:,:), intent(out), pointer :: iroff -! + logical, dimension(:,:), intent(out), pointer :: 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 - i1 = 2**ipx-2+2**(ipx-ip)*(ib-1)+1 - i2 = 2**ipx-2+2**(ipx-ip)*ib + 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