cos_rpdederive_interface.f90

Go to the documentation of this file.

module cos_rpdederive_interface
  interface
  subroutine cos_rpdederive(dydx,K,x,y,c0,c1,d,o)
    use cos_precision
    implicit none
    integer,intent(in)          :: k
    real(kind=dp),intent(inout) :: dydx(k)
    real(kind=dp),intent(in)    :: x(k),y(k)
    integer,intent(in)          :: c0,c1,d,o
    real(kind=dp),dimension(K)  :: dx,dy,k0,k1
    real(kind=dp)               :: y0,y1,x1,x0
  end subroutine cos_rpdederive
  end interface
end module cos_rpdederive_interface
!*******************************************************************************
subroutine cos_rpdederive(dydx,K,x,y,c0,c1,d,o)

  use cos_precision
  implicit none
  integer,intent(in)          :: k
  real(kind=dp),intent(inout) :: dydx(k)
  real(kind=dp),intent(in)    :: x(k),y(k)
  integer,intent(in)          :: c0,c1,d,o

  real(kind=dp),dimension(K)  :: dx,dy,k0,k1
  real(kind=dp)               :: y0,y1,x1,x0
  integer                     :: i
  
  ! space dimension 
  !K  = size(x) 
  do i=1,k-1
     k0(i) = i
     k1(i) = i+1
  enddo
  !
  ! calculation of  dx
  !
  x1 = 2.0 * x(k) - x(k-1)
  x0 = 2.0 * x(1) - x(2)

  do i=3,k
     dx(i-1) = (x(i)-x(i-2)) / 2.0
  enddo
  dx(k)   = (x1 - x(k-1))/2.0
  dx(1)   = (x(2) - x0)/2.0
  !write(*,*) '-------------dans pdederive--------------'
  !print*,"dx(101)= ", dx(101),' x1=',x1,' x(K)=',x(K),K
  !print*,"y(101)= ",y(101)
  !print*,"o= ",o
  !write(*,*) '---------------------------'
  ! 
  ! dy calculation
  !
  if (o.eq.1) then

     if (c1 .ne. 2.0) then
        y1 = 2.0 * y(k) - y(k-1)
        y0 = 2.0 * y(1) - y(2)
     else
        y1 = 4.0 * y(k) - 6.0 * y(k-1) + 4.0 * y(k-2) - y(k-3)
        y0 = 4.0 * y(1) - 6.0 * y(2) + 4.0 * y(3) - y(4)
     endif
     do i=3,k
        dy(i-1) = y(i)-y(i-2)
     enddo

     dy(k)   = y1   - y(k-1)
     dy(1)   = y(2) - y0

  else

     if (c1 .ne. 2.0) then
        y1 = y(k-1) 
        y0 = y(2)
     else
        y1 = 4.0 * y(k) - 6.0 * y(k-1) + 4.0 * y(k-2) - y(k-3)
        y0 = 4.0 * y(1) - 6.0 * y(2) + 4.0 * y(3) - y(4)
     endif

     do i=2,k-1
        dy(i) = y(i+1)-2.0*y(i)+y(i-1)
     enddo
     dy(k)   = y1   - 2.0*y(k) + y(k-1)
     dy(1)   = y(2) - 2.0*y(1) + y0

  endif
  
  do i=1,k
     if (o.eq. 1) then
        dydx(i)  = dy(i) / (2.0 * dx(i)) 
     else
        dydx(i)  = dy(i) / (dx(i) ** 2.0)
     endif
  enddo

  if (c1 .eq. 0) then
     dydx(k) = c1 * dydx(k)
  endif
  if (c0 .eq. 0) then
     dydx(1) = c0 * dydx(1)
  endif

  !print*,"dx   = ",dx(100:101)
  !print*,"dy   = ",dy(100:101)
  !print*,"dydx = ",dydx(100:101)
end subroutine cos_rpdederive