improve error handling in the gray_main routine

- rename errocodes → gray_errors

- restructure the errors into a `error_spec` type

- make the list of errors easily extensible

- rewrite the `print_errn`, `print_errhcd` (now `print_err_raytracing`,
  `print_err_ecrh_cd`) subroutine to handle arbitrary errors

- add functions to easily manipulate errors
  (`raise_error`, `has_error`, `is_critical`)

- remove print statements from quadpack

- log all errors to stderr using the logger module

src/dispersion.f90

@@ -226,7 +226,8 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox,        &
   !
   ! Reference: https://doi.org/10.13182/FST08-A1660
   use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
-
+  use gray_errors, only : gray_error, warmdisp_convergence, warmdisp_result, &
+                          raise_error
   implicit none
 
   ! subroutine arguments
@@ -246,15 +247,8 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox,        &
 
   ! Outputs
 
-  ! error code:
+  ! GRAY error code
-  !   0: converged within `max_iters`
+  integer(kind=gray_error), intent(out) :: error
-  !   1: failed to converge, returned fallback value
-  !   2: failed to converge, returned last value
-  !   4: converged to invalid value
-  !   5: failed to converge, fallback is also invalid
-  !   6: failed to converge, last value is invalid
-  ! Note: invalid means either IEEE 754 NaN or ±Infinity.
-  integer, intent(out) :: error
   ! orthogonal refractive index: N⊥
   complex(wp_), intent(out) :: Npr
   ! polarization vector: ε (local cartesian frame, z∥B̅)
@@ -262,7 +256,7 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox,        &
 
   ! Parameters
 
-  ! switch for the dielectic tensor model:
+  ! switch for the dielectric tensor model:
   !  1: weakly relativistic
   !  2: fully relativistic (faster variant)
   !  3: fully relativistic (slower variant)
@@ -315,7 +309,7 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox,        &
   Npr2_err = 1
 
   ! Compute the N⊥-independent part of the dielectric tensor
-  call dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl)
+  call dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl, error)
 
   iterate_to_solve_disperion: do i = 1, max_iters
     Npr_prev = sqrt(Npr2_prev)
@@ -433,23 +427,23 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox,        &
   if (i > max_iters) then
     if (fallback) then
       ! first try failed, use fallback
-      error = 1
+      error = raise_error(error, warmdisp_convergence, 0)
       Npr2 = Npr2_fallback
     else
       ! first try failed, no retry
-      error = 2
+      error = raise_error(error, warmdisp_convergence, 1)
     end if
   end if
 
   ! Catch NaN and ±Infinity values
   if (.not. ieee_is_finite(Npr2%re) .or. .not. ieee_is_finite(Npr2%im)) then
     Npr2 = 0
-    error = error + 4
+    error = raise_error(error, warmdisp_result, 0)
   end if
 
   if (Npr2%re < 0 .and. Npr2%im < 0) then
     Npr2 = sqrt(0.5_wp_)  ! Lorenzo, what should we have here??
-    error = 99
+    error = raise_error(error, warmdisp_result, 1)
   end if
 
   ! Final result
@@ -482,11 +476,12 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox,        &
 end subroutine warmdisp
 
 
-subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl)
+subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl, error)
   ! Returns the N⊥-independent parts of the dielectric tensor, namely
   ! ε₃₃⁰ and the 3x3 tensors ε̅^(l) for l=1,nlarmor.
   !
   ! Reference: https://doi.org/10.13182/FST08-A1660
+  use gray_errors, only : gray_error
 
   implicit none
 
@@ -503,7 +498,7 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl)
 
   ! Parameters
 
-  ! switch for the dielectic tensor model:
+  ! switch for the dielectric tensor model:
   !  1: weakly relativistic
   !  2: fully relativistic (faster variant)
   !  3: fully relativistic (slower variant)
@@ -518,6 +513,8 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl)
   complex(wp_), intent(out) :: e330
   ! tensors ε̅_ij^(l), defined in eq. 10
   complex(wp_), intent(out) :: epsl(3,3,nlarmor)
+  ! GRAY error code
+  integer(kind=gray_error), intent(out) :: error
 
   ! local variables
   integer :: l, n
@@ -536,10 +533,10 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl)
   if (model == 1) then
     ! Compute the (combinations of) Shkarofsky functions F_q(n)
     ! in which the functions Q⁺_hl(n), Q⁻_hl(n) can be expressed
-    call fsup(nlarmor, Y, Npl, mu, Fp, Fm)
+    call fsup(nlarmor, Y, Npl, mu, Fp, Fm, error)
   else
     ! Compute the real and imaginary parts of the integrals
-    ! defining the functions Q_hl(n), Q_hl(n)
+    ! defining the functions Q_hl(n), Q_hl(-n)
     block
       real(wp_) :: cr, ci, cmxw
       cmxw = 1 + 15/(8*mu) + 105/(128 * mu**2)
@@ -552,7 +549,8 @@ subroutine dielectric_tensor(X, Y, mu, Npl, model, nlarmor, e330, epsl)
         call hermitian(rr, Y, mu, Npl, cr, model, nlarmor)
       else
         ! slow variant
-        call hermitian_2(rr, Y, mu, Npl, cr, model, nlarmor)
+        call hermitian_2(rr, Y, mu, Npl, cr, model, nlarmor, error)
+        if (error /= 0) error = 2
       end if
       ! imaginary part
       call antihermitian(ri, Y, mu, Npl, ci, nlarmor)
@@ -868,8 +866,9 @@ end subroutine hermitian
 !
 !
 !
-subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm)
+subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm,error)
-  use quadpack, only : dqagsmv !dqagimv
+  use gray_errors, only : gray_error, dielectric_tensor, raise_error
+  use quadpack,    only : dqagsmv
   implicit none
 ! local constants
   integer,parameter :: lw=5000,liw=lw/4,npar=7
@@ -878,6 +877,7 @@ subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm)
   integer :: lrm,fast
   real(wp_) :: yg,mu,npl,cr
   real(wp_), dimension(-lrm:lrm,0:2,0:lrm) :: rr
+  integer(kind=gray_error), intent(out) :: error
 ! local variables
   integer :: n,m,ih,k,n1,nn,llm,neval,ier,last,ihmin
   integer, dimension(liw) :: iw
@@ -919,9 +919,9 @@ subroutine hermitian_2(rr,yg,mu,npl,cr,fast,lrm)
       if(n.eq.0.and.m.eq.0) ihmin=2
       do ih=ihmin,2
         apar(7) = real(ih,wp_)
-!        call dqagimv(fhermit,bound,2,apar,npar,epsa,epsr,resfh,
         call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa,epsr,resfh, &
                     epp,neval,ier,liw,lw,last,iw,w)
+        if (ier /= 0) error = raise_error(error, dielectric_tensor, 1)
         rr(n,ih,m) = resfh
       end do
     end do
@@ -1216,13 +1216,16 @@ subroutine expinit
 end subroutine expinit
 
 
-pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm)
+pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm, error)
+  use gray_errors, only : gray_error, dielectric_tensor, raise_error
+
   implicit none
 
   ! subroutine arguments
   integer, intent(in) :: lrm
   real(wp_), intent(in) :: yg, npl, mu
   complex(wp_), intent(out), dimension(0:lrm,0:2) :: cefp, cefm
+  integer(kind=gray_error), intent(out) :: error
 
   ! local constants
   real(wp_), parameter :: apsicr=0.7_wp_
@@ -1256,7 +1259,9 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm)
       y0=phim
     end if
     call zetac (xp,yp,zrp,zip,iflag)
+    if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0)
     call zetac (xm,ym,zrm,zim,iflag)
+    if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0)
 !
     czp=dcmplx(zrp,zip)
     czm=dcmplx(zrm,zim)
@@ -1273,6 +1278,7 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm)
       cphi=phim
       if(alpha.lt.0) cphi=-im*phim
       call zetac (x0,y0,zr0,zi0,iflag)
+      if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0)
       cz0=dcmplx(zr0,zi0)
       cdz0=2.0_wp_*(one-cphi*cz0)
       cf32=cdz0
@@ -1319,7 +1325,9 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm)
         y0=phim
       end if
       call zetac (xp,yp,zrp,zip,iflag)
+      if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0)
       call zetac (xm,ym,zrm,zim,iflag)
+      if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0)
 !
       czp=dcmplx(zrp,zip)
       czm=dcmplx(zrm,zim)
@@ -1336,6 +1344,7 @@ pure subroutine fsup(lrm, yg, npl, mu, cefp, cefm)
         cphi=phim
         if(alpha.lt.0) cphi=-im*phim
         call zetac (x0,y0,zr0,zi0,iflag)
+        if (iflag /= 0) error = raise_error(error, dielectric_tensor, 0)
         cz0=dcmplx(zr0,zi0)
         cdz0=2.0_wp_*(one-cphi*cz0)
         cf32=cdz0

src/eccd.f90

@@ -176,7 +176,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 gray_errors, only : fpp_integration, fcur_integration, raise_error
     use quadpack, only : dqagsmv
     implicit none
 ! local constants
@@ -258,7 +258,7 @@ contains
                    epp,neval,ier,liw,lw,last,iw,w)
              
       if (ier.gt.0) then
-        ierr=ibset(ierr,pcdfp)
+        ierr=raise_error(ierr, fpp_integration)
         return
       end if
 
@@ -302,7 +302,7 @@ contains
               if (abs(resji).lt.1.0e-10_wp_) then
                 resji=0.0_wp_
               else
-                ierr=ibset(ierr,pcdfc+iokhawa+i)
+                ierr=raise_error(ierr, fcur_integration, iokhawa + i)
                 return
               end if
             end if

src/errcodes.f90

@@ -1,83 +0,0 @@
-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(inout) :: istop
-!     if(ibits(ierr,pnpl, lnpl )>1  .or.    & ! N// too large
-!        ibits(ierr,palph,lalph)==1) then      ! alpha < 0
-!       istop = 1
-!    else
-!      istop = 0
-!     end if
-
-!     if(ibits(ierr,pnpl, lnpl )>1) istop = 1  ! N// too large
-    if(ibits(ierr,palph,lalph)==1) istop = 1 ! alpha < 0
-    
-  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**pnprre,'  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_core.f90

@@ -15,7 +15,7 @@ contains
                              print_parameters
     use beams,        only : xgygcoeff, launchangles2n
     use beamdata,     only : pweight, rayi2jk
-    use errcodes,     only : check_err, print_errn, print_errhcd
+    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 pec,          only : pec_init, spec, postproc_profiles, dealloc_pec, &
@@ -62,7 +62,7 @@ contains
     ! i: integration step, jk: global ray index
     integer :: i, jk
 
-    integer :: iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt
+    integer :: iox,nharm,nhf,nnd,iokhawa,ierrn,ierrhcd,index_rt
     integer :: ip,ib,iopmin,ipar,iO
     integer :: igrad_b,istop_pass,nbeam_pass,nlim
     logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff
@@ -271,7 +271,6 @@ contains
           if(igrad_b == 1) call gradi_upd(yw,ak0,xc,du1,gri,ggri)
 
           error = 0
-          istop = 0                                                                        ! stop flag for current beam
           iopmin = 10
 
           ! =========== rays loop BEGIN ===========
@@ -287,7 +286,7 @@ contains
             ! update global error code and print message
             if(ierrn/=0) then
               error = ior(error,ierrn)
-              call print_errn(ierrn,i,anpl)
+              call print_err_raytracing(ierrn, i, anpl)
             end if
 
             ! check entrance/exit plasma/wall
@@ -379,7 +378,7 @@ contains
                   ierrn,igrad_b)                                                           ! * update derivatives after reflection
               if(ierrn/=0) then                                                            ! * update global error code and print message
                 error = ior(error,ierrn)
-                call print_errn(ierrn,i,anpl)
+                call print_err_raytracing(ierrn, i, anpl)
               end if
 
               if(jk == 1 .and. ip == 1) then                                               ! * 1st pass, polarization angles at reflection for central ray
@@ -437,11 +436,11 @@ contains
                                 alpha, didp, nharm, nhf, iokhawa, ierrhcd)
                 anprre = Npr_warm%re
                 anprim = Npr_warm%im
-              end block
                 if(ierrhcd /= 0) then
                   error = ior(error, ierrhcd)
-                call print_errhcd(ierrhcd,i,anprre,anprim,alpha)
+                  call print_err_ecrh_cd(ierrhcd, i, Npr_warm, alpha)
                 end if
+              end block
             else
               tekev=zero
               alpha=zero
@@ -498,12 +497,10 @@ contains
               call print_projxyzt(stv,yw,0)
           end if
 
-          ! check for any error code and stop if necessary
-          call check_err(error,istop)
           ! test whether further trajectory integration is unnecessary
           call vmaxmin(tau1+tau0+lgcpl1,params%raytracing%nray,taumn,taumx)                ! test on tau + coupling
 
-          if(istop == 1) then                                                              ! stop propagation for current beam
+          if(is_critical(error)) then                                                      ! stop propagation for current beam
             istop_pass = istop_pass +1                                                     ! * +1 non propagating beam
             if(ip < params%raytracing%ipass) call turnoffray(0,ip,ib,iroff)                ! * de-activate derived beams
             exit
@@ -1029,7 +1026,7 @@ contains
      bv,xg,yg,derxg,anpl,anpr,ddr,ddi,dersdst,derdnm,error,igrad)
 ! Compute right-hand side terms of the ray equations (dery)
 ! used after full R-K step and grad(S_I) update
-    use errcodes, only : pnpl
+    use gray_errors, only : raise_error, large_npl
     implicit none
 ! arguments
     real(wp_), dimension(3), intent(in) :: xv,anv
@@ -1053,13 +1050,12 @@ contains
     call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, &
                     dery,anpl,anpr,ddr,ddi,dersdst,derdnm)
 
-    error=0
-    if(  abs(anpl) > anplth1) then
     if (abs(anpl) > anplth2) then
-        error=ibset(error,pnpl+1)
+      error = raise_error(error, large_npl, 1)
+    else if (abs(anpl) > anplth1) then
+      error = raise_error(error, large_npl, 0)
     else
-        error=ibset(error,pnpl)
+      error = 0
-      end if
     end if
   end subroutine ywppla_upd
 
@@ -1711,7 +1707,7 @@ contains
     use equilibrium,          only : sgnbphi
     use dispersion,           only : harmnumber, warmdisp
     use eccd,                 only : setcdcoeff, eccdeff, fjch0, fjch, fjncl
-    use errcodes,             only : palph
+    use gray_errors,          only : negative_absorption, raise_error
     use magsurf_data,         only : fluxval
 
     implicit none
@@ -1814,7 +1810,7 @@ contains
     end block
 
     if (alpha < 0) then
-      error = ibset(error, palph)
+      error = raise_error(error, negative_absorption)
       return
     end if
 

src/gray_errors.f90

@@ -0,0 +1,208 @@
+! This module contains the error codes handled by the gray_main subroutine
+!
+! An error is actually a 32bit integer bitmask that can encode several
+! error cases. The individual errors can be raised using the `raise_error`
+! function and combined simply with the intrinsic `ior` function.
+module gray_errors
+
+  use logger, only : log_error
+  use, intrinsic :: iso_fortran_env, only : int32
+
+  implicit none
+
+  ! The type of a GRAY error is `integer(kind=gray_error)`
+  integer, parameter :: gray_error = int32
+
+  ! Specification of a GRAY error
+  type error_spec
+    integer       :: offset    ! bitmask offset
+    integer       :: subcases  ! number of error subcases (max: 10)
+    character(32) :: mod       ! module raising the error
+    character(32) :: proc      ! procedure raising the error
+    character(64) :: msg(0:9)  ! error message (one for each subcase)
+  end type
+
+  ! macros used for defining errors
+# define _ ,
+# define list(x) reshape([character(64) :: x], [10], [''])
+# define after(x) x%offset + x%subcases
+
+  ! All GRAY errors
+
+  type(error_spec), parameter :: large_npl =        &
+    error_spec(offset=0, subcases=2,                &
+               mod='gray_core', proc='gray_main',   &
+               msg=list('N∥ is too large (>0.99)'_
+                        'N∥ is too large (>1.05)'))
+
+  type(error_spec), parameter :: dielectric_tensor =  &
+    error_spec(offset=after(large_npl), subcases=2,   &
+               mod='gray_core', proc='gray_main',     &
+               msg=list('ε tensor, overflow in `fsup`'_
+                        'ε tensor, integration error in `hermitian_2`'))
+
+  type(error_spec), parameter :: warmdisp_convergence =      &
+    error_spec(offset=after(dielectric_tensor), subcases=2,  &
+               mod='dispersion', proc='warmdisp',            &
+               msg=list('failed to converge, returned fallback value'_
+                        'failed to converge, returned last value'))
+
+  type(error_spec), parameter :: warmdisp_result =              &
+    error_spec(offset=after(warmdisp_convergence), subcases=2,  &
+               mod='dispersion', proc='warmdisp',               &
+               msg=list('final N⊥² is NaN or ±Infinity'_
+                        'final N⊥² in 3rd quadrant'))
+
+  type(error_spec), parameter :: negative_absorption =     &
+    error_spec(offset=after(warmdisp_result), subcases=1,  &
+               mod='gray_core', proc='alpha_effj',         &
+               msg=list('negative absorption coeff.'))
+
+  type(error_spec), parameter :: fpp_integration =             &
+    error_spec(offset=after(negative_absorption), subcases=1,  &
+               mod='eccd', proc='eccdeff',                     &
+               msg=list('fpp integration error'))
+
+  type(error_spec), parameter :: fcur_integration =        &
+    error_spec(offset=after(fpp_integration), subcases=3,  &
+               mod='eccd', proc='eccdeff',                 &
+               msg=list('fcur integration error (no trapping)'_
+                        'fcur integration error (1st trapping region)'_
+                        'fcur integration error (2st trapping region)'))
+
+  ! Errors occuring during raytracing
+  type(error_spec), parameter :: raytracing_errors(*) = [large_npl]
+
+  ! Errors occuring during absorption and current drive computations
+  type(error_spec), parameter :: ecrh_cd_errors(*) = &
+    [dielectric_tensor,    &
+     warmdisp_convergence, &
+     warmdisp_result,      &
+     negative_absorption,  &
+     fpp_integration,      &
+     fcur_integration]
+
+contains
+
+  pure function is_critical(error)
+    ! Checks whether critical errors have occurred
+
+    implicit none
+
+    ! subroutines arguments
+    integer(kind=gray_error), intent(in) :: error
+    logical                              :: is_critical
+
+    is_critical = has_error(error, negative_absorption)
+  end function is_critical
+
+
+  pure function has_error(error, spec)
+    ! Checks whether the `error` bitmask contains the error given by `spec`
+    implicit none
+
+    ! function arguments
+    integer(kind=gray_error), intent(in) :: error
+    type(error_spec),         intent(in) :: spec
+    logical                              :: has_error
+
+    has_error = ibits(error, spec%offset, spec%subcases) == 1
+  end function has_error
+
+
+  pure function raise_error(error, spec, subcase)
+    ! Raise the bits of error `spec` (with optional `subcase` number)
+    ! in the `error` bitmask.
+    implicit none
+
+    ! function arguments
+    integer(kind=gray_error), intent(in) :: error
+    type(error_spec),         intent(in) :: spec
+    integer,                  intent(in), optional :: subcase
+    integer(kind=gray_error)             :: raise_error
+
+    raise_error = ibset(error, spec%offset &
+                               + merge(subcase, 0, present(subcase)))
+  end function raise_error
+
+
+  subroutine print_err_raytracing(error, step, Npl)
+    ! Pretty prints raytracing errors
+    !
+    ! The error and some context (integration step, N∥)
+    ! is logged to the stderr using the logger module.
+
+    use const_and_precisions, only : wp_
+    implicit none
+
+    ! subroutines arguments
+    integer,   intent(in) :: error, step
+    real(wp_), intent(in) :: Npl
+
+    ! local variables
+    integer :: slice        ! a slice of the bitmask
+    character(256) :: line  ! formatted log line
+    type(error_spec) :: spec
+    integer :: i, j
+
+    ! format specifier of the log line
+    character(*), parameter :: fmt = &
+      '(a,": ","error=",g0," N∥=",g0.3," step=",g0)'
+
+    ! iterate on the known errors
+    do i = 1, size(raytracing_errors)
+      spec = raytracing_errors(i)
+      slice = ibits(error, spec%offset, spec%subcases)
+      if (slice == 0) cycle
+
+      ! iterate on the subcases
+      do j = 0, spec%subcases - 1
+        if (ibits(slice, j, 1) == 0) cycle
+        write(line, fmt) trim(spec%msg(j)), slice * 2**j, Npl, step
+        call log_error(line, mod=spec%mod, proc=spec%proc)
+      end do
+    end do
+  end subroutine print_err_raytracing
+
+
+  subroutine print_err_ecrh_cd(error, step, Npr, alpha)
+    ! Pretty prints ECRH & CD errors
+    !
+    ! The error and some context (integration step, N⊥, α)
+    ! is logged to the stderr using the logger module.
+
+    use const_and_precisions, only : wp_
+    implicit none
+
+    ! subroutines arguments
+    integer,       intent(in) :: error, step
+    complex(wp_),  intent(in) :: Npr
+    real(wp_),     intent(in) :: alpha
+
+    ! local variables
+    integer :: slice        ! a slice of the bitmask
+    character(256) :: line  ! formatted log line
+    type(error_spec) :: spec
+    integer :: i, j
+
+    ! format specifier of the log line
+    character(*), parameter :: fmt = &
+      '(a,": ", "error=",g0, " N⊥=",f0.0,sp,f0.0,"i", " α=",g0, " step=",g0)'
+
+    ! iterate on the known errors
+    do i = 1, size(ecrh_cd_errors)
+      spec = ecrh_cd_errors(i)
+      slice = ibits(error, spec%offset, spec%subcases)
+      if (slice == 0) cycle
+
+      ! iterate on the subcases
+      do j = 0, spec%subcases - 1
+        if (ibits(slice, j, 1) == 0) cycle
+        write(line, fmt) trim(spec%msg(j)), slice * 2**j, &
+                         Npr%re, Npr%im, alpha, step
+        call log_error(line, mod=spec%mod, proc=spec%proc)
+      end do
+    end do
+  end subroutine print_err_ecrh_cd
+
+end module gray_errors

src/quadpack.f90

@@ -195,7 +195,6 @@ contains
       lvl = 0
     end if
     if(ier==6) lvl = 1
-    if(ier/=0) print*,'habnormal return from dqags',ier,lvl
   end subroutine dqags
 
   subroutine dqagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval, &
@@ -1358,7 +1357,6 @@ contains
         neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last)
 !
     end if
-    if(ier/=0) print*,'habnormal return from dqagi'
   end subroutine dqagi
 
   subroutine dqagie(f,bound,inf,epsabs,epsrel,limit,result,abserr, &
@@ -2248,7 +2246,6 @@ contains
       lvl = 0
     end if
     if(ier.eq.6) lvl = 1
-    if(ier.ne.0) print*,'habnormal return from dqaps',ier,lvl
   end subroutine dqagp
 
   subroutine dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,      &
@@ -3016,7 +3013,6 @@ contains
       lvl = 0
     end if
     if(ier==6) lvl = 1
-    if(ier/=0) print*,'habnormal return from dqags',ier,lvl
   end subroutine dqagsmv
 
   subroutine dqagsemv(f,a,b,apar,np,epsabs,epsrel,limit,result,abserr,neval, &
@@ -3865,7 +3861,6 @@ contains
         neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last)
 !
     end if
-    if(ier/=0) print*,'habnormal return from dqagi'
   end subroutine dqagimv
 
   subroutine dqagiemv(f,bound,inf,apar,np,epsabs,epsrel,limit,result,abserr, &