Line data    Source code

       1             :       MODULE m_dujdr
       2             :       CONTAINS
       3             : 
       4           0 :       SUBROUTINE dujdr(jmtd,jri,rmsh,dx,ub,j2,b,l,lmaxd,dub)
       5             : 
       6             :       USE m_constants
       7             :       USE m_difcub
       8             : 
       9             :       IMPLICIT NONE
      10             :       
      11             :       INTEGER, INTENT(IN) :: jri,jmtd,lmaxd
      12             :       REAL, INTENT(IN) :: rmsh(jmtd),dx
      13             :       REAL, INTENT(IN) :: ub(jmtd,2)       ! u(b2)
      14             :       REAL, INTENT(IN) :: j2(0:lmaxd,jmtd) ! j_l(b2*r)
      15             :       REAL, INTENT(IN) :: b                ! b2
      16             :       INTEGER, INTENT(IN) :: l           ! l of sph. Bessel j2
      17             : 
      18             :       REAL, INTENT(OUT) :: dub(jmtd,2)
      19           0 :       REAL :: xi,t(jri,2)
      20             :       INTEGER :: i,j
      21             : 
      22             :       ! derivatives d/dr for large and small component of q
      23           0 :       DO i=1,jri
      24           0 :        t(i,:) = ub(i,:) / rmsh(i) * j2(l,i)
      25             :       ENDDO
      26             : 
      27           0 :       DO j = 1, 2
      28             :          ! derivative at 1st point
      29           0 :          dub(1,j) = difcub( rmsh(1),t(1,j),rmsh(1) )
      30             : 
      31             :          ! derivative at 2nd...(jri-2)th point
      32           0 :          DO i = 2, jri-2
      33           0 :             dub(i,j) = difcub( rmsh(i-1),t(i-1,j),rmsh(i) )
      34             :          ENDDO
      35             : 
      36             :          ! derivative at last two points
      37           0 :          dub(jri-1,j) = difcub( rmsh(jri-3),t(jri-3,j),rmsh(jri-1) )
      38           0 :          dub(jri,j) = difcub( rmsh(jri-3),t(jri-3,j),rmsh(jri) )
      39             :       ENDDO
      40             : 
      41           0 :       DO i=1,jri
      42           0 :        dub(i,:) = dub(i,:)*rmsh(i)
      43             :       ENDDO
      44             : 
      45             :       ! complete d/dr (ub*j2) = ub'j2 + ub j2' with sph. Bessel func. j
      46             :       ! rule:  j'_{l}(ar) = a*j_{l-1}(ar) - (l+1)/r*j_{l}(ar)
      47             : c      IF(l.ne.0) THEN
      48             : c       DO i=1,jri
      49             : c        xi = rmsh(i)
      50             : c        dub(i,:) = dub(i,:) * j2(l,i) * xi
      51             : c     >           + ub(i,:) *( j2(l-1,i)*b - (l+1)/xi*j2(l,i) ) 
      52             : c       ENDDO
      53             : c      ELSE
      54             : c       DO i=1,jri
      55             : c        xi = rmsh(i)
      56             : c        dub(i,:) = dub(i,:) * j2(l,i) * xi
      57             : c     >           - ub(i,:) * j2(1,i) * b
      58             : c       ENDDO
      59             : c      ENDIF
      60             :       
      61           0 :       END SUBROUTINE dujdr
      62             :       END MODULE m_dujdr