anchor.diff:                  epsilon = EPSILON_DEFAULT
               end if
   
  +            if (present(info_message_inverval_in)) then
  +                info_message_inverval = info_message_inverval_in
  +            else
  +                info_message_inverval = STATUS_PRINT_INTERVAL_DEFAULT
  +            end if
  +
               ! Now solve the model
  -            call valfniteration
  -            is_solved = .TRUE.
  +            call valfniteration(is_solved)
  +            if (present(status)) status = is_solved
  +
           end subroutine ncegm_solve
   
  -        ! Public module getter methods
  -        function ncegm_getPolicy_d() result(d)
  -            real(dp), dimension(glen_a,glen_s,glen_z) :: d
  -            call requireSolved
  -            d = policy_d
  -        end function
  -
  -        function ncegm_getPolicy_c() result(c)
  -            real(dp), dimension(glen_a,glen_s,glen_z) :: c
  -            call requireSolved
  -            c = policy_c
  -        end function
  -
  -        function ncegm_getValueFunction() result(v)
  -            real(dp), dimension(glen_a,glen_s,glen_z) :: v
  -            call requireSolved
  -            v = valuef
  -        end function
  -
  -        ! EGM algorithm functions
  -
  -        subroutine valfniteration()
  +        ! *******************************************************************************************************
  +        ! ** EGM core subroutines
  +        ! *******************************************************************************************************
  +
  +        subroutine valfniteration(status)
               integer                                          :: index_a,index_d,index_s,index_z,index_s_prime,iter
               real(dp), dimension(glen_a,glen_s,glen_z)        :: valuef_next,dvaluef,dvaluef_next,exp_valuef,exp_dvaluef
  -            real(dp), dimension(glen_a,glen_d,glen_s,glen_z) :: policy_c_con, policy_a_prime_con, valuef_next_con
               integer, dimension(glen_a,glen_s,glen_z)         :: d_index_choice
               real(dp)                                         :: sup_norm_diff
               real(dp), dimension(glen_a)                      :: m_end, c_end, valuef_end
  -            integer                                          :: k,d_max
  -            real(dp)                                         :: d,s,z,t1,t2
  -
  -            ! TODO: set and handle policy_a_prime_con
  -
  -            ! 1. Compute derivative of the initial value function and m-grid
  +            integer                                          :: k,d_max,i1
  +            real(dp)                                         :: d,s,z,c,t1,t2
  +            logical, intent(out)                             :: status
  +
  +            ! 1. Compute derivative of the initial value function and the m-grid
               valuef = model%V_initial
   
               do index_z=1, glen_z
  
positive.code:                             d = model%d_grid(d_max)
                              valuef_next(index_a,index_s,index_z) = valuef_next_con(index_a,d_max,index_s,index_z)
                              if (associated(model%dGamma)) then
                                  dvaluef_next(index_a,index_s,index_z) = dGamma_sc(model%a_grid(index_a),model%d_grid(d_max),model%s_grid(index_s),model%z_grid(index_z))*dF(policy_c_con(index_a,d_max,index_s,index_z),d,s,z)
                              end if
                          end do
                          if (.NOT. associated(model%dGamma)) then
  
negative.code:             glen_s = size(model%s_grid)
              glen_z = size(model%z_grid)
              glen_d = size(model%d_grid)
  
              ! (Re)allocate all arrays
              if (allocated(valuef)) deallocate(valuef, policy_c, policy_aprime, policy_d, m_grid,policy_c_con, policy_aprime_con, valuef_next_con)
              allocate(valuef(glen_a,glen_s,glen_z), policy_c(glen_a,glen_s,glen_z), policy_aprime(glen_a,glen_s,glen_z), policy_d(glen_a,glen_s,glen_z), m_grid(glen_a,glen_d,glen_s,glen_z))
              allocate(policy_c_con(glen_a,glen_d,glen_s,glen_z), policy_aprime_con(glen_a,glen_d,glen_s,glen_z), valuef_next_con(glen_a,glen_d,glen_s,glen_z))
  
              ! Ready to solve model
              is_setup = .TRUE.
              status = .TRUE.
              is_solved = .FALSE.
          end subroutine ncegm_setup
  
          ! **
          ! * Starts the solution algorithm of the model set up by a previous call to ncegm_setup(). The algorithm terminates if a sufficient convergence is achieved
          ! * or if the maximum number of iterations is achieved. The input parameters allow to change parameters of the algorithm. The output parameter status can be used
          ! * to check if the value functions converged. All parameters are optional.
          ! *
          ! * Input:
          ! *     - max_iter_in (optional): The maximum number of value function iterations. If the value functions have not converged, the algorithm will stop after max_iter_in iterations.
          ! *     - epsilon_in (optional): The supremum norm of the difference of two subsequent value functions below which the convergence is deemed sufficient and the algorithm is stopped.
          ! *     - info_message_interval_in (optional): This value specifys the amount of iterations after which a status is printed to the output showing the current convergence.
          ! *
          ! * Output:
          ! *     - status (optional): a flag indicating if the value functions have converged
          ! **
          subroutine ncegm_solve(max_iter_in, epsilon_in, info_message_inverval_in, status)
              integer, intent(in), optional          :: max_iter_in,info_message_inverval_in
              real(dp), intent(in), optional         :: epsilon_in
              logical, intent(out), optional         :: status
  
              if (.NOT.is_setup) return
  
  