subroutine msplin(nbreak,nmax,npts,dx,ca, &
     nbc1,alpha, &
     nbcn,beta, &
     lunerr,ierr)
  use iso_c_binding, only: fp => c_double
  !
  !
  !       bob mccann, 16-jul-87
  !
  !
  !       --------
  !       comments
  !       --------
  !
  !
  !	the input variables are the same as for spl1dr_nb except for the
  !   addition of nbreak which is the index of the knot which is to be
  !   multiply defined in order to "break" the spline and thus allow the
  !   first and higher order derivatives to become discontinuous there.
  !
  !
  !       -------------------
  !       global declarations
  !       -------------------
  !
  !
  !       include 'msplin.blk/nolist' !communicate with test program
  !
  !
  !       ---------------
  !       input variables
  !       ---------------
  !
  !
  !============
  implicit none
  integer :: nbreak,ibcn,ibc1
  !============
  real(fp) :: zlope0,delta,f0,f0pj,f0p2j,f0p3j,zlope1,f1,f1mj,f1m2j
  real(fp) :: f1m3j,zfj,zfjm1,zfjm2,zfjm3,zdflft,zfjp1,zfjp2,zfjp3
  real(fp) :: zdfrt
  !============
  integer ::	nmax		!the second dim. of c
  integer ::	npts		!the number of functional values
  !
  real(fp) ::	dx		!the constant step size
  real(fp) ::	ca(4,nmax)	!ca(1,j) are specified
  !				!ca(2-4,j) are returned
  integer ::	nbc1, nbcn	!left and right boundary
  !					!condition specifiers
  !					!=(first,second,or third)
  real(fp) ::	alpha, beta	!specified left or right side
  !				!derivative
  !
  integer ::	lunerr		!unit number for error messages
  !
  !
  !       ----------------
  !       output variables
  !       ----------------
  !
  !
  !	ca(2-4,j) are returned with a discontinuity at j=nbreak
  !
  integer ::	ierr		!=0 => ok, error otherwise
  !
  !
  !       ------------------
  !       local declarations
  !       ------------------
  !
  !
  !***
  !
  !
  !       ----------------------
  !       0.1     initialization
  !       ----------------------
  !
  !
  if(nbreak+3.gt.npts .or. nbreak-3.lt.1) then
     ierr = 1
     write(lunerr,80101) npts, nbreak
     return
  endif
  !
  !
  !       --------------------------------------------------------------------
  !       1.0     compute the boundary conditions at the break (multiple node)
  !       --------------------------------------------------------------------
  !
  zfj   = ca(1,nbreak)
  zfjm1 = ca(1,nbreak-1)
  zfjm2 = ca(1,nbreak-2)
  zfjm3 = ca(1,nbreak-3)
  !
  zdflft = zlope1(dx,zfj,zfjm1,zfjm2,zfjm3)
  !
  zfjp1 = ca(1,nbreak+1)
  zfjp2 = ca(1,nbreak+2)
  zfjp3 = ca(1,nbreak+3)
  !
  zdfrt = zlope0(dx,zfj,zfjp1,zfjp2,zfjp3)
  !
  !
  !       -----------------------------------------
  !       2.0     make the appropriate splidr calls
  !       -----------------------------------------
  !
  !
  !
  ibcn = 1    !specify first derivative b.c. on right
  call spl1dr_nb(nmax,nbreak,dx,ca, &
       nbc1,alpha, &  !the user specifies this b.c. &
       ibcn,zdflft,  &!this one in computed here &
       lunerr,ierr)
  !
  ibc1 = 1    !specify first derivative b.c. on left
  call spl1dr_nb(npts-nbreak+1,npts-nbreak+1,dx,ca(1:4,nbreak:npts), &
       ibc1,zdfrt,&   !this one is computed here &
       nbcn,beta,&    !the user specifies this b.c. &
       lunerr,ierr)
  !
  return
80101 format(' msplin error: nbreak within 3 of edge. npts=',i5, &
           ' nbreak=',i5)
end subroutine msplin
