NADCON5-ng  0.0.2
NADCON5 Next Generation Documentation
regrd2.f
Go to the documentation of this file.
1 c> \ingroup core
2 c> \if MANPAGE
3 c> \page regrd2
4 c> \endif
5 c>
6 c> Part of the NADCON5 \ref core , regrid data
7 c>
8 c> Regrid gridded data using biquadratic interpolation
9 c>
10 c> ### Program arguments
11 c> Arguments are newline terminated and read from standard input
12 c>
13 c> When run from the command line, the program prints a prompt for each argument
14 c>
15 c> They are enumerated here
16 c> \param infile Input Master Grid File Name
17 c> \param outfile Output Regrid File Name
18 c> \param nrow Number of rows in new Grid (Lat)
19 c> \param ncol Number of cols in new Grid (Lon)
20 c>
21 c> ### Program Inputs:
22 c>
23 c> - `lin` Input File
24 c>
25 c> ### Program Outputs:
26 c>
27 c> - `lout` Output File
28 c>
29  program regrd2
30 
31 *** regrid gridded data using biquadratic interpolation
32 
33  implicit double precision(a-h,o-z)
34  parameter(len=6000*40000)
35  real*4 z(len),zrec(40000)
36  integer h(len),hrec(40000)
37  character*68 fname
38  equivalence(h(1),z(1)),(hrec(1),zrec(1))
39  common/gstuff/glamn,dgla,glamx,glomn,dglo,glomx,nla,nlo,nclip
40 
41  common/static/z,zrec
42 
43  lin =1
44  lout=2
45 
46  write(*,*) 'program regrd2'
47  write(*,*)
48 
49  write(*,2)
50  2 format('Enter master grid file: ',$)
51  read(*,1) fname
52  1 format(a)
53  open(lin,file=fname,status='old',form='unformatted')
54  read(lin) glamn,glomn,dgla,dglo,nla,nlo,ikind
55  glamx=glamn+dgla*(nla-1)
56  glomx=glomn+dglo*(nlo-1)
57  if(nla*nlo.gt.len) stop 12345
58  if(nla.lt.4.or.nlo.lt.4) stop 54321
59 
60 *** display header contents
61 
62  write(*,*) 'kind=',ikind
63  write(*,*) 'LAT min=',glamn
64  write(*,*) ' del=',dgla,' # lat=',nla
65  write(*,*) ' max=',glamx
66  write(*,*) 'LON min=',glomn
67  write(*,*) ' del=',dglo,' # lon=',nlo
68  write(*,*) ' max=',glomx
69  write(*,*)
70 
71  write(*,3)
72  3 format('Enter regridded output file:',$)
73  read(*,1) fname
74  open(lout,file=fname,status='new',form='unformatted')
75 
76 *** following disabled
77 
78 * 15 write(*,5)
79 * 5 format('Enter minimum latitude of new grid: ',$)
80 * read(*,*) glamn2
81 * if(glamn2.lt.glamn.or.glamn2.gt.glamx) go to 15
82 
83 * 16 write(*,6)
84 * 6 format('Enter maximum latitude of new grid: ',$)
85 * read(*,*) glamx2
86 * if(glamx2.lt.glamn2.or.glamx2.gt.glamx) go to 16
87 
88 * 17 write(*,7)
89 * 7 format('Enter minimum longitude of new grid: ',$)
90 * read(*,*) glomn2
91 * if(glomn2.lt.glomn.or.glomn2.gt.glomx) go to 17
92 
93 * 18 write(*,8)
94 * 8 format('Enter maximum longitude of new grid: ',$)
95 * read(*,*) glomx2
96 * if(glomx2.lt.glomn2.or.glomx2.gt.glomx) go to 18
97 
98 * write(*,4)
99 * 4 format('Enter delta latitude (in minutes) : ',$)
100 * read(*,*) dglam
101 
102 * write(*,9)
103 * 9 format('Enter delta longitude (in minutes): ',$)
104 * read(*,*) dglom
105 
106 * dgla2=dglam/60.d0
107 * dglo2=dglom/60.d0
108 * nla2=idnint((glamx2-glamn2)/dgla2)+1
109 * nlo2=idnint((glomx2-glomn2)/dglo2)+1
110 
111 *** foregoing disabled
112 *** following clause added
113 
114  93 write(*,94)
115  94 format('Enter new number of rows (latitude) : ',$)
116  read(*,*) nla2
117  if(nla2.lt.0) go to 93
118 
119  98 write(*,99)
120  99 format('Enter new number of columns (longitude): ',$)
121  read(*,*) nlo2
122  if(nlo2.lt.0) go to 98
123 
124  glamn2=glamn
125  glamx2=glamx
126  glomn2=glomn
127  glomx2=glomx
128  dgla2=(glamx2-glamn2)/(nla2-1)
129  dglo2=(glomx2-glomn2)/(nlo2-1)
130 
131 *** back to standard code
132 
133 *** input the grid
134 
135  if(ikind.eq.0) then
136  call r1(lin,h,nla,nlo)
137  else
138  call r2(lin,z,nla,nlo)
139  endif
140  close(lin)
141 
142 *** do the interpolation
143 
144  write(lout) glamn2,glomn2,dgla2,dglo2,nla2,nlo2,ikind
145 
146  if(ikind.eq.0) then
147  do 10 ir=1,nla2
148  gla=glamn2+(ir-1)*dgla2
149  do 20 ic=1,nlo2
150  glo=glomn2+(ic-1)*dglo2
151  call bquad1(h,gla,glo,ival)
152  if(ival.ge.2147483647) stop 11111
153  hrec(ic)=ival
154  20 continue
155  write(lout) (hrec(i),i=1,nlo2)
156  10 continue
157  else
158  do 30 ir=1,nla2
159  gla=glamn2+(ir-1)*dgla2
160  do 40 ic=1,nlo2
161  glo=glomn2+(ic-1)*dglo2
162  call bquad2(z,gla,glo,val)
163  if(val.ge.1.d30) stop 22222
164  zrec(ic)=val
165  40 continue
166  write(lout) (zrec(i),i=1,nlo2)
167  30 continue
168  endif
169 
170  stop
171  end
172  subroutine bquad1(h,gla,glo,ival)
173 
174 *** compute biquadratic interpolation, integer version
175 *** logic not tested for '0/360 meridian wraparound'
176 
177  implicit double precision(a-h,o-z)
178  integer h(nla,nlo)
179  real*4 x,y,fx0,fx1,fx2,val,qterp1,qterp2
180  external qterp1,qterp2
181  common/gstuff/glamn,dgla,glamx,glomn,dglo,glomx,nla,nlo,nclip
182 
183  if(gla.lt.glamn.or.gla.gt.glamx.or.
184  * glo.lt.glomn.or.glo.gt.glomx) then
185  nclip=nclip+1
186  ival=2147483647
187  else
188 
189 *** within grid boundaries, get indicies in the grid
190 
191  ix=(glo-glomn)/dglo+1.d0
192  ix1=ix+1
193  ix2=ix+2
194 
195 *** check if edge collision
196 
197  1 if(ix2.gt.nlo)then
198  ix =ix -1
199  ix1=ix1-1
200  ix2=ix2-1
201  go to 1
202  endif
203 
204  x=(glo-dglo*(ix-1)-glomn)/dglo
205 
206 *** move grid to get nearer center of 3 x 3
207 
208  if(x.lt.0.5.and.ix.gt.1)then
209  ix =ix -1
210  ix1=ix1-1
211  ix2=ix2-1
212  x=x+1.
213  endif
214  if(x.lt.0..or.x.gt.2.) stop 55555
215 
216 *** now do it for y
217 
218  jy=(gla-glamn)/dgla+1.d0
219  jy1=jy+1
220  jy2=jy+2
221  2 if(jy2.gt.nla)then
222  jy =jy -1
223  jy1=jy1-1
224  jy2=jy2-1
225  go to 2
226  endif
227  y=(gla-dgla*(jy-1)-glamn)/dgla
228  if(y.lt.0.5.and.jy.gt.1)then
229  jy =jy -1
230  jy1=jy1-1
231  jy2=jy2-1
232  y=y+1.
233  endif
234  if(y.lt.0..or.y.gt.2.) stop 33333
235 
236  fx0=qterp1(x, h(jy ,ix ), h(jy ,ix1), h(jy ,ix2))
237  fx1=qterp1(x, h(jy1,ix ), h(jy1,ix1), h(jy1,ix2))
238  fx2=qterp1(x, h(jy2,ix ), h(jy2,ix1), h(jy2,ix2))
239  val=qterp2(y, fx0 , fx1 , fx2 )
240  ival=nint(val)
241  endif
242 
243  return
244  end
245  subroutine bquad2(z,gla,glo,val)
246 
247 *** compute biquadratic interpolation, floating point version
248 *** logic not tested for '0/360 meridian wraparound'
249 
250  implicit double precision(a-h,o-z)
251  real*4 z(nla,nlo)
252  real*4 x,y,fx0,fx1,fx2,qterp2
253  external qterp2
254  common/gstuff/glamn,dgla,glamx,glomn,dglo,glomx,nla,nlo,nclip
255 
256  if(gla.lt.glamn.or.gla.gt.glamx.or.
257  * glo.lt.glomn.or.glo.gt.glomx) then
258  nclip=nclip+1
259  val=1.d30
260  else
261 
262 *** within grid boundaries, get indicies in the grid
263 
264  ix=(glo-glomn)/dglo+1.d0
265  ix1=ix+1
266  ix2=ix+2
267 
268 *** check if edge collision
269 
270  1 if(ix2.gt.nlo)then
271  ix =ix -1
272  ix1=ix1-1
273  ix2=ix2-1
274  go to 1
275  endif
276 
277  x=(glo-dglo*(ix-1)-glomn)/dglo
278 
279 *** move grid to get nearer center of 3 x 3
280 
281  if(x.lt.0.5.and.ix.gt.1)then
282  ix =ix -1
283  ix1=ix1-1
284  ix2=ix2-1
285  x=x+1.
286  endif
287  if(x.lt.0..or.x.gt.2.) stop 66666
288 
289 *** now do it for y
290 
291  jy=(gla-glamn)/dgla+1.d0
292  jy1=jy+1
293  jy2=jy+2
294  2 if(jy2.gt.nla)then
295  jy =jy -1
296  jy1=jy1-1
297  jy2=jy2-1
298  go to 2
299  endif
300  y=(gla-dgla*(jy-1)-glamn)/dgla
301  if(y.lt.0.5.and.jy.gt.1)then
302  jy =jy -1
303  jy1=jy1-1
304  jy2=jy2-1
305  y=y+1.
306  endif
307  if(y.lt.0..or.y.gt.2.) stop 44444
308 
309  fx0= qterp2(x, z(jy ,ix ), z(jy ,ix1), z(jy ,ix2))
310  fx1= qterp2(x, z(jy1,ix ), z(jy1,ix1), z(jy1,ix2))
311  fx2= qterp2(x, z(jy2,ix ), z(jy2,ix1), z(jy2,ix2))
312  val=dble(qterp2(y, fx0 , fx1 , fx2 ))
313  endif
314 
315  return
316  end
317  real function qterp1(x,if0,if1,if2)
318 
319 *** linear quadratic interpolation from equally spaced values
320 *** uses newton-gregory forward polynomial
321 *** x ranges from 0 thru 2. (thus s = x)
322 
323 *** forward differences
324 
325  idf0 =if1 -if0
326  idf1 =if2 -if1
327  id2f0=idf1-idf0
328 
329  qterp1=if0 + x*idf0 + 0.5*x*(x-1.)*id2f0
330 
331  return
332  end
333  real function qterp2(x,f0,f1,f2)
334 
335 *** linear quadratic interpolation from equally spaced values
336 *** uses newton-gregory forward polynomial
337 *** x ranges from 0 thru 2. (thus s = x)
338 
339 *** forward differences
340 
341  df0 =f1 -f0
342  df1 =f2 -f1
343  d2f0=df1-df0
344 
345  qterp2=f0 + x*df0 + 0.5*x*(x-1.)*d2f0
346 
347  return
348  end
349  subroutine r1(lin,h,nla,nlo)
350 
351 *** read records south to north (elements are west to east)
352 
353  implicit double precision(a-h,o-z)
354  integer h(nla,nlo)
355 
356  do 1 i=1,nla
357  1 read(lin) (h(i,j),j=1,nlo)
358 
359  return
360  end
361  subroutine r2(lin,z,nla,nlo)
362 
363 *** read records south to north (elements are west to east)
364 
365  implicit double precision(a-h,o-z)
366  real*4 z(nla,nlo)
367 
368  do 1 i=1,nla
369  1 read(lin) (z(i,j),j=1,nlo)
370 
371  return
372  end
bquad2
subroutine bquad2(z, gla, glo, val)
Definition: regrd2.f:246
bquad1
subroutine bquad1(h, gla, glo, ival)
Definition: regrd2.f:173
qterp2
real function qterp2(x, f0, f1, f2)
Definition: regrd2.f:334
qterp1
real function qterp1(x, if0, if1, if2)
Definition: regrd2.f:318
r2
subroutine r2(lin, z, nla, nlo)
Definition: regrd2.f:362
regrd2
program regrd2
Part of the NADCON5 NADCON5 Core Library , regrid data.
Definition: regrd2.f:29
r1
subroutine r1(lin, h, nla, nlo)
Definition: regrd2.f:350