diff --git a/src/neptuneClass.f90 b/src/neptuneClass.f90
index e1f8c61..e77fb8a 100644
--- a/src/neptuneClass.f90
+++ b/src/neptuneClass.f90
@@ -515,6 +515,8 @@ subroutine initialize_input_array(this)
         call this%set_input(parName=C_OPT_HARMONICS, valType='boolean', initFlag=.true.)
         call this%set_input(parName=C_OPT_PROGRESS, valType='boolean', initFlag=.true.)
         call this%set_input(parName=C_OPT_PROGRESS, val = toString(this%write_progress), set = .true.)
+        call this%set_input(parName=C_BOUNDARY_CHECK, valType='boolean', initFlag=.true.)
+        call this%set_input(parName=C_BOUNDARY_CHECK, val = "true", set = .true.)
 
         if(isControlled()) then
             call checkOut(csubid)
@@ -1424,7 +1426,7 @@ integer function setNeptuneVar_char(        &
                  C_OUTPUT_ACO,     C_OUTPUT_OSC,     C_OUTPUT_CSV,     C_OUTPUT_GLL,      &
                  C_COV_PROP,       C_COV_SUN,        C_COV_MOON,       C_COV_DRAG,        &
                  C_COV_SRP,        C_OUTPUT_VAR_ECI, C_OPT_HARMONICS, C_OPT_PROGRESS,     &
-                 C_OUTPUT_VAR_UVW, C_OUTPUT_COV_ECI, C_OUTPUT_COV_UVW, C_OUTPUT_AMN)
+                 C_OUTPUT_VAR_UVW, C_OUTPUT_COV_ECI, C_OUTPUT_COV_UVW, C_OUTPUT_AMN, C_BOUNDARY_CHECK)
 
                 if(trim(adjustl(val)) == C_ON .or. trim(adjustl(toLowercase(val))) == 'true') then
                     itemp = SWITCHED_ON
@@ -1510,7 +1512,10 @@ integer function setNeptuneVar_char(        &
                     case(C_MANEUVERS)
                         call this%set_input(parName=C_MANEUVERS, val = toString(ltemp), set=.true.)
                         call this%derivatives_model%setPertSwitch(this%gravity_model, PERT_MANEUVERS, ltemp)
-
+                  !** CD/CR boundary check
+                  case(C_BOUNDARY_CHECK)
+                    call this%set_input(parName=C_BOUNDARY_CHECK, val = toString(ltemp), set=.true.)
+                    call this%satellite_model%setBoundaryCheck(ltemp)
                   !** covariance matrix propagation
                   case(C_COV_PROP)
                         call this%set_input(parName=C_COV_PROP, val=val, set=.true.)
diff --git a/src/propagatorschnittstelle.f90 b/src/propagatorschnittstelle.f90
index 99a0f84..f2530d5 100644
--- a/src/propagatorschnittstelle.f90
+++ b/src/propagatorschnittstelle.f90
@@ -100,6 +100,7 @@ subroutine OPI_Plugin_init(propagator) bind(c, name="OPI_Plugin_init")
   call OPI_Module_createProperty(propagator, "ecef_states_out", "OFF")
   call OPI_Module_createProperty(propagator, "mean_elements_out", "ON")
   call OPI_Module_createProperty(propagator, "return_set_matrix", "OFF")    ! instead of covariance
+  call OPI_Module_createProperty(propagator, "boundary_check", "ON")    ! Check strict boundaries for CR/CD
 
   do i = 1, 20
     !** manoeuvre --> should be temp, this is a little bit an overkill
@@ -686,6 +687,21 @@ function OPI_Plugin_propagate(propagator, data, julian_day, dt) result(opi_error
 
     end if
 
+    !** Strict checking of CD and CR coefficientes in Cannon Ball Mode:
+    !** on: CD>0, CR>0, <2.
+    !** off: CD>=0, CR unbounded
+    write(temp_string,*) OPI_Module_getPropertyString(propagator,"boundary_check")
+    !** set the object properties
+    ierr = neptune_instance%setNeptuneVar("BOUNDARY_CHECK", temp_string)
+    !** check error
+    if (ierr .ne. 0) then
+        slam_error = t%check_slam_error()
+        if (t%has_to_return()) return
+        if (slam_error) then
+            call resetError()
+        end if
+    endif
+
     !** Allow saving the states, if wanted
     store_data = .false.
     if (.not. create_cheby) then