NADCON5-ng  0.0.1
NADCON5 Next Generation
bicubic.f
Go to the documentation of this file.
1 c> \ingroup core
2 c> Subroutine to perform a 2-D cubic ("bicubic") interpolation
3 c>
4 c> Performs interpolation at location "xla,xlo" off of grid "z", whose
5 c> header information is the standard ".b" header
6 c> information with additional inputs
7 c>
8 c> \param[in] z Input Grid
9 c> \param[in] glamn minimum latitude (real*8 decimal degrees)
10 c> \param[in] glomn minimum longitude (real*8 decimal degrees)
11 c> \param[in] dla latitude spacing (real*8 decimal degrees)
12 c> \param[in] dlo longitude spacing (real*8 decimal degrees)
13 c> \param[in] nla number of lat rows (integer*4)
14 c> \param[in] nlo number of lon cols (integer*4)
15 c> \param[in] maxla actual dimensioned size of "z" in rows
16 c> \param[in] maxlo actual dimensioned size of "z" in cols
17 c> \param[in] xla lat of pt for interpolation (real*8 dec. deg)
18 c> \param[in] xlo lon of pt for interpolation (real*8 dec. deg)
19 c> \param[out] val interpolated value (real*8)
20 c>
21 c> ### Method:
22 c> Fit a 4x4 window over the random point. Unless the point
23 c> is less than one grid spacing from the edge of the grid,
24 c> it will fall in the inner 2x2 cell, and the 4x4 cell will
25 c> be centered upon that.
26 c>
27 c> Thus, if our point of interest is the asterisk, the window
28 c> will look like this:
29 c>
30 c> . . . .
31 c> . . . .
32 c> . .* . .
33 c> . . . .
34 c>
35  subroutine bicubic(z,glamn,glomn,dla,dlo,
36  *nla,nlo,maxla,maxlo,xla,xlo,val)
37 
38 c - Data is assumed real*4
39  implicit real*8 (a-h,o-z)
40  real*4 z(maxla,maxlo)
41 
42 c - Internal use
43  real*4 x,y,fx0,fx1,fx2,fx3,cubterp
44 
45 c - Find which row should be LEAST when fitting
46 c - a 4x4 window around xla.
47  difla = (xla - glamn)
48  ratla = difla / dla
49  jla = int(ratla)
50 c - Fix any edge overlaps
51  if(jla.lt.1)jla = 1
52  if(jla.gt.nla-3)jla = nla-3
53 
54 c - Find which column should be LEAST when fitting
55 c - a 4x4 window around xlo.
56  diflo = (xlo - glomn)
57  ratlo = diflo / dlo
58  jlo = int(ratlo)
59 c - Fix any edge overlaps
60  if(jlo.lt.1)jlo = 1
61  if(jlo.gt.nlo-3)jlo = nlo-3
62 
63 c - In the range of 0(westernmost) to
64 c - 3(easternmost) col, where is our
65 c - random lon value? That is, x must
66 c - be real and fall between 0 and 3.
67 c - (Note, unless it is near an edge,
68 c - it will always be between 1 and 2)
69 
70  x=(xlo-dlo*(jlo-1)-glomn)/dlo
71 
72  if(x.lt.0.0)then
73  write(6,100)x
74  stop
75  endif
76  100 format(
77  *'FATAL in bicubic: x<0 : ',f20.10)
78 
79  if(x.gt.3.0)then
80  write(6,101)x
81  stop
82  endif
83  101 format(
84  *'FATAL in bicubic: x>3 : ',f20.10)
85 
86 c - In the range of 0(southernmost) to
87 c - 3(northernmost) row, where is our
88 c - random lat value? That is, x must
89 c - be real and fall between 0 and 3.
90 c - (Note, unless it is near an edge,
91 c - it will always be between 1 and 2)
92 
93  y=(xla-dla*(jla-1)-glamn)/dla
94 
95  if(y.lt.0.0)then
96  write(6,102)y
97  stop
98  endif
99  102 format(
100  *'FATAL in bicubic: y<0 : ',f20.10)
101 
102  if(y.gt.3.0)then
103  write(6,103)y
104  stop
105  endif
106  103 format(
107  *'FATAL in bicubic: y>3 : ',f20.10)
108 
109 c - Now do the interpolation. First, build a 1-D cubic function
110 c - east-west in the southermost row and interpolate to longitude
111 c - "xlo" (at "x" for 0<x<3). Then do it in the 2nd row,
112 c - then again in the 3rd, and finally the 4th and most northern row.
113 c - The last step is to fit a 1-D cubic function north-south at the
114 c - four previous interpolations, but this time to interpolate to
115 c - latitude "xla" (at "y" for 0<y<3). Obviously we
116 c - could reverse the order, doing 4 north-south cubics
117 c - followed by one east-east and we'd get the same answer.
118 
119 
120  fx0=cubterp(x,z(jla ,jlo ),z(jla ,jlo+1),
121  * z(jla ,jlo+2),z(jla ,jlo+3))
122 
123  fx1=cubterp(x,z(jla+1,jlo ),z(jla+1,jlo+1),
124  * z(jla+1,jlo+2),z(jla+1,jlo+3))
125 
126  fx2=cubterp(x,z(jla+2,jlo ),z(jla+2,jlo+1),
127  * z(jla+2,jlo+2),z(jla+2,jlo+3))
128 
129  fx3=cubterp(x,z(jla+3,jlo ),z(jla+3,jlo+1),
130  * z(jla+3,jlo+2),z(jla+3,jlo+3))
131 
132  val=dble(cubterp(y,fx0,fx1,fx2,fx3))
133 
134 
135 c write(6,1001)glamn,xla,jla,y
136 c write(6,1002)glomn,xlo,jlo,x
137 c write(6,1003)z(jla+3,jlo ),z(jla+3,jlo+1),
138 c * z(jla+3,jlo+2),z(jla+3,jlo+3)
139 c write(6,1003)z(jla+2,jlo ),z(jla+2,jlo+1),
140 c * z(jla+2,jlo+2),z(jla+2,jlo+3)
141 c write(6,1003)z(jla+1,jlo ),z(jla+1,jlo+1),
142 c * z(jla+1,jlo+2),z(jla+1,jlo+3)
143 c write(6,1003)z(jla ,jlo ),z(jla ,jlo+1),
144 c * z(jla ,jlo+2),z(jla ,jlo+3)
145 c write(6,1004)fx0,fx1,fx2,fx3,val
146 c pause
147 
148  1001 format(f14.10,1x,f14.10,1x,i8,1x,f8.5)
149  1002 format(f14.10,1x,f14.10,1x,i8,1x,f8.5)
150  1003 format(4(f15.8,1x))
151  1004 format(
152  *'fx0 = ',f15.8,/,
153  *'fx1 = ',f15.8,/,
154  *'fx2 = ',f15.8,/,
155  *'fx3 = ',f15.8,/,
156  *'val = ',f15.8)
157 
158 
159 
160  return
161  end
162  include 'cubterp.f'
subroutine bicubic(z, glamn, glomn, dla, dlo, nla, nlo, maxla, maxlo, xla, xlo, val)
Subroutine to perform a 2-D cubic ("bicubic") interpolation.
Definition: bicubic.f:37