diff --git a/Makefile b/Makefile index 929579a..9723231 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,7 @@ EXE=gray # Objects list MAINOBJ=main.o OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \ - dierckx.o dispersion.o eccd.o eierf.o graycore.o gray-externals.o \ + dierckx.o dispersion.o eccd.o eierf.o errcodes.o graycore.o gray-externals.o \ gray_params.o equilibrium.o magsurf_data.o math.o minpack.o numint.o \ pec.o polarization.o quadpack.o reflections.o simplespline.o utils.o @@ -29,7 +29,7 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) main.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \ graycore.o gray_params.o reflections.o graycore.o: const_and_precisions.o beamdata.o beams.o coreprofiles.o \ - dispersion.o equilibrium.o gray-externals.o gray_params.o \ + dispersion.o equilibrium.o errcodes.o gray-externals.o gray_params.o \ pec.o polarization.o reflections.o utils.o gray-externals.o: const_and_precisions.o beams.o coreprofiles.o dierckx.o \ dispersion.o eccd.o gray_params.o \ @@ -41,9 +41,11 @@ conical.o: const_and_precisions.o coreprofiles.o: const_and_precisions.o dierckx.o gray_params.o simplespline.o \ utils.o dierckx.o: const_and_precisions.o -dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o -eccd.o: const_and_precisions.o conical.o magsurf_data.o dierckx.o numint.o +dispersion.o: const_and_precisions.o eierf.o errcodes.o math.o quadpack.o +eccd.o: const_and_precisions.o conical.o dierckx.o errcodes.o magsurf_data.o \ + numint.o eierf.o: const_and_precisions.o +errcodes.o: const_and_precisions.o gray_params.o: const_and_precisions.o utils.o equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o \ utils.o gray_params.o diff --git a/src/beamdata.f90 b/src/beamdata.f90 index c4573fa..4357276 100644 --- a/src/beamdata.f90 +++ b/src/beamdata.f90 @@ -8,16 +8,17 @@ module beamdata contains - subroutine init_rtr(rtrparam,ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) + subroutine init_rtr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) use gray_params, only : rtrparam_type use const_and_precisions, only : zero,half,two implicit none type(rtrparam_type), intent(in) :: rtrparam real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & - gri,psjki,tauv,alphav,ppabs,dids,ccci + gri,psjki,ppabs,ccci real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri - real(wp_), dimension(:), intent(out), allocatable :: p0jk + real(wp_), dimension(:), intent(out), allocatable :: tau0,alphaabs0, & + dids0,ccci0,p0jk complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt integer, dimension(:), intent(out), allocatable :: iiv @@ -40,8 +41,8 @@ contains nstep=rtrparam%nstep - call alloc_beam(ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) + call alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) end subroutine init_rtr function rayi2jk(i) result(jk) @@ -100,53 +101,56 @@ contains end if end function rayjk2i - subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) + subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) implicit none real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & - gri,psjki,tauv,alphav,ppabs,dids,ccci + gri,psjki,ppabs,ccci real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri - real(wp_), dimension(:), intent(out), allocatable :: p0jk + real(wp_), dimension(:), intent(out), allocatable :: tau0,alphaabs0, & + dids0,ccci0,p0jk complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt integer, dimension(:), intent(out), allocatable :: iiv - call dealloc_beam(ywork,ypwork,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,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), tauv(nray,nstep), alphav(nray,nstep), & - ppabs(nray,nstep), dids(nray,nstep), ccci(nray,nstep), & - p0jk(nray), ext(nray), eyt(nray), iiv(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,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) + subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) implicit none real(wp_), dimension(:,:), intent(out), allocatable :: ywork,ypwork, & - gri,psjki,tauv,alphav,ppabs,dids,ccci + gri,psjki,ppabs,ccci real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri - real(wp_), dimension(:), intent(out), allocatable :: p0jk + real(wp_), dimension(:), intent(out), allocatable :: tau0,alphaabs0, & + dids0,ccci0,p0jk complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt integer, dimension(:), intent(out), allocatable :: iiv - if (allocated(ywork)) deallocate(ywork) - if (allocated(ypwork)) deallocate(ypwork) - if (allocated(xc)) deallocate(xc) - if (allocated(du1)) deallocate(du1) - if (allocated(gri)) deallocate(gri) - if (allocated(ggri)) deallocate(ggri) - if (allocated(psjki)) deallocate(psjki) - if (allocated(tauv)) deallocate(tauv) - if (allocated(alphav)) deallocate(alphav) - if (allocated(ppabs)) deallocate(ppabs) - if (allocated(dids)) deallocate(dids) - if (allocated(ccci)) deallocate(ccci) - if (allocated(p0jk)) deallocate(p0jk) - if (allocated(ext)) deallocate(ext) - if (allocated(eyt)) deallocate(eyt) - if (allocated(iiv)) deallocate(iiv) + if (allocated(ywork)) deallocate(ywork) + if (allocated(ypwork)) deallocate(ypwork) + if (allocated(xc)) deallocate(xc) + if (allocated(du1)) deallocate(du1) + if (allocated(gri)) deallocate(gri) + if (allocated(ggri)) deallocate(ggri) + if (allocated(psjki)) deallocate(psjki) + if (allocated(ppabs)) deallocate(ppabs) + if (allocated(ccci)) deallocate(ccci) + if (allocated(tau0)) deallocate(tau0) + if (allocated(alphaabs0)) deallocate(alphaabs0) + if (allocated(dids0)) deallocate(dids0) + if (allocated(ccci0)) deallocate(ccci0) + if (allocated(p0jk)) deallocate(p0jk) + if (allocated(ext)) deallocate(ext) + if (allocated(eyt)) deallocate(eyt) + if (allocated(iiv)) deallocate(iiv) end subroutine dealloc_beam subroutine pweight(p0,p0jk) diff --git a/src/dispersion.f90 b/src/dispersion.f90 index f1010d3..e4d87d1 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -183,36 +183,25 @@ subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez) if(i.gt.imxx.and.imxx.gt.1) then if (imx.lt.0) then - write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence & - &disabled.',e12.5)") xg,yg,sqrt(abs(npr2)),npl + err=1 imxx=1 else - write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence & - &failed.',e12.5)") xg,yg,sqrt(abs(npr2)),npl + err=2 exit end if else exit end if - print*,yg,imx,imxx end do ! ! if(i.gt.imx) print*,' i>imx ',yg,errnpr,i ! - if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).eq.huge(one).or. & + if(dble(sqrt(npr2)).lt.zero.or.npr2.ne.npr2.or.abs(npr2).ge.huge(one).or. & abs(npr2).le.tiny(one)) then - write(*,"(' X =',f7.4,' Y =',f7.4,' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(npr2)) npr2=czero - err=99 + err=err+4 end if -! if(dble(npr2).lt.zero) then -! npr2=zero -! print*,' Y =',yg,' npr2 < 0' -! err=99 -! end if -! -! write(11,99) yg,dble(npr2),dimag(npr2),nprf**2,dble(i) ! npr=sqrt(npr2) nprr=dble(npr) diff --git a/src/eccd.f90 b/src/eccd.f90 index ed71c1c..309e944 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -168,6 +168,7 @@ contains ithn,cst2,fcur,eccdpar,effjcd,iokhawa,ierr) use const_and_precisions, only : pi,qesi=>e_,mesi=>me_, & vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_ + use errcodes, only : pcdfp,pcdfc use quadpack, only : dqagsmv implicit none ! local constants @@ -247,7 +248,7 @@ contains epp,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) then - ierr=90 + ierr=ibset(ierr,pcdfp) return end if @@ -291,7 +292,7 @@ contains if (abs(resji).lt.1.0e-10_wp_) then resji=0.0_wp_ else - ierr=91+iokhawa+i + ierr=ibset(ierr,pcdfc+iokhawa+i) return end if end if diff --git a/src/errcodes.f90 b/src/errcodes.f90 new file mode 100644 index 0000000..f3a2d94 --- /dev/null +++ b/src/errcodes.f90 @@ -0,0 +1,80 @@ +module errcodes + implicit none + integer, parameter :: pnpl = 0, lnpl = 2 ! N// too large (2 thresholds) + integer, parameter :: pconv = pnpl + lnpl, lconv = 2 ! Disp. convergence (disabled/failed) + integer, parameter :: pnprre = pconv + lconv, lnprre = 1 ! Re(Nperp)<0 + integer, parameter :: palph = pnprre+ lnprre, lalph = 1 ! alpha<0 + integer, parameter :: pcdfp = palph + lalph, lcdfp = 1 ! fpp integration + integer, parameter :: pcdfc = pcdfp + lcdfp, lcdfc = 3 ! fcur integration (no trap/trap 1/trap 2) +contains + + subroutine check_err(ierr,istop) + implicit none +! arguments + integer, intent(in) :: ierr + integer, intent(out) :: istop + + if(ibits(ierr,pnpl, lnpl )==2 .or. & ! N// too large + ibits(ierr,palph,lalph)==1) then ! alpha < 0 + istop = 1 + else + istop = 0 + end if + end subroutine check_err + + + subroutine print_errn(ierr,i,anpl) + use const_and_precisions, only : wp_ + implicit none +! arguments + integer, intent(in) :: ierr,i + real(wp_), intent(in) :: anpl +! local variables + integer :: ierrs + + ierrs = ibits(ierr,pnpl,lnpl) + if(ierrs/=0) print*,i,' IERR = ', ierrs*2**pnpl,' N// = ',anpl + end subroutine print_errn + + + subroutine print_errhcd(ierr,i,anprre,anprim,alpha) + use const_and_precisions, only : wp_ + implicit none +! arguments + integer, intent(in) :: ierr,i + real(wp_), intent(in) :: anprre,anprim,alpha +! local variables + integer :: ierrs + + ierrs=ibits(ierr,pconv,lconv) + if(ierrs==1) then + print*,i,' IERR = ', ierrs*2**pconv,' Nwarm = ',anprre,anprim, & + ': convergence disabled.' + else if(ierrs==2) then + print*,i,' IERR = ', ierrs*2**pconv,' Nwarm = ',anprre,anprim, & + ': convergence failed.' + end if + + ierrs=ibits(ierr,pnprre,lnprre) + if(ierrs/=0) & + print*,i,' IERR = ', ierrs*2**pconv,' Nwarm = ',anprre,anprim, & + ': Re(Nwarm)<0 or Nwarm**2 invalid.' + + ierrs=ibits(ierr,palph,lalph) + if(ierrs/=0) & + print*,i,' IERR = ', ierrs*2**palph,' alpha = ',alpha + + ierrs=ibits(ierr,pcdfp,lcdfp) + if(ierrs/=0) & + print*,i,' IERR = ', ierrs*2**pcdfp,' fpp integration error' + + ierrs=ibits(ierr,pcdfc,lcdfc) + if(ibits(ierrs,0,1)/=0) & + print*,i,' IERR = ', ierrs*2**pcdfc,' fcur integration error (no trapping)' + if(ibits(ierrs,1,1)/=0) & + print*,i,' IERR = ', ierrs*2**pcdfc,' fcur integration error (1st trapping region)' + if(ibits(ierrs,2,1)/=0) & + print*,i,' IERR = ', ierrs*2**pcdfc,' fcur integration error (2nd trapping region)' + end subroutine print_errhcd + +end module errcodes \ No newline at end of file diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 05e1cf5..da34d94 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -47,13 +47,15 @@ module gray_params integer, save :: ipec,nnd contains - - subroutine read_inputs(filenm,antctrl,eqparam,rwall,prfparam,outparam,unit) - use const_and_precisions, only : wp_ + + subroutine read_params(filenm,rtrparam,hcdparam,antctrl,eqparam,rwall, & + prfparam,outparam,unit) use utils, only : get_free_unit implicit none ! arguments character(len=*), intent(in) :: filenm + type(rtrparam_type), intent(out) :: rtrparam + type(hcdparam_type), intent(out) :: hcdparam type(antctrl_type), intent(out) :: antctrl type(eqparam_type), intent(out) :: eqparam real(wp_), intent(out) :: rwall @@ -70,77 +72,7 @@ contains end if open(u,file=filenm,status= 'old',action='read') -! alpha0, beta0 (cartesian) launching angles - read(u,*) antctrl%alpha, antctrl%beta -! p0mw injected power (MW) - read(u,*) antctrl%power -! abs(iox)=1/2 OM/XM -! psipol0,chipol0 polarization angles at the antenna (if iox<0) - read(u,*) antctrl%iox, antctrl%psi, antctrl%chi - -! ibeam=0 :read data for beam as above -! ibeam=1 :read data from file simple astigmatic beam -! ibeam=2 :read data from file general astigmatic beam - read(u,*) antctrl%ibeam - read(u,*) antctrl%filenm - -! iequil=0 :vacuum -! iequil=1 :analytical equilibrium -! iequil=2 :read eqdsk - read(u,*) eqparam%iequil - read(u,*) eqparam%filenm -! icocos :index for equilibrium from COCOS - O. Sauter Feb 2012 -! ipsinorm :0 standard EQDSK format, 1 format Portone summer 2004 - read(u,*) eqparam%icocos, eqparam%ipsinorm, eqparam%idesc, eqparam%ifreefmt -! ixp=0,-1,+1 : no X point , bottom/up X point -! ssplps : spline parameter for psi interpolation - read(u,*) eqparam%ixp, eqparam%ssplps !, eqparam%ssplf - eqparam%ssplf=0.01_wp_ -! signum of toroidal B and I -! factb factor for magnetic field (only for numerical equil) -! scaling adopted: beta=const, qpsi=const, nustar=const - read(u,*) eqparam%sgnb, eqparam%sgni, eqparam%factb - - read(u,*) rwall - -! iprof=0 :analytical density and temp. profiles -! iprof>0 :numerical density and temp. profiles - read(u,*) prfparam%iprof, prfparam%irho ! irho=0,1,2 -> num profiles vs rhot,rhop,psin - read(u,*) prfparam%filenm -! psbnd value of psi ( > 1 ) of density boundary - read(u,*) prfparam%psnbnd !, prfparam%sspld - prfparam%sspld=0.001_wp_ -! iscal ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling -! factT factn factor for Te&ne scaling - read(u,*) prfparam%factte, prfparam%factne, prfparam%iscal - -! ipec=0/1 :pec profiles grid in psi/rhop -! nrho :number of grid steps for pec profiles +1 - read(u,*) outparam%ipec, outparam%nrho -! istpr0 projection step = dsdt*istprj -! istpl0 plot step = dsdt*istpl - read(u,*) outparam%istpr, outparam%istpl - close(u) - end subroutine read_inputs - - subroutine read_params(filenm,rtrparam,hcdparam,unit) - use utils, only : get_free_unit - implicit none -! arguments - character(len=*), intent(in) :: filenm - type(rtrparam_type), intent(out) :: rtrparam - type(hcdparam_type), intent(out) :: hcdparam - integer, intent(in), optional :: unit -! local variables - integer :: u - - if (present(unit)) then - u=unit - else - u = get_free_unit() - end if - open(u,file=filenm,status= 'old',action='read') - +! ========================================================================== ! nrayr number of rays in radial direction ! nrayth number of rays in angular direction ! rwmax normalized maximum radius of beam power @@ -157,7 +89,7 @@ contains ! nstep maximum number of integration steps ! idst=0/1/2 0 integration in s, 1 integr. in ct, 2 integr. in Sr read(u,*) rtrparam%dst, rtrparam%nstep, rtrparam%idst - +! ========================================================================== ! iwarm=0 :no absorption and cd ! iwarm=1 :weakly relativistic absorption ! iwarm=2 :relativistic absorption, n<1 asymptotic expansion @@ -169,6 +101,57 @@ contains ! ieccd 0/1 NO/YES ECCD calculation ieccd>0 different CD models read(u,*) hcdparam%ieccd +! ========================================================================== +! alpha0, beta0 (cartesian) launching angles + read(u,*) antctrl%alpha, antctrl%beta +! p0mw injected power (MW) + read(u,*) antctrl%power +! abs(iox)=1/2 OM/XM +! psipol0,chipol0 polarization angles at the antenna (if iox<0) + read(u,*) antctrl%iox, antctrl%psi, antctrl%chi + +! ibeam=0 :read data for beam as above +! ibeam=1 :read data from file simple astigmatic beam +! ibeam=2 :read data from file general astigmatic beam + read(u,*) antctrl%ibeam + read(u,*) antctrl%filenm +! ========================================================================== +! iequil=0 :vacuum +! iequil=1 :analytical equilibrium +! iequil=2 :read eqdsk + read(u,*) eqparam%iequil + read(u,*) eqparam%filenm +! icocos :index for equilibrium from COCOS - O. Sauter Feb 2012 +! ipsinorm :0 standard EQDSK format, 1 format Portone summer 2004 + read(u,*) eqparam%icocos, eqparam%ipsinorm, eqparam%idesc, eqparam%ifreefmt +! ixp=0,-1,+1 : no X point , bottom/up X point +! ssplps : spline parameter for psi interpolation + read(u,*) eqparam%ixp, eqparam%ssplps !, eqparam%ssplf + eqparam%ssplf=0.01_wp_ +! signum of toroidal B and I +! factb factor for magnetic field (only for numerical equil) +! scaling adopted: beta=const, qpsi=const, nustar=const + read(u,*) eqparam%sgnb, eqparam%sgni, eqparam%factb +! ========================================================================== + read(u,*) rwall +! ========================================================================== +! iprof=0 :analytical density and temp. profiles +! iprof>0 :numerical density and temp. profiles + read(u,*) prfparam%iprof, prfparam%irho ! irho=0,1,2 -> num profiles vs rhot,rhop,psin + read(u,*) prfparam%filenm +! psbnd value of psi ( > 1 ) of density boundary + read(u,*) prfparam%psnbnd !, prfparam%sspld + prfparam%sspld=0.001_wp_ +! iscal ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling +! factT factn factor for Te&ne scaling + read(u,*) prfparam%factte, prfparam%factne, prfparam%iscal +! ========================================================================== +! ipec=0/1 :pec profiles grid in psi/rhop +! nrho :number of grid steps for pec profiles +1 + read(u,*) outparam%ipec, outparam%nrho +! istpr0 projection step = dsdt*istprj +! istpl0 plot step = dsdt*istpl + read(u,*) outparam%istpr, outparam%istpl close(u) end subroutine read_params diff --git a/src/graycore.f90 b/src/graycore.f90 index 7d3c460..0ea12f8 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -18,6 +18,7 @@ contains use beamdata, only : pweight, print_projxyzt, rayi2jk use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl, & zbinf, zbsup + use errcodes, only : check_err, print_errn, print_errhcd use magsurf_data, only : flux_average use beamdata, only : init_rtr, dealloc_beam, nray, nstep, dst use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & @@ -47,23 +48,22 @@ contains integer, intent(out) :: ierr ! local variables - real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ real(wp_), parameter :: taucr = 12._wp_ real(wp_), dimension(:), allocatable :: rhotn real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre - real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,zeff,dersdst,derdnm,st,st0 - real(wp_) :: tau0,alphaabs0,dids0,ccci0 - real(wp_) :: tau,pow,ddr,ddi,taumn,taumx + real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0 + real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava real(wp_), dimension(3) :: xv,anv0,anv real(wp_), dimension(:,:), allocatable :: yw,ypw,gri real(wp_), dimension(:,:,:), allocatable :: xc,du1,ggri - integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,index_rt=1 + integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1 logical :: ins_pl, somein, allout - real(wp_), dimension(:,:), allocatable :: psjki,tauv,alphav,ppabs,dids,ccci + real(wp_), dimension(:,:), allocatable :: psjki,ppabs,ccci + real(wp_), dimension(:), allocatable :: tau0,alphaabs0,dids0,ccci0 real(wp_), dimension(:), allocatable :: p0jk complex(wp_), dimension(:), allocatable :: ext, eyt integer, dimension(:), allocatable :: iiv @@ -103,8 +103,8 @@ contains call xgygcoeff(fghz,ak0,bres,xgcn) call launchangles2n(alpha0,beta0,xv0,anv0) - call init_rtr(rtrp,yw,ypw,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) + call init_rtr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) if(iwarm > 1) call expinit @@ -127,7 +127,7 @@ contains iox=iox0 sox=-1.0_wp_ if(iox==2) sox=1.0_wp_ - call vectinit(psjki,tauv,alphav,ppabs,dids,ccci,iiv) + call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) call ic_gb(xv0,anv0,ak0,w1,w2,ri1,ri2,phiw,phir,yw,ypw,xc,du1,gri,ggri) psipol=psipol0 @@ -138,9 +138,7 @@ contains st0 = zero if(nray>1) call print_projxyzt(st0,yw,0) ! iproj=0 ==> nfilp=8 -! test if at least part of the beam has entered the plsama - somein = .false. - istop = 0 + somein = .false. ! becomes true if at least part of the beam enters the plasma ! beam/ray propagation do i=1,nstep @@ -153,53 +151,38 @@ contains ! update position and grad if(igrad == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri) -! test if the beam is completely out of the plsama - allout = .true. + allout = .true. ! becomes false if at least part of the beam is inside the plsama + ierr = 0 do jk=1,nray ! compute derivatives with updated gradient and local plasma values xv = yw(1:3,jk) anv = yw(4:6,jk) call ywppla_upd(xv,anv,gri(:,jk),ggri(:,:,jk),sox,bres,xgcn,ypw(:,jk), & - psinv,dens,btot,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm) - - if( abs(anpl) > anplth1) then - if(abs(anpl) <= anplth2) then - ierr=97 -! igrad=0 - else - ierr=98 - istop=1 - end if - else - ierr=0 + psinv,dens,btot,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn) +! update global error code and print message + if (ierrn/=0) then + ierr = ior(ierr,ierrn) + call print_errn(ierrn,i,anpl) end if - if(i==1) then - tau0=zero - alphaabs0=zero - dids0=zero - ccci0=zero - else - tau0=tauv(jk,i-1) - alphaabs0=alphav(jk,i-1) - dids0=dids(jk,i-1) - ccci0=ccci(jk,i-1) - end if zzm = xv(3)*0.01_wp_ ins_pl = (psinv>=zero .and. psinv=zbinf .and. zzm<=zbsup) +! test if the beam is completely out of the plsama allout = allout .and. .not.ins_pl +! test if at least part of the beam has entered the plsama somein = somein .or. ins_pl ! compute ECRH&CD - if(ierr==0 .and. iwarm>0 .and. ins_pl .and. tau0<=taucr) then -! print*,i,jk,rayi2jk(jk),psinv,zzm,anpl + if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. tau0(jk)<=taucr) then tekev=temp(psinv) - if (ieccd> 0) zeff=fzeff(psinv) - call alpha_effj(psinv,xg,yg,dens,tekev,zeff,ak0,bres,derdnm, & - anpl,anpr,sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierr) + call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & + sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd) + if (ierrhcd/=0) then + ierr = ior(ierr,ierrhcd) + call print_errhcd(ierrhcd,i,anprre,anprim,alpha) + end if else tekev=zero - zeff=zero alpha=zero didp=zero anprim=zero @@ -210,39 +193,22 @@ contains end if if(nharm>0) iiv(jk)=i -! full storage required only for psjki,ppabs,ccci -! (jk,i) indexing can be removed from tauv,alphav,dids -! adding (jk) indexing to alphaabs0,tau0,dids0,ccci0 psjki(jk,i) = psinv ! computation of optical depth tau, dP/ds, P(s), dI/ds, I(s) - tau=tau0+0.5_wp_*(alpha+alphaabs0)*dersdst*dst - tauv(jk,i)=tau - alphav(jk,i)=alpha + tau=tau0(jk)+0.5_wp_*(alphaabs0(jk)+alpha)*dersdst*dst pow=p0jk(jk)*exp(-tau) !*exp(-tau1v(jk)) ppabs(jk,i)=p0jk(jk)-pow - dids(jk,i)=didp*pow*alpha - ccci(jk,i)=ccci0+0.5_wp_*(dids0+dids(jk,i))*dersdst*dst + dids=didp*pow*alpha + ccci(jk,i)=ccci0(jk)+0.5_wp_*(dids0(jk)+dids)*dersdst*dst + + tau0(jk)=tau + alphaabs0(jk)=alpha + dids0(jk)=dids + ccci0(jk)=ccci(jk,i) call print_output(i,jk,st,p0jk(jk)/p0,xv,psinv,btot,ak0,anpl,anpr, & - anprim,dens,tekev,alpha,tau,dids(jk,i),nhf,iokhawa, & - index_rt,ddr,ddi) - -! print error code - select case (ierr) - case(97) !+1 - print*,i,jk,' IERR = ', ierr,' N// = ',anpl - case(98) !+2 - print*,i,jk,' IERR = ', ierr,' N// = ',anpl - case(99) !+1*4 - print*,i,jk,' IERR = ', ierr,' Nwarm = ',anprre,anprim - case(94) !+2*4 - print*,i,jk,' IERR = ', ierr,' alpha < 0' - case(90) !+1*16 - print*,i,jk,' IERR = ', ierr,' fpp integration error' - case(91:93) !+2..4*16 - print*,i,jk,' IERR = ', ierr,' fcur integration error' - end select + anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi) end do @@ -251,15 +217,19 @@ contains if(nray > 1) call print_projxyzt(st,yw,0) end if -! test if trajectory integration must be stopped - call vmaxmin(tauv(:,i),nray,taumn,taumx) - if ((taumn > taucr) .or. (somein .and. allout)) then - pabs = sum(ppabs(:,i)) - icd = sum(ccci(:,i)) - istop = 1 - end if +! check for any error code and stop if necessary + call check_err(ierr,istop) +! test whether further trajectory integration is unnecessary + call vmaxmin(tau0,nray,taumn,taumx) + if ((taumn > taucr) .or. (somein .and. allout)) istop = 1 + if(istop == 1) exit end do + +! compute total absorbed power and driven current + if (i>nstep) i=nstep + pabs = sum(ppabs(:,i)) + icd = sum(ccci(:,i)) ! ======= main loop END ====== ! ======= post-proc BEGIN ====== @@ -289,8 +259,8 @@ contains ! ======= post-proc END ====== ! ======= free memory BEGIN ====== - call dealloc_beam(yw,ypw,xc,du1,gri,ggri, & - psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv) + call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, & + tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv) ! call unset_eqspl ! call unset_q ! call unset_rhospl @@ -301,11 +271,12 @@ contains end subroutine gray - subroutine vectinit(psjki,tauv,alphav,ppabs,dids,ccci,iiv) + subroutine vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv) use const_and_precisions, only : wp_, zero implicit none ! arguments - real(wp_), dimension(:,:), intent(out) :: psjki,tauv,alphav,ppabs,dids,ccci + real(wp_), dimension(:,:), intent(out) :: psjki,ppabs,ccci + real(wp_), dimension(:), intent(out) :: tau0,alphaabs0,dids0,ccci0 integer, dimension(:), intent(out) :: iiv !! common/external functions/variables ! integer :: jclosest @@ -317,13 +288,14 @@ contains ! anwcl(1:3)=0.0_wp_ ! xwcl(1:3)=0.0_wp_ - psjki = zero - tauv = zero - alphav = zero - ppabs = zero - dids = zero - ccci = zero - iiv = 1 + psjki = zero + ppabs = zero + ccci = zero + tau0 = zero + alphaabs0 = zero + dids0 = zero + ccci0 = zero + iiv = 1 end subroutine vectinit @@ -687,10 +659,10 @@ contains subroutine ywppla_upd(xv,anv,dgr,ddgr,sox,bres,xgcn,dery,psinv,dens,btot, & - xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm) + xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierr) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update -! use gray_params, only : igrad + use errcodes, only : pnpl implicit none ! arguments real(wp_), dimension(3), intent(in) :: xv,anv @@ -700,18 +672,27 @@ contains real(wp_), dimension(6), intent(out) :: dery real(wp_), intent(out) :: psinv,dens,btot,xg,yg,anpl,anpr real(wp_), intent(out) :: ddr,ddi,dersdst,derdnm + integer, intent(out) :: ierr ! local variables real(wp_) :: gr2,ajphi real(wp_), dimension(3) :: dgr2,bv,derxg,deryg real(wp_), dimension(3,3) :: derbv + real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ -! if(igrad == 1) then - gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2 - dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3)) -! end if + gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2 + dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3)) call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg,ajphi) call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & dery,anpl,anpr,ddr,ddi,dersdst,derdnm) + + ierr=0 + if( abs(anpl) > anplth1) then + if(abs(anpl) > anplth2) then + ierr=ibset(ierr,pnpl+1) + else + ierr=ibset(ierr,pnpl) + end if + end if end subroutine ywppla_upd @@ -1192,18 +1173,20 @@ contains end subroutine disp_deriv - subroutine alpha_effj(psinv,xg,yg,dens,tekev,zeff,ak0,bres,derdnm,anpl,anpr, & + subroutine alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, & sox,anprre,anprim,alpha,didp,nhmin,nhmax,iokhawa,ierr) use const_and_precisions, only : wp_,zero,pi,mc2=>mc2_ use gray_params, only : iwarm,ilarm,ieccd,imx + use coreprofiles, only : fzeff use equilibrium, only : sgnbphi use dispersion, only : harmnumber, warmdisp use eccd, only : setcdcoeff,eccdeff,fjch0,fjch,fjncl + use errcodes, only : palph use magsurf_data, only : fluxval implicit none ! arguments real(wp_),intent(in) ::psinv,ak0,bres - real(wp_),intent(in) :: xg,yg,tekev,dens,zeff,anpl,anpr,derdnm,sox + real(wp_),intent(in) :: xg,yg,tekev,dens,anpl,anpr,derdnm,sox real(wp_),intent(out) :: anprre,anprim,alpha,didp integer, intent(out) :: nhmin,nhmax,iokhawa integer, intent(out) :: ierr @@ -1211,9 +1194,9 @@ contains real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ ! local variables real(wp_) :: rbavi,rrii,rhop - integer :: lrm,ithn + integer :: lrm,ithn,ierrcd real(wp_) :: amu,ratiovgr,rbn,rbx - real(wp_) :: cst2,bmxi,bmni,fci + real(wp_) :: zeff,cst2,bmxi,bmni,fci real(wp_), dimension(:), allocatable :: eccdpar real(wp_) :: effjcd,effjcdav,akim,btot complex(wp_) :: ex,ey,ez @@ -1239,13 +1222,14 @@ contains ratiovgr=2.0_wp_*anpr/derdnm!*vgm alpha=2.0_wp_*akim*ratiovgr if(alpha: effjcdav [A m/W ] if(ieccd>0) then ! current drive computation + zeff=fzeff(psinv) ithn=1 if(lrm>nhmin) ithn=2 rhop=sqrt(psinv) @@ -1259,18 +1243,19 @@ contains ! cohen model call setcdcoeff(zeff,rbn,rbx,cst2,eccdpar) call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierr) + ithn,cst2,fjch,eccdpar,effjcd,iokhawa,ierrcd) case(2) ! no trapping call setcdcoeff(zeff,cst2,eccdpar) call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierr) + ithn,cst2,fjch0,eccdpar,effjcd,iokhawa,ierrcd) case default ! neoclassical model call setcdcoeff(zeff,rbx,fci,amu,rhop,cst2,eccdpar) call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax, & - ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierr) + ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd) end select + ierr=ierr+ierrcd !deallocate(eccdpar) effjcdav=rbavi*effjcd didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii) diff --git a/src/main.f90 b/src/main.f90 index 3eb97e9..0b857ba 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,7 +1,7 @@ program gray_main use const_and_precisions, only : wp_,one use graycore, only : gray - use gray_params, only : read_inputs,read_params, antctrl_type,eqparam_type, & + use gray_params, only : read_params, antctrl_type,eqparam_type, & prfparam_type,outparam_type,rtrparam_type,hcdparam_type use beams, only : read_beam0, read_beam1, read_beam2 use equilibrium, only : read_equil_an,read_eqdsk,change_cocos,eq_scal, & @@ -32,8 +32,7 @@ program gray_main real(wp_) :: rwallm, rmxm, r0m, z0m, dzmx ! ======= read parameters BEGIN ======= - call read_inputs('graynew.data',antp,eqp,rwallm,prfp,outp) - call read_params('gray_params.data',rtrp,hcdp) + call read_params('gray_params.data',rtrp,hcdp,antp,eqp,rwallm,prfp,outp) ! ======= read parameters END ======= ! ======= read input data BEGIN ======= diff --git a/src/pec.f90 b/src/pec.f90 index 61d3eb6..cee1a67 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -102,7 +102,7 @@ contains do jk=1,nray ii=iiv(jk) if (ii < nstep ) then - if(psjki(jk,ii+1) /= zero) ii=ii+1 + if(psjki(jk,ii+1) /= zero) ii=ii+1 !!! CHECK end if xxi=zero ypt=zero