diff --git a/SOURCE/radiation_model_mod.f90 b/SOURCE/radiation_model_mod.f90
index a7bdd404b8feb3c7f9324acaac9209f0e46a1d08..710f20fc0a185ef757acaf78535240b897991ad5 100644
--- a/SOURCE/radiation_model_mod.f90
+++ b/SOURCE/radiation_model_mod.f90
@@ -6925,7 +6925,7 @@
        TYPE(c_ptr)                               ::  gridsurf_rma_p !< allocated c pointer
        INTEGER(iwp)                              ::  minfo          !< MPI RMA window info handle
 #endif
-
+       REAL(wp), PARAMETER                       ::  eps_lad = 1E-10_wp !< epsilon for value comparison
 !
 !--     precalculate face areas for different face directions using normal vector
         DO d = 0, nsurf_type
@@ -6957,13 +6957,12 @@
 !
 !--                Find topography top index
                    k_topo = topo_top_ind(j,i,0)
-
+                   npcbl = npcbl + COUNT( lad_s(:,j,i) > eps_lad )
                    DO k = nzt+1, 0, -1
-                       IF ( lad_s(k,j,i) /= 0.0_wp )  THEN
+                       IF ( lad_s(k,j,i) > eps_lad )  THEN
 !--                        we are at the top of the pcs
                            pct(j,i) = k + k_topo
                            pch(j,i) = k
-                           npcbl = npcbl + pch(j,i)
                            EXIT
                        ENDIF
                    ENDDO
@@ -7082,9 +7081,11 @@
                    k_topo = topo_top_ind(j,i,0)
 
                    DO k = k_topo + 1, pct(j,i)
-                       ipcgb = ipcgb + 1
-                       gridpcbl(k,j,i) = ipcgb
-                       pcbl(:,ipcgb) = (/ k, j, i /)
+                       IF ( lad_s(k - k_topo, j, i) > eps_lad ) THEN
+                          ipcgb = ipcgb + 1
+                          gridpcbl(k,j,i) = ipcgb
+                          pcbl(:,ipcgb) = (/ k, j, i /)
+                       END IF
                    ENDDO
                ENDDO
            ENDDO