! COLNEW_SOLVER, Release 1, a fortran 90 wrapper for the updated fortran 90 colnew software.
! The following colnew_wrapper module is hidden from the user program and calls the colnew software.
! The user program "uses" this module as an intermediate step in the call to the colnew module.
!------------------------------------------------------------------------------------------------------------


        MODULE COLNEW_WRAPPER

!-----------------------------------------------------------------------------------------------------------
!        Public routines inside this module include (for descriptions go to the beginning of each routine):
!   
!        COLNEW_INIT
!        COLNEW_SOLVER
!        COLNEW EVAL
!        COLNEW_TERMINATE
!        COLNEW_CONTINUE
!
!        Private routines inside this module include:
!
!        CHECK_STAT
!        COLNEW_LINSPACE
!
!        Private routines that are link and called in this module are:
!
!        COLNEW which is also link to the f77 blas routines.
!        APPSLN Colnew's solution evaluation routine at any point x.
!        
!------------------------------------------------------------------------------------------------------------
         
         USE COLNEW_M

      
         IMPLICIT NONE

!-----------------------------------------------------------------------------------------------------------

         PRIVATE ! All routines in this module are declared private except for the ones below which are
                 ! available to the main program.

         PUBLIC::COLNEW_INIT,COLNEW_SOLVER,COLNEW_TERMINATE, &
                 COLNEW_CONTINUE,COLNEW_EVAL

!----------------------------------------------------------------------------------------------------------

         ! Set up structured type colnew_sol that will store solution and associated information.
         TYPE,PUBLIC::COLNEW_SOL

              INTEGER::NODE ! Number of differential equations
              INTEGER::LEFTBC ! Number of left boundary conditions
              INTEGER::NSOL ! Number of solution derivative components equal to the sum of the orders of the
                            ! differential equations.
              INTEGER,DIMENSION(:),POINTER::ORDERS ! The orders of each differential equation.
              INTEGER::NPTS ! The number of mesh points
              INTEGER::NFIXPTS ! The number of additional fixed points besides the left and right endpoints.
              INTEGER::INFO ! Indicator of success or failure
              INTEGER::MXNSUB ! Maximum number of subintervals
              DOUBLE PRECISION, DIMENSION(:),POINTER::X ! The current mesh.
              INTEGER, DIMENSION(:),POINTER::ISPACE ! Integer work spaces needed for the colnew software
              DOUBLE PRECISION,DIMENSION(:),POINTER::FSPACE ! real work space needed for the colnew software.
    

         END TYPE COLNEW_SOL

!----------------------------------------------------------------------------------------------------------

         CONTAINS  ! The subroutine and functions to follow.

!----------------------------------------------------------------------------------------------------------

         SUBROUTINE COLNEW_EVAL(SOL,T,Z,FSUB)
            ! Evaluate the solution at the point X. The solution is returned in Z.
            ! All info required for the evaluation of the solution is in SOL. FSUB
            ! is required if the superinterpolant is used inside the colnew module.
            !--------------------------------------------------------------------------------------------
            !    CALLED BY: user's program
            !    CALLS: APPSLN (f90 program in Colnew module)
            !--------------------------------------------------------------------------------------------
            ! Input arguments
            TYPE(COLNEW_SOL)::SOL ! Current solution info
            DOUBLE PRECISION :: T ! Any point on the interval

            ! Output arguments 
            DOUBLE PRECISION :: Z(SOL%NSOL) ! Solution approximation
          

            EXTERNAL FSUB
            !-------------------------------------------------------------------------------------------

            CALL APPSLN(T,Z,SOL%FSPACE,SOL%ISPACE,FSUB)
                          

         END SUBROUTINE COLNEW_EVAL

!--------------------------------------------------------------------------------------------------------


         FUNCTION COLNEW_SOLVER(SOL_INIT,FSUB,&
                  BCSUB,GUESS, &
                  DFDY,DBCDY,NONLINEAR,NUM_COLL_POINTS,TRACE, &
                  PROBLEM_TYPE,FIXPOINTS, &
                  TTOL,MESH_HALVING)  RESULT(SOL)
                  ! This function is called by the user's main program to solve the boundary value problem.
                  ! It assumes that COLNEW_INIT has been called previously to form the solution structure
                  ! SOL_INIT. The ODEs and boundary conditions are made available through the user provided
                  ! subroutines,FSUB and BCSUB, and the corresponding partial derivatives are optionally 
                  ! available through the user provided subroutines, DFDY and DBCDY. If they are not 
                  ! available, divided differences are used in the routines DFDY_DIVDIFF and DBCDY_DIVDIFF.
                  ! COLNEW_SOLVER first sets up the call to the COLNEW software. The GUESS routine is also
                  ! optional. The COLNEW package then
                  ! returns an approximate solution in the Z vector of NSOL components where NSOL is the sum
                  ! of the orders of the original ODEs. 
                  ! Here is a description of the optional parameters in the argument list. The descriptions
                  ! can also be found in the main program:
                  !----------------------------------------------------------------------------------------
                  ! Other optional arguments with there keywords are included here. The default 
                  ! for each argument is indicated beside the corresponding value.
                  !
                  ! 1) nonlinear = 0  if the problem is linear 
                  !              = 1  if the problem is nonlinear (default)
                  ! 2) num_coll_points = number of collocation points per subinterval (=k)
                  !                      where max(orders).le.k.le 7. ???? ... ask why colnew 
                  !                      accepts k.lt.max(orders)?????
                  !                      k = 3 (default)
                  ! 3) trace = output control
                  !          = -1 for full diagnostic printout (default)
                  !          = 0 for selected printout
                  !          = 1 for no printout
                  ! 4) problem_type = 0 if the problem is regular (default)
                  !                 = 1 if the first relax factor is =rstart, and the nonlinear
                  !                     iteration does not rely on past convergence (use for only
                  !                     extra sensitive nonlinear problems)
                  !                 = 2 if we are to return immediately upon a) two successive
                  !                     nonconvergences, or b) after obtaining an error estimate 
                  !                     for the first time.
                  ! 5) fixpoints = an array of size "nfixpts" containing the points, other than
                  !                the left and right boundary points, which are to be included 
                  !                in every mesh.
                  !                fixpoints is an empty array by default
                  ! 6) ttol = tolerance applied to all components of the solution vector.
                  !           ttol=1.0d-5 is the default.
                  ! 7) mesh_halving = 1 if initial mesh is provided by the user and adaptive
                  !                   mesh selection is to be done (default)
                  !                 = 2 if initial mesh is provided by the user but in addition
                  !                     NO adaptive mesh selection is to be done.
                  !-----------------------------------------------------------------------------------------
                  ! CALLED BY: user program
                  ! CALLS : COLNEW (f90 program),CHECK_STAT
                  !----------------------------------------------------------------------------------------- 
                  ! Input Arguments (see descriptions above)
                  TYPE(COLNEW_SOL)::SOL_INIT

     

                  INTEGER,OPTIONAL::NONLINEAR
                  INTEGER,OPTIONAL::NUM_COLL_POINTS
                  INTEGER,OPTIONAL::TRACE
                  INTEGER,OPTIONAL::PROBLEM_TYPE
                  DOUBLE PRECISION,DIMENSION(:),OPTIONAL::FIXPOINTS
                  DOUBLE PRECISION,OPTIONAL::TTOL
                  INTEGER,OPTIONAL::MESH_HALVING
        
           
                  ! User supplied routines
                  EXTERNAL FSUB,BCSUB,DFDY,DBCDY,GUESS
                  OPTIONAL :: DFDY,DBCDY,GUESS

                  ! Output arguments

                  TYPE(COLNEW_SOL)::SOL !Updated solution information

                  ! Local variables (most of which needed for call to colnew)
                  INTEGER :: NCOMP  ! Number of ODEs
                  INTEGER :: M(SOL_INIT%NODE) ! Orders of the differential equations
                  INTEGER :: IPAR(11) ! 11 Components need for the call to colnew to be described 
                                      ! further below
                  INTEGER :: IFLAG ! Status of computation
                  INTEGER :: I ! Loop counter
                  INTEGER :: LTOL(SOL_INIT%NSOL) ! Tolerance label
                  DOUBLE PRECISION::ALEFT,ARIGHT ! Left and right end points
                  DOUBLE PRECISION::X ! Mesh ppoint
                  DOUBLE PRECISION::ZETA(SOL_INIT%NSOL) ! Location of boundary conditions: aright or aleft
                  DOUBLE PRECISION::TOL(SOL_INIT%NSOL) ! Tolerance
                  DOUBLE PRECISION,DIMENSION(SOL_INIT%NFIXPTS)::FIXPNT ! Array of extra fixed points
                  INTEGER::IER ! Dynamic array location error flag

                  ! Set up the arguments for the call to colnew.
                  ! Also, one can find  comments for the arguments in the colnew software.
                  !--------------------------------------------------------------------------------------




                  NCOMP=SOL_INIT%NODE
                  M=SOL_INIT%ORDERS
                  ALEFT=SOL_INIT%X(1)
                  ARIGHT=SOL_INIT%X(SOL_INIT%NPTS)

                  ! Store the location of each boundary point in the colnew argument zeta.
                  DO I=1,SOL_INIT%LEFTBC
                    ZETA(I)=ALEFT
                  ENDDO

                  DO I=(SOL_INIT%LEFTBC+1),SOL_INIT%NSOL
                    ZETA(I)=ARIGHT
                  ENDDO

                  IF(PRESENT(NONLINEAR))THEN
                    IPAR(1)=NONLINEAR
                  ELSE
                    IPAR(1)=1
                  ENDIF

                  IF(PRESENT(NUM_COLL_POINTS))THEN
                    IPAR(2)=NUM_COLL_POINTS
                  ELSE
                    IPAR(2)=3
                  ENDIF

                  ! Number of subintervals in the initial mesh.
                  IPAR(3)=(SOL_INIT%NPTS)-1

                  IPAR(4)=SOL_INIT%NSOL

                  IPAR(5)=SOL_INIT%MXNSUB*(4+3*SOL_INIT%NSOL+ & 
                         (5+IPAR(2)*SOL_INIT%NSOL)*  &
                         ((IPAR(2)+1)*SOL_INIT%NSOL)+4*SOL_INIT%NSOL**2)


                  IPAR(6)=SOL_INIT%MXNSUB*(3+(IPAR(2)+1)*SOL_INIT%NSOL)

                  ! Dynamically store memory for the colnew work arrays fspace and ispace 
                  ALLOCATE(SOL%FSPACE(3*IPAR(5)), &
                           SOL%ISPACE(3*IPAR(6)),STAT=IER)
                  CALL CHECK_STAT(IER)

                 IF(PRESENT(TRACE))THEN
                   IPAR(7)=TRACE
                 ELSE
                   IPAR(7)=-1
                 ENDIF

                 IF(PRESENT(MESH_HALVING))THEN
                   IPAR(8)=MESH_HALVING
                 ELSE
                   IPAR(8)=1
                 ENDIF

                 DO I=2,(SOL_INIT%NPTS)-1
                     SOL%FSPACE(I)=SOL_INIT%X(I)
                 ENDDO

                 IF(PRESENT(GUESS))THEN
                   IPAR(9)=1
                 ELSE
                   IPAR(9)=0
                 ENDIF
  
                 IF(PRESENT(PROBLEM_TYPE))THEN
                   IPAR(10)=PROBLEM_TYPE
                 ELSE
                   IPAR(10)=0
                 ENDIF

                 IF(PRESENT(FIXPOINTS))THEN
                   IPAR(11)=SOL_INIT%NFIXPTS
                   FIXPNT=FIXPOINTS
                 ELSE
                   IPAR(11)=0
                 ENDIF

                 DO I=1,SOL_INIT%NSOL
                   LTOL(I)=I
                 ENDDO
    
                 ! The tolerance level is the same for all solution and derivative values
                 IF(PRESENT(TTOL))THEN 
                   DO I=1,SOL_INIT%NSOL
                      TOL(I)=TTOL
                   ENDDO
                 ELSE
                   DO I=1,SOL_INIT%NSOL
                      TOL(I)=1.0D-5
                   ENDDO
                 ENDIF

    
                 ! Make the call to the Colnew subroutine inside the COLNEW module.

 
                   IF(PRESENT(DFDY))THEN
                      IF(PRESENT(DBCDY))THEN

                       CALL COLNEW(NCOMP,M,ALEFT,ARIGHT,ZETA,IPAR,LTOL,&
                            TOL,FIXPNT,SOL%ISPACE,SOL%FSPACE,IFLAG, &
                            FSUB,BCSUB,GUESS,DFSUB=DFDY,DGSUB=DBCDY)

                      ELSE

                       CALL COLNEW(NCOMP,M,ALEFT,ARIGHT,ZETA,IPAR,LTOL,&
                            TOL,FIXPNT,SOL%ISPACE,SOL%FSPACE,IFLAG, &
                            FSUB,BCSUB,GUESS,DFSUB=DFDY)

                      ENDIF

                    ELSE

                      IF(PRESENT(DBCDY))THEN

                       CALL COLNEW(NCOMP,M,ALEFT,ARIGHT,ZETA,IPAR,LTOL,&
                            TOL,FIXPNT,SOL%ISPACE,SOL%FSPACE,IFLAG,&
                            FSUB,BCSUB,GUESS,DGSUB=DBCDY)

                      ELSE

                       CALL COLNEW(NCOMP,M,ALEFT,ARIGHT,ZETA,IPAR,LTOL,&
                            TOL,FIXPNT,SOL%ISPACE,SOL%FSPACE,IFLAG,&
                            FSUB,BCSUB,GUESS)

                      ENDIF

                   ENDIF
                 



                 ! Load up the SOL fields for return back to the users main program.

                 SOL%NODE=SOL_INIT%NODE
                 SOL%LEFTBC=SOL_INIT%LEFTBC
                 SOL%NSOL=SOL_INIT%NSOL
                 SOL%NPTS=(SOL%ISPACE(1)) + 1
                 SOL%NFIXPTS=SOL_INIT%NFIXPTS
                 SOL%INFO=IFLAG
                 SOL%MXNSUB=SOL_INIT%MXNSUB
    

                 ALLOCATE(SOL%ORDERS(SOL_INIT%NODE), &
                          SOL%X(SOL%ISPACE(1)+1), &
                          STAT=IER)
                 CALL CHECK_STAT(IER)

                 SOL%ORDERS=SOL_INIT%ORDERS

                 DO I=1,SOL%ISPACE(1)+1
                  SOL%X(I)=SOL%FSPACE(I)
                 ENDDO

                

       END FUNCTION COLNEW_SOLVER

!-----------------------------------------------------------------------------------------------      


       FUNCTION COLNEW_INIT(NODE,LEFTBC,ORDERS,X,NFIXPTS, &
          MAX_NUM_SUBINTERVALS) RESULT(SOL)
          ! The major purpose of this subroutine is to initialize and store the initial mesh.
          ! The initial guess for the solution vector is not determined here. It is determined
          ! by making the subroutine GUESS one of the arguments in the call to colnew_solver.
          ! Other fields such as NODE, NSOL, LEFTBC, ORDERS, and NFIXPTS are loaded here.
          ! MAX_NUM_SUBINTERVALS is the maximum number of subintervals for the final mesh and an 
          ! optional argument. It can be set by the user.
          !--------------------------------------------------------------------------------------
          ! CALLED BY: User's program
          ! Calls to CHECK_STAT, COLNEW_LINSPACE
          !--------------------------------------------------------------------------------------
          ! Input arguments:         
          INTEGER :: NODE,LEFTBC
          DOUBLE PRECISION, DIMENSION(:) :: X ! An array containing the initial mesh
          INTEGER,DIMENSION(:)::ORDERS
          INTEGER,OPTIONAL::MAX_NUM_SUBINTERVALS,NFIXPTS

          ! Function output
          TYPE(COLNEW_SOL) :: SOL 

          !Local variables
          INTEGER :: NPTS,IER,NSOL ! Number of mesh points and memory allocation error flag. 
          LOGICAL :: DEFAULT_MESH ! Flag to signal default mesh.
          !----------------------------------------------------------------------------
       
           NSOL=SUM(ORDERS)

          !Set npts and check endpoints
          NPTS = SIZE(X)
          IF (X(1) >= X(NPTS)) THEN
            PRINT *,'THE INTERVAL [A,B] MUST HAVE A < B.'
            STOP
          END IF

          !   if npts = 2, use default of 10 equally spaced points.      
          DEFAULT_MESH = (NPTS == 2)
          IF (DEFAULT_MESH) NPTS = 6
    
          !   Create storage for allocatable fields of sol called x and orders.
          ALLOCATE(SOL%X(NPTS),SOL%ORDERS(NODE),STAT=IER)
          CALL CHECK_STAT(IER)

          !   Load x field
          IF (DEFAULT_MESH) THEN
            SOL%X = COLNEW_LINSPACE(X(1),X(2),NPTS)
          ELSE
            SOL%X = X
          END IF
 
          !   load other fields.
          SOL%NPTS = NPTS
          SOL%NODE = NODE
          SOL%NSOL = NSOL
          SOL%LEFTBC = LEFTBC
          SOL%INFO = 0
          SOL%ORDERS=ORDERS

          IF(PRESENT(MAX_NUM_SUBINTERVALS))THEN
             SOL%MXNSUB=MAX_NUM_SUBINTERVALS
          ELSE
             SOL%MXNSUB=50000
          ENDIF

          IF(PRESENT(NFIXPTS))THEN
             SOL%NFIXPTS=NFIXPTS
          ELSE
             SOL%NFIXPTS=0
          ENDIF

       END FUNCTION COLNEW_INIT
  
!--------------------------------------------------------------------------------------------------------

       SUBROUTINE CHECK_STAT(IER)
        ! Check status of an attempted array allocation or deallocation by
        ! checking error flag ier. print error condition and stop if 
        ! allocation/deallocation fails.
        !---------------------------------------------------------------------
        !   CALLED BY: COLNEW_SOLVER, COLNEW_INIT, COLNEW_TERMINATE
        !---------------------------------------------------------------------
        !   Input argument:
        INTEGER :: IER
        !---------------------------------------------------------------------
        IF (IER /= 0) THEN
          PRINT *,'A STORAGE ALLOCATION OR DEALLOCATION ERROR OCCURRED.'
          STOP
        END IF
        
       END SUBROUTINE CHECK_STAT  

!--------------------------------------------------------------------------------------------------------  

       FUNCTION COLNEW_LINSPACE(A,B,NPTS) RESULT(X)
        ! Generate a uniform mesh of npts points between a and b (inclusive);
        ! return mesh in array x, locations 1 to npts.
        !---------------------------------------------------------------------
        !   CALLED BY: COLNEW_INIT
        !---------------------------------------------------------------------
        !   Input arguments: 
        INTEGER :: NPTS
        DOUBLE PRECISION :: A,B
    
        ! Function result
        DOUBLE PRECISION ::  X(NPTS)

        ! Local variables
        INTEGER :: K ! Loop index
        DOUBLE PRECISION :: DELTA ! Subinterval size
        INTEGER :: MNPTS 
        !----------------------------------------------------------------------
        !   Mesh must have at least two points.
        MNPTS = MAX(NPTS,2)

        DELTA = (B - A) / (MNPTS - 1)
        X = (/ (A + DELTA*(K-1), K = 1,MNPTS) /)

       END FUNCTION COLNEW_LINSPACE

!---------------------------------------------------------------------------------------------------------

       SUBROUTINE COLNEW_TERMINATE(SOL)
        ! Colnew_terminate deallocates the allocateable array fields of sol, namely
        ! x, ispace, fspace, orders
        !--------------------------------------------------------------------------
        !     CALLED BY: USER'S PROGRAM
        !     CALLS: CHECK_STAT
        !--------------------------------------------------------------------------
        !   input/output argument:
        TYPE(COLNEW_SOL) :: SOL

        !   LOCAL VARIABLES:    
        INTEGER :: IER
        !----------------------------------------------------------------------------


        DEALLOCATE(SOL%ISPACE,STAT=IER)
!
        CALL CHECK_STAT(IER)
! 
        DEALLOCATE(SOL%FSPACE,STAT=IER)
!
        CALL CHECK_STAT(IER)
!
        DEALLOCATE(SOL%ORDERS,STAT=IER)
!!
        CALL CHECK_STAT(IER)
 
        DEALLOCATE(SOL%X,STAT=IER)

        CALL CHECK_STAT(IER)

       END SUBROUTINE COLNEW_TERMINATE

!----------------------------------------------------------------------------------------------------------


         FUNCTION COLNEW_CONTINUE(SOL_INIT,FSUB,&
                  BCSUB, &
                  DFDY,DBCDY,GUESS,NONLINEAR,NUM_COLL_POINTS,TRACE, &
                  PROBLEM_TYPE,FIXPOINTS, &
                  TTOL,MESH_HALVING)  RESULT(SOL)
                  ! This routine is used to carry out simple continuation after colnew is first called inside
                  ! colnew_solver. This function is called by the user's main program and must have the same
                  ! set of optional parameters that were used in the original call to colnew_solver. 
                  ! This function has almost the except same code as the colnew_solver routine.
                  ! Two important differences are that ipar(3) equals the number of subintervals of the previous
                  ! mesh and ipar(9) is set to 2 to "remind" colnew that the previous mesh and initial guess are
                  ! already stored in fspace.
                  ! Here is a description of the optional parameters in the argument list. The descriptions
                  ! can also be found in the main program:
                  !----------------------------------------------------------------------------------------
                  ! Other optional arguments with there keywords are included here. The default 
                  ! for each argument is indicated beside the corresponding value.
                  !
                  ! 1) nonlinear = 0  if the problem is linear 
                  !              = 1  if the problem is nonlinear (default)
                  ! 2) num_coll_points = number of collocation points per subinterval (=k)
                  !                      where max(orders).le.k.le 7. ???? ... ask why colnew 
                  !                      accepts k.lt.max(orders)?????
                  !                      k = 3 (default)
                  ! 3) trace = output control
                  !          = -1 for full diagnostic printout (default)
                  !          = 0 for selected printout
                  !          = 1 for no printout
                  ! 4) problem_type = 0 if the problem is regular (default)
                  !                 = 1 if the first relax factor is =rstart, and the nonlinear
                  !                     iteration does not rely on past convergence (use for only
                  !                     extra sensitive nonlinear problems)
                  !                 = 2 if we are to return immediately upon a) two successive
                  !                     nonconvergences, or b) after obtaining an error estimate 
                  !                     for the first time.
                  ! 5) fixpoints = an array of size "nfixpts" containing the points, other than
                  !                the left and right boundary points, which are to be included 
                  !                in every mesh.
                  !                fixpoints is an empty array by default
                  ! 6) ttol = tolerance applied to all components of the solution vector.
                  !           ttol=1.0d-5 is the default.
                  ! 7) mesh_halving = 1 if initial mesh is provided by the user and adaptive
                  !                   mesh selection is to be done (default)
                  !                 = 2 if initial mesh is provided by the user but in addition
                  !                     NO adaptive mesh selection is to be done.
                  !-----------------------------------------------------------------------------------------
                  ! CALLED BY: user program
                  ! CALLS : COLNEW (f90 program),CHECK_STAT
                  !----------------------------------------------------------------------------------------- 
                  ! Input Arguments (see descriptions above)
                  TYPE(COLNEW_SOL)::SOL_INIT
                  INTEGER,OPTIONAL::NONLINEAR
                  INTEGER,OPTIONAL::NUM_COLL_POINTS
                  INTEGER,OPTIONAL::TRACE
                  INTEGER,OPTIONAL::PROBLEM_TYPE
                  DOUBLE PRECISION,DIMENSION(:),OPTIONAL::FIXPOINTS
                  DOUBLE PRECISION,OPTIONAL::TTOL
                  INTEGER,OPTIONAL::MESH_HALVING
        

                  DOUBLE PRECISION :: HR1

                  INTEGER :: HR2,HR3


                  real :: sum_cpu



           
                  ! User supplied routines
                  EXTERNAL FSUB,BCSUB,DFDY,DBCDY,GUESS
                  OPTIONAL :: DFDY,DBCDY,GUESS

                  ! Output arguments

                  TYPE(COLNEW_SOL)::SOL !Updated solution information

                  ! Local variables (most of which needed for call to colnew)
                  INTEGER :: NCOMP  ! Number of ODEs
                  INTEGER :: M(SOL_INIT%NODE) ! Orders of the differential equations
                  INTEGER :: IPAR(11) ! 11 Components need for the call to colnew to be described 
                                      ! further below
                  INTEGER :: IFLAG ! Status of computation
                  INTEGER :: I ! Loop counter
                  INTEGER :: LTOL(SOL_INIT%NSOL) ! Tolerance label
                  DOUBLE PRECISION::ALEFT,ARIGHT ! Left and right end points
                  DOUBLE PRECISION::X ! Mesh ppoint
                  DOUBLE PRECISION::ZETA(SOL_INIT%NSOL) ! Location of boundary conditions: aright or aleft
                  DOUBLE PRECISION::TOL(SOL_INIT%NSOL) ! Tolerance
                  DOUBLE PRECISION,DIMENSION(SOL_INIT%NFIXPTS)::FIXPNT ! Array of extra fixed points
                  INTEGER::IER ! Dynamic array location error flag

                  ! Set up the arguments for the call to colnew.
                  ! Also, one can find  comments for the arguments in the colnew software.
                  !--------------------------------------------------------------------------------------


                  NCOMP=SOL_INIT%NODE
                  M=SOL_INIT%ORDERS
                  ALEFT=SOL_INIT%X(1)
                  ARIGHT=SOL_INIT%X(SOL_INIT%NPTS)


                  ! Store the location of each boundary point in the colnew argument zeta.
                  DO I=1,SOL_INIT%LEFTBC
                    ZETA(I)=ALEFT
                  ENDDO

                  DO I=(SOL_INIT%LEFTBC+1),SOL_INIT%NSOL
                    ZETA(I)=ARIGHT
                  ENDDO

                  IF(PRESENT(NONLINEAR))THEN
                    IPAR(1)=NONLINEAR
                  ELSE
                    IPAR(1)=1
                  ENDIF

                  IF(PRESENT(NUM_COLL_POINTS))THEN
                    IPAR(2)=NUM_COLL_POINTS
                  ELSE
                    IPAR(2)=3
                  ENDIF

                  ! Number of subintervals in the initial mesh.
                  IPAR(3)=SOL_INIT%ISPACE(1)

                  IPAR(4)=SOL_INIT%NSOL

                  IPAR(5)=SOL_INIT%MXNSUB*(4+3*SOL_INIT%NSOL+ & 
                         (5+IPAR(2)*SOL_INIT%NSOL)*  &
                         ((IPAR(2)+1)*SOL_INIT%NSOL)+4*SOL_INIT%NSOL**2)

                  IPAR(6)=SOL_INIT%MXNSUB*(3+(IPAR(2)+1)*SOL_INIT%NSOL)

                 IF(PRESENT(TRACE))THEN
                   IPAR(7)=TRACE
                 ELSE
                   IPAR(7)=-1
                 ENDIF

                 IF(PRESENT(MESH_HALVING))THEN
                   IPAR(8)=MESH_HALVING
                 ELSE
                   IPAR(8)=1
                 ENDIF

                 IPAR(9)=2
  
                 IF(PRESENT(PROBLEM_TYPE))THEN
                   IPAR(10)=PROBLEM_TYPE
                 ELSE
                   IPAR(10)=0
                 ENDIF

                 IF(PRESENT(FIXPOINTS))THEN
                   IPAR(11)=SOL_INIT%NFIXPTS
                   FIXPNT=FIXPOINTS
                 ELSE
                   IPAR(11)=0
                 ENDIF

                 DO I=1,SOL_INIT%NSOL
                   LTOL(I)=I
                 ENDDO
    
                 ! The tolerance level is the same for all solution and derivative values
                 IF(PRESENT(TTOL))THEN 
                   DO I=1,SOL_INIT%NSOL
                      TOL(I)=TTOL
                   ENDDO
                 ELSE
                   DO I=1,SOL_INIT%NSOL
                      TOL(I)=1.0D-5
                   ENDDO
                 ENDIF
   

                   IF(PRESENT(DFDY))THEN
                      IF(PRESENT(DBCDY))THEN

                       CALL COLNEW(NCOMP,M,ALEFT,ARIGHT,ZETA,IPAR,LTOL,&
                            TOL,FIXPNT,SOL_INIT%ISPACE,SOL_INIT%FSPACE,&
                            IFLAG,FSUB,BCSUB,GUESS,DFSUB=DFDY, &
                            DGSUB=DBCDY)

                      ELSE

                       CALL COLNEW(NCOMP,M,ALEFT,ARIGHT,ZETA,IPAR,LTOL,&
                            TOL,FIXPNT,SOL_INIT%ISPACE,SOL_INIT%FSPACE,&
                            IFLAG,FSUB,BCSUB,GUESS,DFSUB=DFDY)

                      ENDIF

                    ELSE

                      IF(PRESENT(DBCDY))THEN


                       CALL COLNEW(NCOMP,M,ALEFT,ARIGHT,ZETA,IPAR,LTOL,&
                            TOL,FIXPNT,SOL_INIT%ISPACE,SOL_INIT%FSPACE,&
                            IFLAG,FSUB,BCSUB,GUESS,DGSUB=DBCDY)


                      ELSE


                       CALL COLNEW(NCOMP,M,ALEFT,ARIGHT,ZETA,IPAR,LTOL,&
                            TOL,FIXPNT,SOL_INIT%ISPACE,SOL_INIT%FSPACE,&
                            IFLAG,FSUB,BCSUB,GUESS)

                      ENDIF

                   ENDIF


                 ! Load up the SOL fields for return back to the users main program.


                 SOL%NODE=SOL_INIT%NODE
                 SOL%LEFTBC=SOL_INIT%LEFTBC
                 SOL%NSOL=SOL_INIT%NSOL
                 SOL%NPTS=(SOL_INIT%ISPACE(1)) + 1
                 SOL%NFIXPTS=SOL_INIT%NFIXPTS
                 SOL%INFO=IFLAG
                 SOL%MXNSUB=SOL_INIT%MXNSUB
   

                 ALLOCATE(SOL%ORDERS(SOL_INIT%NODE), &
                          SOL%X((SOL_INIT%ISPACE(1))+1), &
                          STAT=IER)
                 CALL CHECK_STAT(IER)

                 SOL%ORDERS=SOL_INIT%ORDERS

                 DO I=1,SOL_INIT%ISPACE(1)+1
                  SOL%X(I)=SOL_INIT%FSPACE(I)
                 ENDDO

         
                  ! Dynamically store memory for the colnew work arrays fspace and ispace 
                  ALLOCATE(SOL%FSPACE(IPAR(5)),SOL%ISPACE(IPAR(6)), &
                           STAT=IER)
                  CALL CHECK_STAT(IER)

                  SOL%FSPACE=SOL_INIT%FSPACE
                  SOL%ISPACE=SOL_INIT%ISPACE


       END FUNCTION COLNEW_CONTINUE
!-----------------------------------------------------------------------------------------------      
       END MODULE COLNEW_WRAPPER
