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