source: trunk/LMDZ.TITAN/libf/phytitan/interp_line.F @ 3000

Last change on this file since 3000 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

File size: 1.8 KB
Line 
1      subroutine interp_line(x1,y1,len1,x2,y2,len2)
2      implicit none
3!-----------------------------------------------------------------------
4!
5!  Purpose: Do a series of linear interpolations
6!  Data sets are organized as vectors (see below)
7!  If x2(:), and abscissa at which interpolation is requiered, lies
8!  outside of the interval covered by x1(:), instead of doing an
9!  extrapolation, y2() is set to the value y1() corresponding to
10!  the nearby x1(:) point
11
12c-----------------------------------------------------------------------
13!  arguments
14!  ---------
15!  inputs:
16      real x1(len1) ! ordered list of abscissas
17      real y1(len1) ! values at x1(:)
18      integer len1  ! length of x1(:) and y1(:)
19      real x2(len2) !ordered list of abscissas at which interpolation is done
20      integer len2  ! length of x2(:) and y2(:)
21!  outputs:
22      real y2(len2) ! interpolated values
23!-----------------------------------------------------------------------
24
25! local variables:
26      integer i,j
27     
28
29      do i=1,len2
30        ! check if x2(i) lies outside of the interval covered by x1()
31        if(((x2(i).le.x1(1)).and.(x2(i).le.x1(len1))).or.
32     &     ((x2(i).ge.x1(1)).and.(x2(i).ge.x1(len1)))) then
33          ! set y2(i) to y1(1) or y1(len1)
34          if (abs(x2(i)-x1(1)).lt.abs(x2(i)-x1(len1))) then
35            ! x2(i) lies closest to x1(1)
36            y2(i)=y1(1)
37          else
38            ! x2(i) lies closest to x1(len1)
39            y2(i)=y1(len1)
40          endif
41
42        else
43        ! find the nearest neigbours and do a linear interpolation
44         do j=1,len1-1
45          if(((x2(i).ge.x1(j)).and.(x2(i).le.x1(j+1))).or.
46     &       ((x2(i).le.x1(j)).and.(x2(i).ge.x1(j+1)))) then
47            y2(i)=((x2(i)-x1(j))/(x1(j+1)-x1(j)))*y1(j+1)+
48     &            ((x2(i)-x1(j+1))/(x1(j)-x1(j+1)))*y1(j)
49          endif
50         enddo
51        endif
52
53      enddo
54
55      end
Note: See TracBrowser for help on using the repository browser.