-
Notifications
You must be signed in to change notification settings - Fork 11
Expand file tree
/
Copy pathlin_interp_ref.F90
More file actions
62 lines (45 loc) · 1.78 KB
/
Copy pathlin_interp_ref.F90
File metadata and controls
62 lines (45 loc) · 1.78 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
module ekat_ut
use iso_c_binding, only: c_int, c_double, c_float
implicit none
#ifdef EKAT_TEST_DOUBLE_PRECISION
integer, parameter :: c_real = c_double
#else
integer, parameter :: c_real = c_float
#endif
contains
subroutine linear_interp_c(x1,x2,y1,y2,km1,km2,ncol) bind(c)
integer(kind=c_int), value, intent(in) :: km1, km2, ncol
real(kind=c_real), intent(in) :: x1(ncol,km1), y1(ncol,km1)
real(kind=c_real), intent(in) :: x2(ncol,km2)
real(kind=c_real), intent(out) :: y2(ncol,km2)
call linear_interp(x1,x2,y1,y2,km1,km2,ncol)
end subroutine linear_interp_c
!==============================================================
! Linear interpolation to get values on various grids
subroutine linear_interp(x1,x2,y1,y2,km1,km2,ncol)
use iso_c_binding
implicit none
integer(kind=c_int), intent(in) :: km1, km2
integer(kind=c_int), intent(in) :: ncol
real(kind=c_real), intent(in) :: x1(ncol,km1), y1(ncol,km1)
real(kind=c_real), intent(in) :: x2(ncol,km2)
real(kind=c_real), intent(out) :: y2(ncol,km2)
integer :: k1, k2, i
do i=1,ncol
do k2=1,km2
if( x2(i,k2) <= x1(i,1) ) then
y2(i,k2) = y1(i,1) + (y1(i,2)-y1(i,1))*(x2(i,k2)-x1(i,1))/(x1(i,2)-x1(i,1))
elseif( x2(i,k2) >= x1(i,km1) ) then
y2(i,k2) = y1(i,km1) + (y1(i,km1)-y1(i,km1-1))*(x2(i,k2)-x1(i,km1))/(x1(i,km1)-x1(i,km1-1))
else
do k1 = 2,km1
if( (x2(i,k2)>=x1(i,k1-1)).and.(x2(i,k2)<x1(i,k1)) ) then
y2(i,k2) = y1(i,k1-1) + (y1(i,k1)-y1(i,k1-1))*(x2(i,k2)-x1(i,k1-1))/(x1(i,k1)-x1(i,k1-1))
endif
enddo ! end k1 loop
endif
enddo ! end k2 loop
enddo ! i loop
return
end subroutine linear_interp
end module ekat_ut