The TELEMAC-MASCARET system  trunk
stirling.f
Go to the documentation of this file.
1 ! *******************
2  SUBROUTINE stirling
3 ! *******************
4 !
5  &(ni,xi,yi,no,xostep,yo)
6 !
7 !***********************************************************************
8 ! ARTEMIS V7P0 June 2014
9 !***********************************************************************
10 !
11 !brief INTERPOLATES A FUNCTION TO DESIRED RESOLUTION.
12 !
13 !history N.DURAND (HRW)
14 !+ March 2001
15 !+
16 !+ Ported to V7P0
17 !
18 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
19 !| NI XI YI | | INPUT NUMBER OF DISCRETE VALUES,
20 !| | | INPUT X AND Y VALUES
21 !| NO XOSTEP YO | | OUTPUT NUMBER OF DISCRETE VALUES,
22 !| | | CORRESPONDING OUTPUT RESOLUTION IN X,
23 !| | | AND OUTPUT Y VALUES
24 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
25 !
27  IMPLICIT NONE
28 !
29 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
30 !
31  INTEGER :: NI,NO
32  DOUBLE PRECISION :: XI(ni),YI(ni)
33  DOUBLE PRECISION :: XOSTEP,XO,YO(no)
34 !
35  INTEGER :: I,IV,IVP
36  DOUBLE PRECISION :: A1,A2,A3,A4,A5,A6
37  DOUBLE PRECISION :: RMIN,R,R2,U
38 !
39  INTRINSIC float, dabs
40 !
41 !-----------------------------------------------------------------------
42  xostep = (xi(ni) - xi(1)) / float(no-1)
43  yo(1) = yi(1)
44  yo(no) = yi(ni)
45 !
46 !-----------------------------------------------------------------------
47 ! PARAMETERS FOR THE PARABOLA FITTING
48 ! OUTSIDE DO AND IF LOOPS FOR OPTIMISATION PURPOSES
49 ! ESTIMATED ONLY ONCE
50  a3 = ( yi(3)-yi(1) )/( xi(3)-xi(1) )/( xi(3)-xi(2) ) -
51  & ( yi(2)-yi(1) )/( xi(2)-xi(1) )/( xi(3)-xi(2) )
52  a2 = ( yi(2)-yi(1) )/( xi(2)-xi(1) ) - a3*( xi(2)+xi(1) )
53  a1 = yi(1) - ( a2 + a3*xi(1) )*xi(1)
54 !
55  a4 = ( yi(ni-2)-yi(ni) )/( xi(ni-2)-xi(ni) )/( xi(ni-2)-xi(ni-1) )
56  & - ( yi(ni-1)-yi(ni) )/( xi(ni-1)-xi(ni) )/( xi(ni-2)-xi(ni-1) )
57  a5 = ( yi(ni-1)-yi(ni) )/( xi(ni-1)-xi(ni) )
58  & - a4*( xi(ni-1)+xi(ni) )
59  a6 = yi(ni) - ( a5 + a4*xi(ni) )*xi(ni)
60 !
61 !-----------------------------------------------------------------------
62 ! LOOP ON THE INTERPOLATION POINTS IN BETWEEN
63  ivp = 2
64  DO i = 2,(no-1)
65  rmin = 1.d0
66  xo = xi(1)+float(i-1)*xostep
67 !
68 !-----------------------------------------------------------------------
69 ! DETERMINES IVP SUCH THAT XO IS CLOSEST TO XI(IVP)
70 ! R<0 FOR XO BEFORE XI(IVP); R>0 FOR XO AFTER XI(IVP)
71  DO iv = ivp,(ni-1)
72  r = ( xo - xi(iv) ) / ( xi(iv+1) - xi(iv) )
73  IF( dabs(r).LE.dabs(rmin) ) THEN
74  rmin = r
75  ivp = iv
76  ENDIF
77  ENDDO ! IV = IVP,(NI-1)
78  r = ( xo - xi(ni) ) / ( xi(ni) - xi(ni-1) )
79  IF( dabs(r).LE.dabs(rmin) ) THEN
80  rmin = r
81  ivp = ni
82  ENDIF
83  r = rmin
84 !
85 !-----------------------------------------------------------------------
86 ! XO CLOSEST TO XI(2) :
87 ! YO COMPUTED FROM PARABOLA USING 1ST 3 POINTS
88  IF( ivp.EQ.2 ) THEN
89  yo(i) = a1 + a2*xo + a3*xo**2
90 !
91 !-----------------------------------------------------------------------
92 ! XO CLOSEST TO XI(3) OR XI(NI-2) :
93 ! 2ND ORDER STIRLING WITH CONSTANT STEP
94  ELSEIF( ivp.EQ.3 .OR. ivp.EQ.ni-2 ) THEN
95  r2 = r**2
96  yo(i) = ( r2-r )*yi(ivp-1)/2.d0 +
97  & ( 1.d0-r2 )*yi(ivp) + ( r2+r )*yi(ivp+1)/2.d0
98 !
99 !-----------------------------------------------------------------------
100 ! XO CLOSEST TO XI(NI-1) OR XI(NI) :
101 ! YO COMPUTED FROM PARABOLA USING LAST 3 POINTS
102  ELSEIF( ivp.EQ.ni-1 .OR. ivp.EQ.ni ) THEN
103  yo(i) = a6 + a5*xo + a4*xo**2
104 !
105 !-----------------------------------------------------------------------
106 ! SUBSEQUENT POINTS
107 ! 4TH ORDER STIRLING WITH CONSTANT STEP
108  ELSE
109  u = r*( (r**2)-1.d0 )/12.d0
110  yo(i) = u*( (r/2.d0)-1.d0 ) * yi(ivp-2) +
111  & (1.d0-r)*(2.d0*u-(r/2.d0)) * yi(ivp-1) +
112  & ( r*(3.d0*u-r) + 1.d0 ) * yi(ivp) +
113  & (1.d0+r)*(-2.d0*u+(r/2.d0)) * yi(ivp+1) +
114  & u*( (r/2.d0)+1.d0 ) * yi(ivp+2)
115  ENDIF
116 !
117 !-----------------------------------------------------------------------
118 ! END OF THE LOOP ON THE INTERPOLATION POINTS
119  ENDDO ! I = 2,(NO-1)
120 !
121  RETURN
122  END
subroutine stirling(NI, XI, YI, NO, XOSTEP, YO)
Definition: stirling.f:7