module blendmod
!$$$   module documentation block
!                .      .    .                                       .
! module:    blendmod        blend function module
!   prgmmr: parrish          org: np22                date: 2012-11-26
!
! abstract: contains routines needed for polynomial blending function.  This
!            function is defined over the interval 0<=x<=1, with value 0 at
!            x=0, and value 1 at x=1. The order parameter iord defines the
!            order of continuity of the function at the endpoints.
!            For n=0, the blend function = x for 0 <= x <=1, and =1 for x>1, =0 for x<0.
!            In this case, the function is continuous, but the derivatives are discontinuous.
!            For n>0, the function and its first n derivatives are continuous at x=0 and x=1.
!            See the comments in subroutine blend for more details.  Subroutine
!            blend was originally written by Jim Purser.
!
! program history log:
!   2012-11-26  parrish
!
! subroutines included:
!  init_blend
!  blend_f
!  blend_df
!  blend

   use kinds, only: r_kind,i_kind

   implicit none

! set default to private
   private
! set subroutines to public
   public :: init_blend
   public :: blend_f
   public :: blend_df
   public :: blend
! set passed variables to public

   integer(i_kind) mm,iblend(0:40)
   real(r_kind) x0,dxx,dnorm
  
   contains

   subroutine init_blend(xbegin,xend,iord,ierror)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_blend
!   prgmmr: parrish          org: np22                date: 2012-11-26
!
! abstract: initialize private module variables so routines blend_f and blend_df generate the
!            appropriate blend function and 1st derivative for blend function order iord
!            and transition range xbegin to xend.  Note that xbegin can be greater 
!            than xend, but xbegin==xend is not allowed.
!
! program history log:
!   2012-11-26  parrish - initial documentation
!
!   input argument list:
!     xbegin - beginning of blend function interval (  blend_f(xbegin)=0 )
!     xend   - end of blend function interval ( blend_f(xend)=1 )
!     iord   - order of blend function
!     ierror - error return:  =0 for successful return, /=0 only if xbegin==xend.
!               NOTE:  xbegin > xend is allowed.
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

     use constants, only: half,one
     implicit none

     integer(i_kind),intent(in):: iord
     real(r_kind),intent(in):: xbegin,xend
     integer(i_kind),intent(inout):: ierror

     real(r_kind) x1,xmid,y,dy

!      ierror = 0 for normal return, /=0 if xend=xbegin

     if(xend==xbegin) then
        ierror=1
        return
     end if
     mm=iord
     x0=xbegin
     x1=xend
     dxx=one/(x1-x0)
     call blend(mm,iblend)
     xmid=x0+half*(x1-x0)
     dnorm=one
     call blend_df(xmid,y,dy)
     dnorm=one/dy

   end subroutine init_blend

   subroutine blend_f(x_in,y)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    blend_f
!   prgmmr: parrish          org: np22                date: 2012-11-26
!
! abstract:  for input value x_in, return blend function value y, where 0 <=y < =1.
!           Before calling blend_f, make sure that init_blend has been called
!           for the desired values of x where the blend function transitions
!           from 0 to 1.
!
! program history log:
!   2012-11-26  parrish - initial documentation
!
!   input argument list:
!     x_in   - input argument value
!
!   output argument list:
!     y      - output value  blend_f(x_in)
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

     use constants, only: zero,one
     implicit none

     real(r_kind),intent(in):: x_in
     real(r_kind),intent(out):: y
     real(r_kind) x
     integer(i_kind) j

     x=(x_in-x0)*dxx
     if(x < zero) then
        y=zero
     elseif(x > one) then
        y=one
     else
        y=iblend(mm)
        do j=mm-1,0,-1
           y=x*y+iblend(j)
        end do
        y=y*x**(mm+1)
     end if

   end subroutine blend_f

   subroutine blend_df(x_in,y,dy)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    blend_df
!   prgmmr: parrish          org: np22                date: 2012-11-26
!
! abstract:  for input value x_in, return blend function in y and 1st derivative
!             of blend function in dy.
!           Before calling blend_df, make sure that init_blend has been called
!           for the desired values of x where the blend function transitions
!           from 0 to 1.
!
! program history log:
!   2012-11-26  parrish - initial documentation
!
!   input argument list:
!     x_in   - input argument value
!
!   output argument list:
!     y      - output value of blend function for input value x_in.
!     dy     - output value of derivative of blend function for input value x_in.
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

     use constants, only: zero,one
     implicit none

     real(r_kind),intent(in):: x_in
     real(r_kind),intent(out):: y,dy
     real(r_kind) x
     integer(i_kind) j

     x=(x_in-x0)*dxx
     if(x < zero) then
        y=zero ; dy=zero
     elseif(x > one) then
        y=one ; dy=zero
     else
        dy=zero
        y=iblend(mm)
        do j=mm-1,0,-1
           dy=x*dy+y
           y=x*y+iblend(j)
        end do
        dy=((mm+one)*y+x*dy)*x**mm
        y=y*x**(mm+1)
     end if
     dy=dy*dnorm

end subroutine blend_df

subroutine blend(n,iblend)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    blend
!   prgmmr: purser           org: w/nmc22     date:1998
!
! abstract: put coefficients for n+1,..,2n+1, into iblend(0),..
!           iblend(n)
!
! program history log:
!   2004-05-13  kleist  documentation
!   2008-04-23  safford - rm unused uses
!   2012-11-26  parrish - move from prewgt.f90 to this module.
!
!   input argument list:
!     n      - number of powers to blend
!
!   output argument list:
!     iblend - blended coefficients
!
! remarks: put the coefficients for powers n+1,..,2n+1, into iblend(0),
!          ..iblend(n),for the "blending polynomial" of continuity-
!          degree n in the interval [0,1].  For example, with n=1, the 
!          blending polynomial has up to 1st derivatives continuous 
!          with y(0)=0, y(1)=1, y'(0)=y'(1)=0, when y(x)=3x^2-2x^3. 
!          Hence iblend={3,-2}
! 
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
  use kinds, only: i_kind
  implicit none

! Declare passed variables
  integer(i_kind)               ,intent(in   ) :: n
  integer(i_kind),dimension(0:n),intent(  out) :: iblend

! Declare local parameters
  integer(i_kind),parameter:: nn=12

! Declare local variables
  integer(i_kind) np,i,j,ib
  integer(i_kind),dimension(0:nn):: ipascal(0:nn)

  if(n>nn)stop
  np=n+1
  do i=0,n
    ipascal(i)=0
  enddo

  ipascal(0)=1
  do i=0,n
     do j=i,1,-1
        ipascal(j)=ipascal(j)-ipascal(j-1)
     enddo
  enddo

  ib=1
  do i=1,n
     ib=(ib*2*(2*i+1))/i
  enddo
  do j=0,n
     iblend(j)=(ib*ipascal(j))/(np+j)
  enddo

  return
end subroutine blend

end module blendmod