diff --git a/src/graycore.f90 b/src/graycore.f90
index 7a1e0b1..9d7b1e6 100644
--- a/src/graycore.f90
+++ b/src/graycore.f90
@@ -512,12 +512,16 @@ contains
         pabs_pass(iox) = pabs_pass(iox) + pabs_beam                                        ! 0D results for current pass, sum on O/X mode beams
         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 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
+        if(ip.lt.ipass .and. iopmin.gt.2) then                                             ! not last pass AND at least one ray re-entered plasma
+          cpl_beam1 = sum(p0ray * exp(-tau0) * cpls(:,iO)/cpl1, MASK=iop.gt.2) / &
+                          sum(p0ray * exp(-tau0), MASK=iop.gt.2)                           ! * average O-mode coupling for next beam (on active rays)
+          cpl_beam2 = one - cpl_beam1                                                      ! * average X-mode coupling for next beam
+          
+          if(iop(1).gt.2) then                                                             ! * central ray O/X-mode coupling for next beam
+            cpl_cbeam1 = cpls(1,iO)/cpl1(1)
+            cpl_cbeam2 = one - cpl_cbeam1
+          end if
+        else                                                                               ! last pass OR no ray re-entered plasma
           cpl_beam1 = zero
           cpl_beam2 = zero
         end if
@@ -531,8 +535,9 @@ contains
         write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs_beam
         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 (ctr ray, O/X):',cpl_cbeam1,cpl_cbeam2       ! central ray coupling for next O/X beams
+          write(*,'(a,2(f9.4,1x))') 'Coupling (average, O/X):',cpl_beam1,cpl_beam2         ! average coupling for next O/X beams (=0 if no ray re-entered plasma)
+          if(iop(1).gt.2) 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, &