NADCON5-ng  0.0.2
NADCON5 Next Generation Documentation
biquad.f
Go to the documentation of this file.
1 c> \ingroup core
2 c> \if MANPAGE
3 c> \page biquad
4 c> \endif
5 c>
6 c> Subroutine to perform a 2-D quadratic ("biquadratic") interpolation
7 c>
8 c> Performs a biquadratic interpolation at location
9 c> `xla,xlo` off of grid `z`, whose
10 c> header information is the standard ".b" header
11 c> information
12 c>
13 c> \param[in] z Input Grid
14 c> \param[in] glamn minimum latitude (real*8 decimal degrees)
15 c> \param[in] glomn minimum longitude (real*8 decimal degrees)
16 c> \param[in] dla latitude spacing (real*8 decimal degrees)
17 c> \param[in] dlo longitude spacing (real*8 decimal degrees)
18 c> \param[in] nla number of lat rows (integer*4)
19 c> \param[in] nlo number of lon cols (integer*4)
20 c> \param[in] maxla actual dimensioned size of "z" in rows
21 c> \param[in] maxlo actual dimensioned size of "z" in cols
22 c> \param[in] xla lat of pt for interpolation (real*8 dec. deg)
23 c> \param[in] xlo lon of pt for interpolation (real*8 dec. deg)
24 c> \param[out] val interpolated value (real*8)
25 c>
26 c> ### Method:
27 c> Fit a 3x3 window over the random point. The closest
28 c> 2x2 points will surround the point. But based on which
29 c> quadrant of that 2x2 cell in which the point falls, the
30 c> 3x3 window could extend NW, NE, SW or SE from the 2x2 cell.
31  subroutine biquad(z,glamn,glomn,dla,dlo,
32  *nla,nlo,maxla,maxlo,xla,xlo,val)
33 
34 c - Subroutine to perform a 2-D quadratic ("biquadratic")
35 c - interpolation at location "xla,xlo" off of grid "z", whose
36 c - header information is the standard ".b" header
37 c - information of
38 c - glamn = minimum latitude (real*8 decimal degrees)
39 c - glomn = minimum longitude (real*8 decimal degrees)
40 c - dla = latitude spacing (real*8 decimal degrees)
41 c - dlo = longitude spacing (real*8 decimal degrees)
42 c - nla = number of lat rows (integer*4)
43 c - nlo = number of lon cols (integer*4)
44 c - maxla = actual dimensioned size of "z" in rows
45 c - maxlo = actual dimensioned size of "z" in cols
46 c - Further input:
47 c - xla = lat of pt for interpolation (real*8 dec. deg)
48 c - xlo = lon of pt for interpolation (real*8 dec. deg)
49 c - Output:
50 c - val = interpolated value (real*8)
51 c -
52 c - Method:
53 c - Fit a 3x3 window over the random point. The closest
54 c - 2x2 points will surround the point. But based on which
55 c - quadrant of that 2x2 cell in which the point falls, the
56 c - 3x3 window could extend NW, NE, SW or SE from the 2x2 cell.
57 
58 c - Data is assumed real*4
59  implicit real*8 (a-h,o-z)
60  real*4 z(maxla,maxlo)
61 
62 c - Internal use
63  real*4 x,y,fx0,fx1,fx2,qterp
64 
65 c - Find which row should be LEAST when fitting
66 c - a 3x3 window around xla.
67  difla = (xla - glamn)
68  ratla = difla / (dla/2.d0)
69  ila = int(ratla)+1
70  if(mod(ila,2).ne.0)then
71  jla = (ila+1)/2 - 1
72  else
73  jla = (ila )/2
74  endif
75 c - Fix any edge overlaps
76  if(jla.lt.1)jla = 1
77  if(jla.gt.nla-2)jla = nla-2
78 
79 c - Find which column should be LEAST when fitting
80 c - a 3x3 window around xlo.
81  diflo = (xlo - glomn)
82  ratlo = diflo / (dlo/2.d0)
83  ilo = int(ratlo)+1
84  if(mod(ilo,2).ne.0)then
85  jlo = (ilo+1)/2 - 1
86  else
87  jlo = (ilo )/2
88  endif
89 c - Fix any edge overlaps
90  if(jlo.lt.1)jlo = 1
91  if(jlo.gt.nlo-2)jlo = nlo-2
92 
93 c - In the range of 0(westernmost) to
94 c - 2(easternmost) col, where is our
95 c - random lon value? That is, x must
96 c - be real and fall between 0 and 2.
97 
98  x=(xlo-dlo*(jlo-1)-glomn)/dlo
99 
100  if(x.lt.0.0)then
101  val = 1.d30
102  write(6,100)x,val
103 c stop
104  return
105  endif
106  100 format(
107  *'FATAL in biquad: x<0 : ',f20.10,/,
108  *' --> Returning with val = ',f40.1)
109 
110  if(x.gt.2.0)then
111  val = 1.d30
112  write(6,101)x,val
113 c stop
114  return
115  endif
116  101 format(
117  *'FATAL in biquad: x>2 : ',f20.10,/,
118  *' --> Returning with val = ',f40.1)
119 
120 c - In the range of 0(southernmost) to
121 c - 2(northernmost) row, where is our
122 c - random lat value? That is, x must
123 c - be real and fall between 0 and 2.
124 
125  y=(xla-dla*(jla-1)-glamn)/dla
126 
127  if(y.lt.0.0)then
128  val = 1.d30
129  write(6,102)y,val
130 c stop
131  return
132  endif
133  102 format(
134  *'FATAL in biquad: y<0 : ',f20.10,/,
135  *' --> Returning with val = ',f40.1)
136 
137  if(y.gt.2.0)then
138  val = 1.d30
139  write(6,103)y,val
140 c stop
141  return
142  endif
143  103 format(
144  *'FATAL in biquad: y>2 : ',f20.10,/,
145  *' --> Returning with val = ',f40.1)
146 
147 c - Now do the interpolation. First, build a paraboloa
148 c - east-west the southermost row and interpolate to longitude
149 c - "xlo" (at "x" for 0<x<2). Then do it in the middle
150 c - row, then finally the northern row. The last step
151 c - is to fit a parabola north-south at the three previous
152 c - interpolations, but this time to interpolate to
153 c - latitude "xla" (at "y" for 0<y<2). Obviously we
154 c - could reverse the order, doing 3 north-south parabolas
155 c - followed by one east-east and we'd get the same answer.
156 
157  fx0=qterp(x,z(jla ,jlo),z(jla ,jlo+1),z(jla ,jlo+2))
158  fx1=qterp(x,z(jla+1,jlo),z(jla+1,jlo+1),z(jla+1,jlo+2))
159  fx2=qterp(x,z(jla+2,jlo),z(jla+2,jlo+1),z(jla+2,jlo+2))
160  val=dble(qterp(y,fx0,fx1,fx2))
161 
162 c goto 888
163 
164 c write(6,*) ' '
165 c write(6,*) ' ---------------------------------'
166 c write(6,*) ' Inside biquad.f'
167 c write(6,*) ' '
168 
169 c write(6,*) ' Southern Row of 3x3 window: '
170 c write(6,*) ' jla ,jlo ,z(jla ,jlo ) = ',
171 c *jla ,jlo ,z(jla ,jlo )
172 c write(6,*) ' jla ,jlo+1,z(jla ,jlo+1) = ',
173 c *jla ,jlo+1,z(jla ,jlo+1)
174 c write(6,*) ' jla ,jlo+2,z(jla ,jlo+2) = ',
175 c *jla ,jlo+2,z(jla ,jlo+2)
176 c write(6,*) ' fx0 = ',fx0
177 
178 c write(6,*) ' Middle Row of 3x3 window: '
179 c write(6,*) ' jla+1,jlo ,z(jla+1,jlo ) = ',
180 c *jla+1,jlo ,z(jla+1,jlo )
181 c write(6,*) ' jla+1,jlo+1,z(jla+1,jlo+1) = ',
182 c *jla+1,jlo+1,z(jla+1,jlo+1)
183 c write(6,*) ' jla+1,jlo+2,z(jla+1,jlo+2) = ',
184 c *jla+1,jlo+2,z(jla+1,jlo+2)
185 c write(6,*) ' fx1 = ',fx1
186 
187 c write(6,*) ' Northern Row of 3x3 window: '
188 c write(6,*) ' jla+2,jlo ,z(jla+2,jlo ) = ',
189 c *jla+2,jlo ,z(jla+2,jlo )
190 c write(6,*) ' jla+2,jlo+1,z(jla+2,jlo+1) = ',
191 c *jla+2,jlo+1,z(jla+2,jlo+1)
192 c write(6,*) ' jla+2,jlo+2,z(jla+2,jlo+2) = ',
193 c *jla+2,jlo+2,z(jla+2,jlo+2)
194 c write(6,*) ' fx2 = ',fx2
195 
196 c write(6,*) ' Final interpolated value = ',val
197 
198 c 888 continue
199 
200 
201  return
202  end
203  include 'qterp.f'
subroutine biquad(z, glamn, glomn, dla, dlo, nla, nlo, maxla, maxlo, xla, xlo, val)
Subroutine to perform a 2-D quadratic ("biquadratic") interpolation.
Definition: biquad.f:33