V3FIT
splinint.f
1  SUBROUTINE splinint(grn, cm, jacob, h, u, u1, w, w1, nk, nots,
2  1 ifunc, nmesh)
3  USE vparams
4  USE vspline
5  IMPLICIT NONE
6 C-----------------------------------------------
7 C D u m m y A r g u m e n t s
8 C-----------------------------------------------
9  INTEGER nots, ifunc, nmesh
10  INTEGER, DIMENSION(*) :: nk
11  REAL(rprec), DIMENSION(nmesh) :: grn, cm
12  REAL(rprec), DIMENSION(nots) :: jacob, h
13  REAL(rprec), DIMENSION(nmesh) :: u, u1, w, w1
14 C-----------------------------------------------
15 C L o c a l V a r i a b l e s
16 C-----------------------------------------------
17  INTEGER :: j, k, nb, ia, ib, k1, ksp1, nmesh1
18  REAL(rprec), DIMENSION(nmesh) :: func
19  REAL(rprec), DIMENSION(nots) :: af, bs
20 C-----------------------------------------------
21 
22 !
23 ! COMPUTES af,bs FACTORS IN
24 ! (ifunc=INTFUN) Int[ GRN(X) * F(X) ds ] or
25 ! (ifunc=INTDER) Int[ GRN(X) * {d(CmF)/ds} ds ]
26 ! = af(k)*f(k) + bs(k)*f''(k)
27 ! WHERE f(k),f''(k) ARE SPLINE COEFFICIENTS OF F(X)
28 !
29 ! NOTE: FOR ifunc = INTDER, the OHS factor in CmF cancels
30 ! THE HS factor in the Integral.
31 ! FOR ifunc = INTFUN, GRN is assumed to be pre-multiplied
32 ! OUTSIDE this routine by HS factor
33 ! ALSO, COMPUTES af(k) + (SUM on i)bs(i)*J(i,k) = jacob(k),
34 ! WHERE J(i,k) = d[g(i)]/d[f(k)], g = f''
35 !
36 ! nk(k): Number of smesh-pts in k-th spline interval
37 ! xknots(k) < smesh <= xknots(k+1), k = 1,nots-1
38 !
39 ! NOTE: The ifunc=INTDER CASE is done by integrating by parts,
40 ! so that the half-point integration (GRN at half mesh pts)
41 ! becomes a full-point integration in Cm*F.
42 !
43  nb = ideriv !Pressure, iota derivatives vanish at origin
44  ksp1 = nots - 1
45  nmesh1 = nmesh - 1
46 
47  IF (ifunc .eq. intder) THEN
48 !
49 ! Integrate by parts (in s), func(1) and func(nmesh) are 'surface terms'
50 !
51  func(1) = -cm(1)*grn(2)
52  func(2:nmesh1) = cm(2:nmesh1)*(grn(2:nmesh1)-grn(3:nmesh1+1))
53  func(nmesh) = cm(nmesh)*grn(nmesh)
54  ELSE
55  func = grn
56  ENDIF
57 
58  af(:nots) = zero
59  bs(:nots) = zero
60 
61  ia = 1
62  DO k = 1, ksp1
63  IF (nk(k) .ne. 0) THEN
64  k1 = k + 1
65  ib = ia + nk(k) - 1
66  DO j = ia,ib
67  af(k) = af(k) + func(j)*w(j)
68  bs(k) = bs(k) + func(j)*u(j)
69  af(k1) = af(k1) + func(j)*w1(j)
70  bs(k1) = bs(k1) + func(j)*u1(j)
71  ENDDO
72  ia = ib + 1
73  ENDIF
74  END DO
75 
76  IF (ib .ne. nmesh) stop 'ib!=nmesh'
77  IF (nb .eq. natur) bs(1) = 0. !Natural boundary conditions
78  bs(nots) = 0.
79  CALL jacprod (bs, h, nots, nb) !Returns bs(i)=bs(j)*J(j,i)
80  jacob(:nots) = af(:nots) + bs(:nots)
81 
82  END SUBROUTINE splinint