anchor.diff:                              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)
  +                                c = policy_c_con(index_a,d_max,index_s,index_z)
  +                                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(c,d,s,z)
                               end if
                           end do
                           if (.NOT. associated(model%dGamma)) then
  
positive.code:                 epsilon = EPSILON_DEFAULT
              end if
  
              ! Now solve the model
              call valfniteration
              is_solved = .TRUE.
          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()
              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
              valuef = model%V_initial
  
              do index_z=1, glen_z
  
negative.code:     implicit none
      private
      public :: ncegm_setup,ncegm_model,ncegm_solve
      public :: ncegm_getPolicy_c,ncegm_getPolicy_aprime,ncegm_getPolicy_d,ncegm_getValueFunction
      public :: ncegm_getConditionalPolicy_c,ncegm_getConditionalPolicy_aprime,ncegm_getConditionalValueFunction
      save
  
      ! *******************************************************************************************************
      ! ** Types and interfaces
      ! *******************************************************************************************************
      interface
          function ReturnFunction(c,d,s,z) result(r)
              use kinds, only: dp
  