anchor.diff:                                  index_s_prime = 1 ! Case: no state variable is used --> dummy index
                               end if
                               call egm_maximize(exp_valuef(:,index_s_prime,index_z),exp_dvaluef(:,index_s_prime,index_z),d,0.0_dp,z,m_end,c_end,valuef_end,k)
  +
                           end if
   
                           do index_s=1,glen_s
                               s=model%s_grid(index_s)
  -                            ! The FOC depend on the state; therefore, it has to be evaluated for every state
  +                            ! The FOC depends on the state; therefore, it has to be evaluated for every state
                               if (.NOT.model%state_independent_foc) then
                                   if (associated(model%Psi)) then
                                       index_s_prime = model%Psi(index_s,index_d,index_z)
  
positive.code:                         do index_a=1,glen_a
                              d_max = d_index_choice(index_a,index_s,index_z)
                              policy_c(index_a,index_s,index_z) = policy_c_con(index_a,d_max,index_s,index_z)
                              policy_d=model%d_grid(d_max)
                          end do
                      end do
                  end do
              end if
  
          end subroutine
  
negative.code:         end function
  
  
  
  
          ! *******************************************************************************************************
          ! ** Private helper subroutines
          ! *******************************************************************************************************
          subroutine requireSolved()
              if (.NOT.is_solved) then
                  print *, "Error: The model must first be solved before obtaining policy and value functions."
                  stop
              end if
          end subroutine
  
          logical function validate_model(m) result(valid)
              type(ncegm_model), intent(in) :: m
              valid = .FALSE.
              if (.NOT. (associated(m%F) .AND. associated(m%dF) .AND. associated(m%d2F))) then
                  write (error_unit,*) "Invalid model: The return function and first and second derivatives with respect to c must be provided."
              elseif (.NOT. associated(m%Gamma)) then
                  write (error_unit,*) "Invalid model: The budget constraint function Gamma must be provided"
              elseif (.NOT. (allocated(m%a_grid) .AND. allocated(m%d_grid))) then
                  write (error_unit,*) "Invalid model: The grids for A and D must be provided."
              else if (allocated(m%s_grid) .NEQV. associated(m%Psi)) then
                  write (error_unit,*) "Invalid model: If the state variable s is used, s_grid and the transition function Psi must be provided."
              else if (allocated(m%z_grid) .NEQV. allocated(m%z_transition)) then
                  write (error_unit,*) "Invalid model: If the stochastic state variable z is used, z_grid and the transition matrix must be provided."
              elseif (.NOT. allocated(m%V_initial)) then
                  write (error_unit,*) "Invalid model: Missing initial guess for value function.s"
              elseif (m%beta<=0 .OR. m%beta>=1) then
                  write (error_unit,*) "Invalid model: Discount factor Beta must be in interval (0,1)"
              else
                  valid = .TRUE.
              end if
          end function
  
  
          ! *******************************************************************************************************
          ! ** Public getter functions for the module client to obtain the output of the algorithm
          ! *******************************************************************************************************
          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_getPolicy_aprime() result(aprime)
              real(dp), dimension(glen_a,glen_s,glen_z) :: aprime
              call requireSolved
              aprime = policy_aprime
          end function
  
          function ncegm_getValueFunction() result(v)
              real(dp), dimension(glen_a,glen_s,glen_z) :: v
              call requireSolved
              v = valuef
          end function
  
          function ncegm_getConditionalPolicy_c() result(c)
              real(dp), dimension(glen_a,glen_d,glen_s,glen_z) :: c
              call requireSolved
              c = policy_c_con
          end function
  
          function ncegm_getConditionalPolicy_aprime() result(aprime)
              real(dp), dimension(glen_a,glen_d,glen_s,glen_z) :: aprime
              call requireSolved
              aprime = policy_aprime_con
          end function
  
          function ncegm_getConditionalValueFunction() result(v)
              real(dp), dimension(glen_a,glen_d,glen_s,glen_z) :: v
              call requireSolved
              v = valuef_next_con
          end function
  
          ! *******************************************************************************************************
          ! ** Convenience functions to simplify access to array- and scalar-valued model functions
          ! *******************************************************************************************************
  
          function F_sc(c,d,s,z)
              real(dp), intent(in)   :: c,d,s,z
              real(dp)               :: F_sc
  