Line data Source code
1 : !----------------------------------------------------------------------------
2 : ! Utility module used for interpolation of aerosol optics table
3 : ! NOTE: Results will be set to table edges for interpolations beyond
4 : ! the edges -- no extropolations
5 : !----------------------------------------------------------------------------
6 : module table_interp_mod
7 : use shr_kind_mod, only: r8=>shr_kind_r8
8 :
9 : implicit none
10 :
11 : private
12 : public :: table_interp
13 : public :: table_interp_wghts
14 : public :: table_interp_calcwghts
15 :
16 : ! overload the interpolation routines
17 : interface table_interp
18 : module procedure interp1d
19 : module procedure interp2d
20 : module procedure interp4d
21 : end interface table_interp
22 :
23 : ! interpolation weights and indices
24 : type :: table_interp_wghts
25 : real(r8) :: wt1
26 : real(r8) :: wt2
27 : integer :: ix1
28 : integer :: ix2
29 : end type table_interp_wghts
30 :
31 : contains
32 :
33 : !--------------------------------------------------------------------------
34 : ! 1-D interpolation
35 : !--------------------------------------------------------------------------
36 0 : pure function interp1d( ncol, nxs, xwghts, tbl ) result(res)
37 :
38 : integer, intent(in) :: ncol ! number of model columns
39 : integer, intent(in) :: nxs ! table size
40 : real(r8), intent(in) :: tbl(nxs) ! table values to be interpolated
41 : type(table_interp_wghts), intent(in) :: xwghts(ncol) ! interpolation weights and indices
42 :
43 : real(r8) :: res(ncol)
44 :
45 : integer :: i
46 :
47 0 : do i = 1,ncol
48 :
49 0 : res(i) = xwghts(i)%wt1*tbl(xwghts(i)%ix1) &
50 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2)
51 :
52 : end do
53 :
54 0 : end function interp1d
55 :
56 : !--------------------------------------------------------------------------
57 : ! 2-D interpolation
58 : !--------------------------------------------------------------------------
59 0 : pure function interp2d( ncoef, ncol, nxs, nys, xwghts, ywghts, tbl ) result(res)
60 :
61 : integer, intent(in) :: ncoef ! number chebyshev coefficients
62 : integer, intent(in) :: ncol ! number of model columns
63 : integer, intent(in) :: nxs ! table x-dimension size
64 : integer, intent(in) :: nys ! table y-dimension size
65 : real(r8), intent(in) :: tbl(ncoef,nxs,nys) ! table values to be interpolated
66 : type(table_interp_wghts), intent(in) :: xwghts(ncol) ! x interpolation weights and indices
67 : type(table_interp_wghts), intent(in) :: ywghts(ncol) ! y interpolation weights and indices
68 :
69 : real(r8) :: res(ncoef,ncol)
70 :
71 0 : real(r8) :: fx(ncoef,2)
72 :
73 : integer :: i
74 :
75 0 : do i = 1,ncol
76 :
77 : ! interp x dir
78 0 : fx(:,1) = xwghts(i)%wt1*tbl(:,xwghts(i)%ix1,ywghts(i)%ix1) & ! @ y1
79 0 : + xwghts(i)%wt2*tbl(:,xwghts(i)%ix2,ywghts(i)%ix1)
80 0 : fx(:,2) = xwghts(i)%wt1*tbl(:,xwghts(i)%ix1,ywghts(i)%ix2) & ! @ y2
81 0 : + xwghts(i)%wt2*tbl(:,xwghts(i)%ix2,ywghts(i)%ix2)
82 :
83 : ! interp y dir
84 0 : res(:,i) = ywghts(i)%wt1*fx(:,1) + ywghts(i)%wt2*fx(:,2)
85 :
86 : end do
87 :
88 0 : end function interp2d
89 :
90 : !--------------------------------------------------------------------------
91 : ! 4-D interpolation
92 : !--------------------------------------------------------------------------
93 0 : pure function interp4d( ncol, nxs, nys, nzs, nts, xwghts, ywghts, zwghts, twghts, tbl ) result(res)
94 :
95 : integer, intent(in) :: ncol ! number of model columns
96 : integer, intent(in) :: nxs ! table x-dimension size
97 : integer, intent(in) :: nys ! table y-dimension size
98 : integer, intent(in) :: nzs ! table z-dimension size
99 : integer, intent(in) :: nts ! table t-dimension size
100 : real(r8), intent(in) :: tbl(nxs,nys,nzs,nts) ! table values to be interpolated
101 : type(table_interp_wghts), intent(in) :: xwghts(ncol) ! x interpolation weights and indices
102 : type(table_interp_wghts), intent(in) :: ywghts(ncol) ! y interpolation weights and indices
103 : type(table_interp_wghts), intent(in) :: zwghts(ncol) ! z interpolation weights and indices
104 : type(table_interp_wghts), intent(in) :: twghts(ncol) ! t interpolation weights and indices
105 :
106 : real(r8) :: res(ncol)
107 :
108 : real(r8) :: fx(8)
109 : real(r8) :: fy(4)
110 : real(r8) :: fz(2)
111 :
112 : integer :: i
113 :
114 0 : do i = 1,ncol
115 :
116 : ! interp x dir
117 0 : fx(1) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix1,zwghts(i)%ix1,twghts(i)%ix1) & ! @ y1, z1, t1
118 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix1,zwghts(i)%ix1,twghts(i)%ix1)
119 0 : fx(2) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix2,zwghts(i)%ix1,twghts(i)%ix1) & ! @ y2, z1, t1
120 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix2,zwghts(i)%ix1,twghts(i)%ix1)
121 :
122 0 : fx(3) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix1,zwghts(i)%ix2,twghts(i)%ix1) & ! @ y1, z2, t1
123 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix1,zwghts(i)%ix2,twghts(i)%ix1)
124 : fx(4) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix2,zwghts(i)%ix2,twghts(i)%ix1) & ! @ y2, z2, t1
125 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix2,zwghts(i)%ix2,twghts(i)%ix1)
126 :
127 0 : fx(5) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix1,zwghts(i)%ix1,twghts(i)%ix2) & ! @ y1, z1, t2
128 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix1,zwghts(i)%ix1,twghts(i)%ix2)
129 : fx(6) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix2,zwghts(i)%ix1,twghts(i)%ix2) & ! @ y2, z1, t2
130 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix2,zwghts(i)%ix1,twghts(i)%ix2)
131 :
132 : fx(7) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix1,zwghts(i)%ix2,twghts(i)%ix2) & ! @ y1, z2, t2
133 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix1,zwghts(i)%ix2,twghts(i)%ix2)
134 : fx(8) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix2,zwghts(i)%ix2,twghts(i)%ix2) & ! @ y2, z2, t2
135 0 : + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix2,zwghts(i)%ix2,twghts(i)%ix2)
136 :
137 : ! interp y dir
138 0 : fy(1) = ywghts(i)%wt1*fx(1) + ywghts(i)%wt2*fx(2) ! @ z1, t1
139 0 : fy(2) = ywghts(i)%wt1*fx(3) + ywghts(i)%wt2*fx(4) ! @ z2, t1
140 0 : fy(3) = ywghts(i)%wt1*fx(5) + ywghts(i)%wt2*fx(6) ! @ z1, t2
141 0 : fy(4) = ywghts(i)%wt1*fx(7) + ywghts(i)%wt2*fx(8) ! @ z2, t2
142 :
143 : ! interp z dir
144 0 : fz(1) = zwghts(i)%wt1*fy(1) + zwghts(i)%wt2*fy(2) ! @ t1
145 0 : fz(2) = zwghts(i)%wt1*fy(3) + zwghts(i)%wt2*fy(4) ! @ t2
146 :
147 : ! interp t dir
148 0 : res(i) = twghts(i)%wt1*fz(1) + twghts(i)%wt2*fz(2)
149 :
150 : end do
151 :
152 0 : end function interp4d
153 :
154 : !--------------------------------------------------------------------------
155 : ! determines interpolation weights and indices for given values at the model columns
156 : !--------------------------------------------------------------------------
157 0 : pure function table_interp_calcwghts( ngrid, xgrid, ncols, xcols ) result(wghts)
158 :
159 : integer, intent(in) :: ngrid ! number of grid point values
160 : real(r8), intent(in) :: xgrid(ngrid) ! grid point values
161 : integer, intent(in) :: ncols ! number of model columns
162 : real(r8), intent(in) :: xcols(ncols) ! values at the model columns
163 :
164 : type(table_interp_wghts) :: wghts(ncols) ! interpolations weights at the model columns
165 :
166 : integer :: i
167 0 : real(r8) :: xs(ncols)
168 :
169 0 : xs(:) = xcols(:)
170 :
171 : ! do not extrapolate beyond the edges of the table
172 0 : where(xs < xgrid(1))
173 0 : xs = xgrid(1)
174 : end where
175 0 : where(xs > xgrid(ngrid))
176 : xs = xgrid(ngrid)
177 : end where
178 :
179 0 : do i = 1,ncols
180 0 : wghts(i)%ix2 = find_index(ngrid,xgrid,xs(i))
181 0 : wghts(i)%ix1 = wghts(i)%ix2 - 1
182 0 : wghts(i)%wt1 = (xgrid(wghts(i)%ix2)-xs(i)) &
183 0 : /(xgrid(wghts(i)%ix2)-xgrid(wghts(i)%ix1))
184 0 : wghts(i)%wt2 = 1._r8 - wghts(i)%wt1
185 : end do
186 :
187 0 : end function table_interp_calcwghts
188 :
189 : ! private methods
190 : !--------------------------------------------------------------------------
191 : !--------------------------------------------------------------------------
192 : ! determines last index of grid vals of which is greater then or equal to
193 : ! value vx
194 : !--------------------------------------------------------------------------
195 0 : pure function find_index( nvals, vals, vx ) result(res)
196 : integer, intent(in) :: nvals
197 : real(r8), intent(in) :: vals(nvals)
198 : real(r8), intent(in) :: vx
199 : integer :: res
200 :
201 : integer :: ndx
202 :
203 0 : res = -1
204 :
205 0 : find_ndx: do ndx = 2, nvals
206 0 : if (vals(ndx)>=vx) then
207 : res = ndx
208 : exit find_ndx
209 : end if
210 : end do find_ndx
211 :
212 0 : end function find_index
213 :
214 0 : end module table_interp_mod
|
