source: LMDZ5/branches/IPSLCM6.0.10/tools/Max_diff_nc_with_lib/NR_util/poly.f90

Last change on this file was 1910, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1860:1909 into testing branch
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 4.8 KB
Line 
1MODULE poly_m
2
3  IMPLICIT NONE
4
5  INTEGER, PARAMETER, private :: NPAR_POLY=8
6
7  INTERFACE poly
8     MODULE PROCEDURE poly_rr,poly_rrv,poly_dd,poly_ddv,&
9          poly_rc,poly_cc,poly_msk_rrv,poly_msk_ddv
10  END INTERFACE
11
12  private poly_rr,poly_rrv,poly_dd,poly_ddv,&
13       poly_rc,poly_cc,poly_msk_rrv,poly_msk_ddv
14
15CONTAINS
16
17  FUNCTION poly_rr(x,coeffs)
18    REAL, INTENT(IN) :: x
19    REAL, DIMENSION(:), INTENT(IN) :: coeffs
20    REAL :: poly_rr
21    REAL :: pow
22    REAL, DIMENSION(:), ALLOCATABLE :: vec
23    INTEGER :: i,n,nn
24    n=size(coeffs)
25    if (n <= 0) then
26       poly_rr=0.0
27    else if (n < NPAR_POLY) then
28       poly_rr=coeffs(n)
29       do i=n-1,1,-1
30          poly_rr=x*poly_rr+coeffs(i)
31       end do
32    else
33       allocate(vec(n+1))
34       pow=x
35       vec(1:n)=coeffs
36       do
37          vec(n+1)=0.0
38          nn=ishft(n+1,-1)
39          vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2)
40          if (nn == 1) exit
41          pow=pow*pow
42          n=nn
43       end do
44       poly_rr=vec(1)
45       deallocate(vec)
46    end if
47  END FUNCTION poly_rr
48  !BL
49  FUNCTION poly_dd(x,coeffs)
50    DOUBLE PRECISION, INTENT(IN) :: x
51    DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: coeffs
52    DOUBLE PRECISION :: poly_dd
53    DOUBLE PRECISION :: pow
54    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: vec
55    INTEGER :: i,n,nn
56    n=size(coeffs)
57    if (n <= 0) then
58       poly_dd=0.0d0
59    else if (n < NPAR_POLY) then
60       poly_dd=coeffs(n)
61       do i=n-1,1,-1
62          poly_dd=x*poly_dd+coeffs(i)
63       end do
64    else
65       allocate(vec(n+1))
66       pow=x
67       vec(1:n)=coeffs
68       do
69          vec(n+1)=0.0d0
70          nn=ishft(n+1,-1)
71          vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2)
72          if (nn == 1) exit
73          pow=pow*pow
74          n=nn
75       end do
76       poly_dd=vec(1)
77       deallocate(vec)
78    end if
79  END FUNCTION poly_dd
80  !BL
81  FUNCTION poly_rc(x,coeffs)
82    COMPLEX, INTENT(IN) :: x
83    REAL, DIMENSION(:), INTENT(IN) :: coeffs
84    COMPLEX :: poly_rc
85    COMPLEX :: pow
86    COMPLEX, DIMENSION(:), ALLOCATABLE :: vec
87    INTEGER :: i,n,nn
88    n=size(coeffs)
89    if (n <= 0) then
90       poly_rc=0.0
91    else if (n < NPAR_POLY) then
92       poly_rc=coeffs(n)
93       do i=n-1,1,-1
94          poly_rc=x*poly_rc+coeffs(i)
95       end do
96    else
97       allocate(vec(n+1))
98       pow=x
99       vec(1:n)=coeffs
100       do
101          vec(n+1)=0.0
102          nn=ishft(n+1,-1)
103          vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2)
104          if (nn == 1) exit
105          pow=pow*pow
106          n=nn
107       end do
108       poly_rc=vec(1)
109       deallocate(vec)
110    end if
111  END FUNCTION poly_rc
112  !BL
113  FUNCTION poly_cc(x,coeffs)
114    COMPLEX, INTENT(IN) :: x
115    COMPLEX, DIMENSION(:), INTENT(IN) :: coeffs
116    COMPLEX :: poly_cc
117    COMPLEX :: pow
118    COMPLEX, DIMENSION(:), ALLOCATABLE :: vec
119    INTEGER :: i,n,nn
120    n=size(coeffs)
121    if (n <= 0) then
122       poly_cc=0.0
123    else if (n < NPAR_POLY) then
124       poly_cc=coeffs(n)
125       do i=n-1,1,-1
126          poly_cc=x*poly_cc+coeffs(i)
127       end do
128    else
129       allocate(vec(n+1))
130       pow=x
131       vec(1:n)=coeffs
132       do
133          vec(n+1)=0.0
134          nn=ishft(n+1,-1)
135          vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2)
136          if (nn == 1) exit
137          pow=pow*pow
138          n=nn
139       end do
140       poly_cc=vec(1)
141       deallocate(vec)
142    end if
143  END FUNCTION poly_cc
144  !BL
145  FUNCTION poly_rrv(x,coeffs)
146    REAL, DIMENSION(:), INTENT(IN) :: coeffs,x
147    REAL, DIMENSION(size(x)) :: poly_rrv
148    INTEGER :: i,n,m
149    m=size(coeffs)
150    n=size(x)
151    if (m <= 0) then
152       poly_rrv=0.0
153    else if (m < n .or. m < NPAR_POLY) then
154       poly_rrv=coeffs(m)
155       do i=m-1,1,-1
156          poly_rrv=x*poly_rrv+coeffs(i)
157       end do
158    else
159       do i=1,n
160          poly_rrv(i)=poly_rr(x(i),coeffs)
161       end do
162    end if
163  END FUNCTION poly_rrv
164  !BL
165  FUNCTION poly_ddv(x,coeffs)
166    DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: coeffs,x
167    DOUBLE PRECISION, DIMENSION(size(x)) :: poly_ddv
168    INTEGER :: i,n,m
169    m=size(coeffs)
170    n=size(x)
171    if (m <= 0) then
172       poly_ddv=0.0d0
173    else if (m < n .or. m < NPAR_POLY) then
174       poly_ddv=coeffs(m)
175       do i=m-1,1,-1
176          poly_ddv=x*poly_ddv+coeffs(i)
177       end do
178    else
179       do i=1,n
180          poly_ddv(i)=poly_dd(x(i),coeffs)
181       end do
182    end if
183  END FUNCTION poly_ddv
184  !BL
185  FUNCTION poly_msk_rrv(x,coeffs,mask)
186    REAL, DIMENSION(:), INTENT(IN) :: coeffs,x
187    LOGICAL, DIMENSION(:), INTENT(IN) :: mask
188    REAL, DIMENSION(size(x)) :: poly_msk_rrv
189    poly_msk_rrv=unpack(poly_rrv(pack(x,mask),coeffs),mask,0.0)
190  END FUNCTION poly_msk_rrv
191  !BL
192  FUNCTION poly_msk_ddv(x,coeffs,mask)
193    DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: coeffs,x
194    LOGICAL, DIMENSION(:), INTENT(IN) :: mask
195    DOUBLE PRECISION, DIMENSION(size(x)) :: poly_msk_ddv
196    poly_msk_ddv=unpack(poly_ddv(pack(x,mask),coeffs),mask,0.0d0)
197  END FUNCTION poly_msk_ddv
198
199END MODULE poly_m
Note: See TracBrowser for help on using the repository browser.