* reduced arrays required for dI/ds,dP/ds integration

* new module for error handling
* input parameters collected in a single file
* fixed uninitialized pabs,icd

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
+dispersion.o: const_and_precisions.o eierf.o errcodes.o math.o quadpack.o
-eccd.o: const_and_precisions.o conical.o magsurf_data.o dierckx.o numint.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

src/beamdata.f90

@@ -8,16 +8,17 @@ module beamdata
 
 contains
 
-  subroutine init_rtr(rtrparam,ywork,ypwork,xc,du1,gri,ggri, &
+  subroutine init_rtr(rtrparam,ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
-                      psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
+                      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, &
+    call alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
-                    psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
+                    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, &
+  subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
-                    psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
+                        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, &
+    call dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
-                    psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
+                      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),                    &
+             xc(3,nrayth,nrayr),  du1(3,nrayth,nrayr),                   &
-             psjki(nray,nstep),  tauv(nray,nstep),  alphav(nray,nstep),  &
+             psjki(nray,nstep),   ppabs(nray,nstep),   ccci(nray,nstep), &
-             ppabs(nray,nstep),  dids(nray,nstep),    ccci(nray,nstep),  &
+             tau0(nray),   alphaabs0(nray), dids0(nray),    ccci0(nray), &
-             p0jk(nray), ext(nray), eyt(nray), iiv(nray))
+             p0jk(nray),         ext(nray),   eyt(nray),      iiv(nray))
   end subroutine alloc_beam
 
 
-  subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri, &
+  subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
-                    psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
+                          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(ywork))     deallocate(ywork)
-    if (allocated(ypwork))  deallocate(ypwork)
+    if (allocated(ypwork))    deallocate(ypwork)
-    if (allocated(xc))      deallocate(xc)
+    if (allocated(xc))        deallocate(xc)
-    if (allocated(du1))     deallocate(du1)
+    if (allocated(du1))       deallocate(du1)
-    if (allocated(gri))     deallocate(gri)
+    if (allocated(gri))       deallocate(gri)
-    if (allocated(ggri))    deallocate(ggri)
+    if (allocated(ggri))      deallocate(ggri)
-    if (allocated(psjki))   deallocate(psjki)
+    if (allocated(psjki))     deallocate(psjki)
-    if (allocated(tauv))    deallocate(tauv)
+    if (allocated(ppabs))     deallocate(ppabs)
-    if (allocated(alphav))  deallocate(alphav)
+    if (allocated(ccci))      deallocate(ccci)
-    if (allocated(ppabs))   deallocate(ppabs)
+    if (allocated(tau0))      deallocate(tau0)
-    if (allocated(dids))    deallocate(dids)
+    if (allocated(alphaabs0)) deallocate(alphaabs0)
-    if (allocated(ccci))    deallocate(ccci)
+    if (allocated(dids0))     deallocate(dids0)
-    if (allocated(p0jk))    deallocate(p0jk)
+    if (allocated(ccci0))     deallocate(ccci0)
-    if (allocated(ext))     deallocate(ext)
+    if (allocated(p0jk))      deallocate(p0jk)
-    if (allocated(eyt))     deallocate(eyt)
+    if (allocated(ext))       deallocate(ext)
-    if (allocated(iiv))     deallocate(iiv)
+    if (allocated(eyt))       deallocate(eyt)
+    if (allocated(iiv))       deallocate(iiv)
   end subroutine dealloc_beam
 
   subroutine pweight(p0,p0jk)

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 &
+        err=1
-         &disabled.',e12.5)") xg,yg,sqrt(abs(npr2)),npl
         imxx=1
       else
-        write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence &
+        err=2
-         &failed.',e12.5)") xg,yg,sqrt(abs(npr2)),npl
         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)

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

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

src/gray_params.f90

@@ -48,12 +48,14 @@ module gray_params
 
 contains
   
-  subroutine read_inputs(filenm,antctrl,eqparam,rwall,prfparam,outparam,unit)
+  subroutine read_params(filenm,rtrparam,hcdparam,antctrl,eqparam,rwall, &
-    use const_and_precisions, only : wp_
+                         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
 

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_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0
-    real(wp_) :: tau0,alphaabs0,dids0,ccci0
+    real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx
-    real(wp_) :: tau,pow,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, &
+    call init_rtr(rtrp,yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
-                  psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
+                  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. ! becomes true if at least part of the beam enters the plasma
-    somein = .false.
-    istop = 0
 ! 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. ! becomes false if at least part of the beam is inside the plsama
-      allout = .true.
+      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)
+           psinv,dens,btot,xg,yg,anpl,anpr,ddr,ddi,dersdst,derdnm,ierrn)
-
+! update global error code and print message
-        if(  abs(anpl) >  anplth1) then
+        if (ierrn/=0) then
-          if(abs(anpl) <= anplth2) then
+          ierr = ior(ierr,ierrn)
-            ierr=97
+          call print_errn(ierrn,i,anpl)
-!            igrad=0
-          else 
-            ierr=98
-            istop=1
-          end if 
-        else
-          ierr=0
         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<one .and. zzm>=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
+        if(ierrn==0 .and. iwarm>0 .and. ins_pl .and. tau0(jk)<=taucr) then
-!          print*,i,jk,rayi2jk(jk),psinv,zzm,anpl
           tekev=temp(psinv)
-          if (ieccd> 0) zeff=fzeff(psinv)
+          call alpha_effj(psinv,xg,yg,dens,tekev,ak0,bres,derdnm,anpl,anpr, &
-          call alpha_effj(psinv,xg,yg,dens,tekev,zeff,ak0,bres,derdnm, &
+             sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierrhcd) 
-             anpl,anpr,sox,anprre,anprim,alpha,didp,nharm,nhf,iokhawa,ierr) 
+          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
+        tau=tau0(jk)+0.5_wp_*(alphaabs0(jk)+alpha)*dersdst*dst
-        tauv(jk,i)=tau
-        alphav(jk,i)=alpha
         pow=p0jk(jk)*exp(-tau) !*exp(-tau1v(jk))
         ppabs(jk,i)=p0jk(jk)-pow
 
-        dids(jk,i)=didp*pow*alpha
+        dids=didp*pow*alpha
-        ccci(jk,i)=ccci0+0.5_wp_*(dids0+dids(jk,i))*dersdst*dst
+        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,  &
+           anprim,dens,tekev,alpha,tau,dids,nhf,iokhawa,index_rt,ddr,ddi)
-           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
 
       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
+! check for any error code and stop if necessary
-      call vmaxmin(tauv(:,i),nray,taumn,taumx)
+      call check_err(ierr,istop)
-      if ((taumn > taucr) .or. (somein .and. allout)) then
+! test whether further trajectory integration is unnecessary
-        pabs  = sum(ppabs(:,i))
+      call vmaxmin(tau0,nray,taumn,taumx)
-        icd   = sum(ccci(:,i))
+      if ((taumn > taucr) .or. (somein .and. allout)) istop = 1
-        istop = 1
+
-      end if
       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, &
+    call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
-                    psjki,tauv,alphav,ppabs,dids,ccci,p0jk,ext,eyt,iiv)
+                      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
+      psjki     = zero
-      tauv = zero
+      ppabs     = zero
-      alphav = zero
+      ccci      = zero
-      ppabs = zero
+      tau0      = zero
-      dids = zero
+      alphaabs0 = zero
-      ccci = zero
+      dids0     = zero
-      iiv = 1
+      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
-      gr2 = dgr(1)**2 + dgr(2)**2 + dgr(3)**2
+    dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
-      dgr2 = 2*(dgr(1)*ddgr(:,1) + dgr(2)*ddgr(:,2) + dgr(3)*ddgr(:,3))
-!    end if
     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<zero) then
-          ierr=94
+          ierr=ibset(ierr,palph)
           return
         end if
           
 !       calcolo della efficienza <j/p>:  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)

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,antp,eqp,rwallm,prfp,outp)
-  call read_params('gray_params.data',rtrp,hcdp)
 ! ======= read parameters  END  =======
 
 ! ======= read input data BEGIN =======

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