| 5 |
michaesp |
1 |
real function int2d(ar,n1,n2,rid,rjd)
|
|
|
2 |
c-----------------------------------------------------------------------
|
|
|
3 |
c Purpose:
|
|
|
4 |
c This subroutine interpolates a 2d-array to an arbitrary
|
|
|
5 |
c location within the grid.
|
|
|
6 |
c Arguments:
|
|
|
7 |
c ar real input surface pressure, define as ar(n1,n2)
|
|
|
8 |
c n1,n2 int input dimensions of ar
|
|
|
9 |
c ri,rj real input grid location to be interpolated to
|
|
|
10 |
c History:
|
|
|
11 |
c-----------------------------------------------------------------------
|
|
|
12 |
|
|
|
13 |
c argument declarations
|
|
|
14 |
integer n1,n2
|
|
|
15 |
real ar(n1,n2), rid,rjd
|
|
|
16 |
|
|
|
17 |
c local declarations
|
|
|
18 |
integer i,j,ip1,jp1,ih,jh
|
|
|
19 |
real frac0i,frac0j,frac1i,frac1j,ri,rj
|
|
|
20 |
|
|
|
21 |
c do linear interpolation
|
|
|
22 |
ri=amax1(1.,amin1(float(n1),rid))
|
|
|
23 |
rj=amax1(1.,amin1(float(n2),rjd))
|
|
|
24 |
ih=nint(ri)
|
|
|
25 |
jh=nint(rj)
|
|
|
26 |
|
|
|
27 |
c Check for interpolation in i
|
|
|
28 |
* if (abs(float(ih)-ri).lt.1.e-3) then
|
|
|
29 |
* i =ih
|
|
|
30 |
* ip1=ih
|
|
|
31 |
* else
|
|
|
32 |
i =min0(int(ri),n1-1)
|
|
|
33 |
ip1=i+1
|
|
|
34 |
* endif
|
|
|
35 |
|
|
|
36 |
c Check for interpolation in j
|
|
|
37 |
* if (abs(float(jh)-rj).lt.1.e-3) then
|
|
|
38 |
* j =jh
|
|
|
39 |
* jp1=jh
|
|
|
40 |
* else
|
|
|
41 |
j =min0(int(rj),n2-1)
|
|
|
42 |
jp1=j+1
|
|
|
43 |
* endif
|
|
|
44 |
|
|
|
45 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
46 |
c no interpolation at all
|
|
|
47 |
int2d=ar(i,j)
|
|
|
48 |
else
|
|
|
49 |
frac0i=ri-float(i)
|
|
|
50 |
frac0j=rj-float(j)
|
|
|
51 |
frac1i=1.-frac0i
|
|
|
52 |
frac1j=1.-frac0j
|
|
|
53 |
int2d = ar(i ,j ) * frac1i * frac1j
|
|
|
54 |
& + ar(i ,jp1) * frac1i * frac0j
|
|
|
55 |
& + ar(ip1,j ) * frac0i * frac1j
|
|
|
56 |
& + ar(ip1,jp1) * frac0i * frac0j
|
|
|
57 |
endif
|
|
|
58 |
end
|
|
|
59 |
real function int2dm(ar,n1,n2,rid,rjd,misdat)
|
|
|
60 |
c-----------------------------------------------------------------------
|
|
|
61 |
c Purpose:
|
|
|
62 |
c This subroutine interpolates a 2d-array to an arbitrary
|
|
|
63 |
c location within the grid. The interpolation includes the
|
|
|
64 |
c testing of the missing data flag 'misdat'.
|
|
|
65 |
c Arguments:
|
|
|
66 |
c ar real input surface pressure, define as ar(n1,n2)
|
|
|
67 |
c n1,n2 int input dimensions of ar
|
|
|
68 |
c ri,rj real input grid location to be interpolated to
|
|
|
69 |
c misdat real input missing data flag (on if misdat<>0)
|
|
|
70 |
c Warning:
|
|
|
71 |
c This routine has not yet been seriously tested
|
|
|
72 |
c History:
|
|
|
73 |
c-----------------------------------------------------------------------
|
|
|
74 |
|
|
|
75 |
c argument declarations
|
|
|
76 |
integer n1,n2
|
|
|
77 |
real ar(n1,n2), rid,rjd, misdat
|
|
|
78 |
|
|
|
79 |
c local declarations
|
|
|
80 |
integer i,j,ip1,jp1,ih,jh
|
|
|
81 |
real frac0i,frac0j,frac1i,frac1j,ri,rj,int2d
|
|
|
82 |
|
|
|
83 |
c check if routine without missing data checking can be called instead
|
|
|
84 |
if (misdat.eq.0.) then
|
|
|
85 |
int2dm=int2d(ar,n1,n2,rid,rjd)
|
|
|
86 |
return
|
|
|
87 |
endif
|
|
|
88 |
|
|
|
89 |
c do linear interpolation
|
|
|
90 |
ri=amax1(1.,amin1(float(n1),rid))
|
|
|
91 |
rj=amax1(1.,amin1(float(n2),rjd))
|
|
|
92 |
ih=nint(ri)
|
|
|
93 |
jh=nint(rj)
|
|
|
94 |
|
|
|
95 |
c Check for interpolation in i
|
|
|
96 |
* if (abs(float(ih)-ri).lt.1.e-3) then
|
|
|
97 |
* i =ih
|
|
|
98 |
* ip1=ih
|
|
|
99 |
* else
|
|
|
100 |
i =min0(int(ri),n1-1)
|
|
|
101 |
ip1=i+1
|
|
|
102 |
* endif
|
|
|
103 |
|
|
|
104 |
c Check for interpolation in j
|
|
|
105 |
* if (abs(float(jh)-rj).lt.1.e-3) then
|
|
|
106 |
* j =jh
|
|
|
107 |
* jp1=jh
|
|
|
108 |
* else
|
|
|
109 |
j =min0(int(rj),n2-1)
|
|
|
110 |
jp1=j+1
|
|
|
111 |
* endif
|
|
|
112 |
|
|
|
113 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
114 |
c no interpolation at all
|
|
|
115 |
int2dm=ar(i,j)
|
|
|
116 |
else
|
|
|
117 |
if ((misdat.eq.ar(i ,j )).or.
|
|
|
118 |
& (misdat.eq.ar(i ,jp1)).or.
|
|
|
119 |
& (misdat.eq.ar(ip1,j )).or.
|
|
|
120 |
& (misdat.eq.ar(ip1,jp1))) then
|
|
|
121 |
int2dm=misdat
|
|
|
122 |
else
|
|
|
123 |
frac0i=ri-float(i)
|
|
|
124 |
frac0j=rj-float(j)
|
|
|
125 |
frac1i=1.-frac0i
|
|
|
126 |
frac1j=1.-frac0j
|
|
|
127 |
int2dm = ar(i ,j ) * frac1i * frac1j
|
|
|
128 |
& + ar(i ,jp1) * frac1i * frac0j
|
|
|
129 |
& + ar(ip1,j ) * frac0i * frac1j
|
|
|
130 |
& + ar(ip1,jp1) * frac0i * frac0j
|
|
|
131 |
endif
|
|
|
132 |
endif
|
|
|
133 |
end
|
|
|
134 |
real function int2dp(ar,n1,n2,rid,rjd)
|
|
|
135 |
c-----------------------------------------------------------------------
|
|
|
136 |
c Purpose:
|
|
|
137 |
c This subroutine interpolates a 2d-array to an arbitrary
|
|
|
138 |
c location within the grid. The 2d-array is periodic in the
|
|
|
139 |
c i-direction: ar(0,.)=ar(n1,.) and ar(1,.)=ar(n1+1,.).
|
|
|
140 |
c Therefore rid can take values in the range 0,...,n1+1
|
|
|
141 |
c Arguments:
|
|
|
142 |
c ar real input surface pressure, define as ar(n1,n2)
|
|
|
143 |
c n1,n2 int input dimensions of ar
|
|
|
144 |
c ri,rj real input grid location to be interpolated to
|
|
|
145 |
c History:
|
|
|
146 |
c-----------------------------------------------------------------------
|
|
|
147 |
|
|
|
148 |
c argument declarations
|
|
|
149 |
integer n1,n2
|
|
|
150 |
real ar(0:n1+1,n2), rid,rjd
|
|
|
151 |
|
|
|
152 |
c local declarations
|
|
|
153 |
integer i,j,ip1,jp1,ih,jh
|
|
|
154 |
real frac0i,frac0j,frac1i,frac1j,ri,rj
|
|
|
155 |
|
|
|
156 |
c do linear interpolation
|
|
|
157 |
ri=amax1(0.,amin1(float(n1+1),rid))
|
|
|
158 |
rj=amax1(1.,amin1(float(n2),rjd))
|
|
|
159 |
ih=nint(ri)
|
|
|
160 |
jh=nint(rj)
|
|
|
161 |
|
|
|
162 |
c Check for interpolation in i
|
|
|
163 |
* if (abs(float(ih)-ri).lt.1.e-5) then
|
|
|
164 |
* i =ih
|
|
|
165 |
* ip1=ih
|
|
|
166 |
* else
|
|
|
167 |
i =min0(int(ri),n1-1)
|
|
|
168 |
ip1=i+1
|
|
|
169 |
* endif
|
|
|
170 |
|
|
|
171 |
c Check for interpolation in j
|
|
|
172 |
* if (abs(float(jh)-rj).lt.1.e-5) then
|
|
|
173 |
* j =jh
|
|
|
174 |
* jp1=jh
|
|
|
175 |
* else
|
|
|
176 |
j =min0(int(rj),n2-1)
|
|
|
177 |
jp1=j+1
|
|
|
178 |
* endif
|
|
|
179 |
|
|
|
180 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
181 |
c no interpolation at all
|
|
|
182 |
int2dp=ar(i,j)
|
|
|
183 |
else
|
|
|
184 |
frac0i=ri-float(i)
|
|
|
185 |
frac0j=rj-float(j)
|
|
|
186 |
frac1i=1.-frac0i
|
|
|
187 |
frac1j=1.-frac0j
|
|
|
188 |
int2dp = ar(i ,j ) * frac1i * frac1j
|
|
|
189 |
& + ar(i ,jp1) * frac1i * frac0j
|
|
|
190 |
& + ar(ip1,j ) * frac0i * frac1j
|
|
|
191 |
& + ar(ip1,jp1) * frac0i * frac0j
|
|
|
192 |
endif
|
|
|
193 |
end
|
|
|
194 |
real function cint2d(ar,n1,n2,n3,rid,rjd,ikd)
|
|
|
195 |
c-----------------------------------------------------------------------
|
|
|
196 |
c Purpose:
|
|
|
197 |
c This subroutine interpolates a 3d-array to an arbitrary
|
|
|
198 |
c location in the horizontal. ikd specifies the level (must
|
|
|
199 |
c be an integer). A bicubic method is applied (following
|
|
|
200 |
c the Numerical Recipes).
|
|
|
201 |
c Arguments:
|
|
|
202 |
c ar real input field, define as ar(n1,n2,n3)
|
|
|
203 |
c n1,n2 int input dimensions of ar
|
|
|
204 |
c rid,rjd real input grid location to be interpolated to
|
|
|
205 |
c ikd int input level for interpolation
|
|
|
206 |
c History:
|
|
|
207 |
c-----------------------------------------------------------------------
|
|
|
208 |
|
|
|
209 |
c argument declarations
|
|
|
210 |
integer n1,n2,n3,ikd
|
|
|
211 |
real ar(n1,n2,n3), rid,rjd
|
|
|
212 |
|
|
|
213 |
c local declarations
|
|
|
214 |
integer i,j,k
|
|
|
215 |
real y(4),y1(4),y2(4),y12(4)
|
|
|
216 |
|
|
|
217 |
c indices of lower left corner of interpolation grid
|
|
|
218 |
|
|
|
219 |
i=int(rid)
|
|
|
220 |
j=int(rjd)
|
|
|
221 |
k=ikd
|
|
|
222 |
|
|
|
223 |
c do not interpolate if i or j are at the boundaries
|
|
|
224 |
|
|
|
225 |
if ((i.eq.1).or.(j.eq.1).or.(i.eq.n1).or.(j.eq.n2)) then
|
|
|
226 |
cint2d=ar(i,j,k)
|
|
|
227 |
return
|
|
|
228 |
endif
|
|
|
229 |
|
|
|
230 |
c define the arrays y,y1,y2,y12 necessary for the bicubic
|
|
|
231 |
c interpolation (following the Numerical Recipes).
|
|
|
232 |
|
|
|
233 |
y(1)=ar(i,j,k)
|
|
|
234 |
y(2)=ar(i+1,j,k)
|
|
|
235 |
y(3)=ar(i+1,j+1,k)
|
|
|
236 |
y(4)=ar(i,j+1,k)
|
|
|
237 |
|
|
|
238 |
y1(1)=-(ar(i-1,j,k)-ar(i+1,j,k))/2.
|
|
|
239 |
y1(2)=-(ar(i,j,k)-ar(i+2,j,k))/2.
|
|
|
240 |
y1(3)=-(ar(i,j+1,k)-ar(i+2,j+1,k))/2.
|
|
|
241 |
y1(4)=-(ar(i-1,j+1,k)-ar(i+1,j+1,k))/2.
|
|
|
242 |
|
|
|
243 |
y2(1)=-(ar(i,j-1,k)-ar(i,j+1,k))/2.
|
|
|
244 |
y2(2)=-(ar(i+1,j-1,k)-ar(i+1,j+1,k))/2.
|
|
|
245 |
y2(3)=-(ar(i+1,j,k)-ar(i+1,j+2,k))/2.
|
|
|
246 |
y2(4)=-(ar(i,j,k)-ar(i,j+2,k))/2.
|
|
|
247 |
|
|
|
248 |
y12(1)=(ar(i+1,j+1,k)-ar(i+1,j-1,k)-ar(i-1,j+1,k)+ar(i-1,j-1,k))/4.
|
|
|
249 |
y12(2)=(ar(i+2,j+1,k)-ar(i+2,j-1,k)-ar(i,j+1,k)+ar(i,j-1,k))/4.
|
|
|
250 |
y12(3)=(ar(i+2,j+2,k)-ar(i+2,j,k)-ar(i,j+2,k)+ar(i,j,k))/4.
|
|
|
251 |
y12(4)=(ar(i+1,j+2,k)-ar(i+1,j,k)-ar(i-1,j+2,k)+ar(i-1,j,k))/4.
|
|
|
252 |
|
|
|
253 |
call bcuint(y,y1,y2,y12,i,j,rid,rjd,cint2d)
|
|
|
254 |
return
|
|
|
255 |
end
|
|
|
256 |
real function int3d(ar,n1,n2,n3,rid,rjd,rkd)
|
|
|
257 |
c-----------------------------------------------------------------------
|
|
|
258 |
c Purpose:
|
|
|
259 |
c This subroutine interpolates a 3d-array to an arbitrary
|
|
|
260 |
c location within the grid.
|
|
|
261 |
c Arguments:
|
|
|
262 |
c ar real input surface pressure, define as ar(n1,n2,n3)
|
|
|
263 |
c n1,n2,n3 int input dimensions of ar
|
|
|
264 |
c ri,rj,rk real input grid location to be interpolated to
|
|
|
265 |
c History:
|
|
|
266 |
c-----------------------------------------------------------------------
|
|
|
267 |
|
|
|
268 |
c argument declarations
|
|
|
269 |
integer n1,n2,n3
|
|
|
270 |
real ar(n1,n2,n3), rid,rjd,rkd
|
|
|
271 |
|
|
|
272 |
c local declarations
|
|
|
273 |
integer i,j,k,ip1,jp1,kp1,ih,jh,kh
|
|
|
274 |
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk
|
|
|
275 |
|
|
|
276 |
c do linear interpolation
|
|
|
277 |
ri=amax1(1.,amin1(float(n1),rid))
|
|
|
278 |
rj=amax1(1.,amin1(float(n2),rjd))
|
|
|
279 |
rk=amax1(1.,amin1(float(n3),rkd))
|
|
|
280 |
ih=nint(ri)
|
|
|
281 |
jh=nint(rj)
|
|
|
282 |
kh=nint(rk)
|
|
|
283 |
|
|
|
284 |
c Check for interpolation in i
|
|
|
285 |
* if (abs(float(ih)-ri).lt.1.e-3) then
|
|
|
286 |
* i =ih
|
|
|
287 |
* ip1=ih
|
|
|
288 |
* else
|
|
|
289 |
i =min0(int(ri),n1-1)
|
|
|
290 |
ip1=i+1
|
|
|
291 |
* endif
|
|
|
292 |
|
|
|
293 |
c Check for interpolation in j
|
|
|
294 |
if (abs(float(jh)-rj).lt.1.e-3) then
|
|
|
295 |
j =jh
|
|
|
296 |
jp1=jh
|
|
|
297 |
else
|
|
|
298 |
j =min0(int(rj),n2-1)
|
|
|
299 |
jp1=j+1
|
|
|
300 |
endif
|
|
|
301 |
|
|
|
302 |
c Check for interpolation in k
|
|
|
303 |
* if (abs(float(kh)-rk).lt.1.e-3) then
|
|
|
304 |
* k =kh
|
|
|
305 |
* kp1=kh
|
|
|
306 |
* else
|
|
|
307 |
k =min0(int(rk),n3-1)
|
|
|
308 |
kp1=k+1
|
|
|
309 |
* endif
|
|
|
310 |
|
|
|
311 |
if (k.eq.kp1) then
|
|
|
312 |
c no interpolation in k
|
|
|
313 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
314 |
c no interpolation at all
|
|
|
315 |
int3d=ar(i,j,k)
|
|
|
316 |
c print *,'int3d 00: ',rid,rjd,rkd,int3d
|
|
|
317 |
else
|
|
|
318 |
c horizontal interpolation only
|
|
|
319 |
frac0i=ri-float(i)
|
|
|
320 |
frac0j=rj-float(j)
|
|
|
321 |
frac1i=1.-frac0i
|
|
|
322 |
frac1j=1.-frac0j
|
|
|
323 |
int3d = ar(i ,j ,k ) * frac1i * frac1j
|
|
|
324 |
& + ar(i ,jp1,k ) * frac1i * frac0j
|
|
|
325 |
& + ar(ip1,j ,k ) * frac0i * frac1j
|
|
|
326 |
& + ar(ip1,jp1,k ) * frac0i * frac0j
|
|
|
327 |
c print *,'int3d 10: ',rid,rjd,rkd,int3d
|
|
|
328 |
endif
|
|
|
329 |
else
|
|
|
330 |
frac0k=rk-float(k)
|
|
|
331 |
frac1k=1.-frac0k
|
|
|
332 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
333 |
c vertical interpolation only
|
|
|
334 |
int3d = ar(i ,j ,k ) * frac1k
|
|
|
335 |
& + ar(i ,j ,kp1) * frac0k
|
|
|
336 |
c print *,'int3d 01: ',rid,rjd,rkd,int3d
|
|
|
337 |
else
|
|
|
338 |
c full 3d interpolation
|
|
|
339 |
frac0i=ri-float(i)
|
|
|
340 |
frac0j=rj-float(j)
|
|
|
341 |
frac1i=1.-frac0i
|
|
|
342 |
frac1j=1.-frac0j
|
|
|
343 |
int3d = ar(i ,j ,k ) * frac1i * frac1j * frac1k
|
|
|
344 |
& + ar(i ,jp1,k ) * frac1i * frac0j * frac1k
|
|
|
345 |
& + ar(ip1,j ,k ) * frac0i * frac1j * frac1k
|
|
|
346 |
& + ar(ip1,jp1,k ) * frac0i * frac0j * frac1k
|
|
|
347 |
& + ar(i ,j ,kp1) * frac1i * frac1j * frac0k
|
|
|
348 |
& + ar(i ,jp1,kp1) * frac1i * frac0j * frac0k
|
|
|
349 |
& + ar(ip1,j ,kp1) * frac0i * frac1j * frac0k
|
|
|
350 |
& + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
|
|
|
351 |
c print *,'int3d 11: ',rid,rjd,rkd,int3d
|
|
|
352 |
endif
|
|
|
353 |
endif
|
|
|
354 |
end
|
|
|
355 |
real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
356 |
c-----------------------------------------------------------------------
|
|
|
357 |
c Purpose:
|
|
|
358 |
c This subroutine interpolates a 3d-array to an arbitrary
|
|
|
359 |
c location within the grid. The interpolation includes the
|
|
|
360 |
c testing of the missing data flag 'misdat'.
|
|
|
361 |
c Arguments:
|
|
|
362 |
c ar real input surface pressure, define as ar(n1,n2,n3)
|
|
|
363 |
c n1,n2,n3 int input dimensions of ar
|
|
|
364 |
c ri,rj,rk real input grid location to be interpolated to
|
|
|
365 |
c misdat real input missing data flag (on if misdat<>0)
|
|
|
366 |
c Warning:
|
|
|
367 |
c This routine has not yet been seriously tested
|
|
|
368 |
c History:
|
|
|
369 |
c-----------------------------------------------------------------------
|
|
|
370 |
|
|
|
371 |
c argument declarations
|
|
|
372 |
integer n1,n2,n3
|
|
|
373 |
real ar(n1,n2,n3), rid,rjd,rkd, misdat
|
|
|
374 |
|
|
|
375 |
c local declarations
|
|
|
376 |
integer i,j,k,ip1,jp1,kp1,ih,jh,kh
|
|
|
377 |
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d
|
|
|
378 |
|
|
|
379 |
c check if routine without missing data checking can be called instead
|
|
|
380 |
if (misdat.eq.0.) then
|
|
|
381 |
int3dm=int3d(ar,n1,n2,n3,rid,rjd,rkd)
|
|
|
382 |
return
|
|
|
383 |
endif
|
|
|
384 |
|
|
|
385 |
c do linear interpolation
|
|
|
386 |
ri=amax1(1.,amin1(float(n1),rid))
|
|
|
387 |
rj=amax1(1.,amin1(float(n2),rjd))
|
|
|
388 |
rk=amax1(1.,amin1(float(n3),rkd))
|
|
|
389 |
ih=nint(ri)
|
|
|
390 |
jh=nint(rj)
|
|
|
391 |
kh=nint(rk)
|
|
|
392 |
|
|
|
393 |
c Check for interpolation in i
|
|
|
394 |
* if (abs(float(ih)-ri).lt.1.e-3) then
|
|
|
395 |
* i =ih
|
|
|
396 |
* ip1=ih
|
|
|
397 |
* else
|
|
|
398 |
i =min0(int(ri),n1-1)
|
|
|
399 |
ip1=i+1
|
|
|
400 |
* endif
|
|
|
401 |
|
|
|
402 |
c Check for interpolation in j
|
|
|
403 |
* if (abs(float(jh)-rj).lt.1.e-3) then
|
|
|
404 |
* j =jh
|
|
|
405 |
* jp1=jh
|
|
|
406 |
* else
|
|
|
407 |
j =min0(int(rj),n2-1)
|
|
|
408 |
jp1=j+1
|
|
|
409 |
* endif
|
|
|
410 |
|
|
|
411 |
c Check for interpolation in k
|
|
|
412 |
* if (abs(float(kh)-rk).lt.1.e-3) then
|
|
|
413 |
* k =kh
|
|
|
414 |
* kp1=kh
|
|
|
415 |
* else
|
|
|
416 |
k =min0(int(rk),n3-1)
|
|
|
417 |
kp1=k+1
|
|
|
418 |
* endif
|
|
|
419 |
|
|
|
420 |
if (k.eq.kp1) then
|
|
|
421 |
c no interpolation in k
|
|
|
422 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
423 |
c no interpolation at all
|
|
|
424 |
if (misdat.eq.ar(i,j,k)) then
|
|
|
425 |
int3dm=misdat
|
|
|
426 |
else
|
|
|
427 |
int3dm=ar(i,j,k)
|
|
|
428 |
endif
|
|
|
429 |
c print *,'int3dm 00: ',rid,rjd,rkd,int3dm
|
|
|
430 |
else
|
|
|
431 |
c horizontal interpolation only
|
|
|
432 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
433 |
& (misdat.eq.ar(i ,jp1,k )).or.
|
|
|
434 |
& (misdat.eq.ar(ip1,j ,k )).or.
|
|
|
435 |
& (misdat.eq.ar(ip1,jp1,k ))) then
|
|
|
436 |
int3dm=misdat
|
|
|
437 |
else
|
|
|
438 |
frac0i=ri-float(i)
|
|
|
439 |
frac0j=rj-float(j)
|
|
|
440 |
frac1i=1.-frac0i
|
|
|
441 |
frac1j=1.-frac0j
|
|
|
442 |
int3dm = ar(i ,j ,k ) * frac1i * frac1j
|
|
|
443 |
& + ar(i ,jp1,k ) * frac1i * frac0j
|
|
|
444 |
& + ar(ip1,j ,k ) * frac0i * frac1j
|
|
|
445 |
& + ar(ip1,jp1,k ) * frac0i * frac0j
|
|
|
446 |
c print *,'int3dm 10: ',rid,rjd,rkd,int3dm
|
|
|
447 |
endif
|
|
|
448 |
endif
|
|
|
449 |
else
|
|
|
450 |
frac0k=rk-float(k)
|
|
|
451 |
frac1k=1.-frac0k
|
|
|
452 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
453 |
c vertical interpolation only
|
|
|
454 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
455 |
& (misdat.eq.ar(i ,j ,kp1))) then
|
|
|
456 |
int3dm=misdat
|
|
|
457 |
else
|
|
|
458 |
int3dm = ar(i ,j ,k ) * frac1k
|
|
|
459 |
& + ar(i ,j ,kp1) * frac0k
|
|
|
460 |
c print *,'int3dm 01: ',rid,rjd,rkd,int3dm
|
|
|
461 |
endif
|
|
|
462 |
else
|
|
|
463 |
c full 3d interpolation
|
|
|
464 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
465 |
& (misdat.eq.ar(i ,jp1,k )).or.
|
|
|
466 |
& (misdat.eq.ar(ip1,j ,k )).or.
|
|
|
467 |
& (misdat.eq.ar(ip1,jp1,k )).or.
|
|
|
468 |
& (misdat.eq.ar(i ,j ,kp1)).or.
|
|
|
469 |
& (misdat.eq.ar(i ,jp1,kp1)).or.
|
|
|
470 |
& (misdat.eq.ar(ip1,j ,kp1)).or.
|
|
|
471 |
& (misdat.eq.ar(ip1,jp1,kp1))) then
|
|
|
472 |
int3dm=misdat
|
|
|
473 |
else
|
|
|
474 |
frac0i=ri-float(i)
|
|
|
475 |
frac0j=rj-float(j)
|
|
|
476 |
frac1i=1.-frac0i
|
|
|
477 |
frac1j=1.-frac0j
|
|
|
478 |
int3dm = ar(i ,j ,k ) * frac1i * frac1j * frac1k
|
|
|
479 |
& + ar(i ,jp1,k ) * frac1i * frac0j * frac1k
|
|
|
480 |
& + ar(ip1,j ,k ) * frac0i * frac1j * frac1k
|
|
|
481 |
& + ar(ip1,jp1,k ) * frac0i * frac0j * frac1k
|
|
|
482 |
& + ar(i ,j ,kp1) * frac1i * frac1j * frac0k
|
|
|
483 |
& + ar(i ,jp1,kp1) * frac1i * frac0j * frac0k
|
|
|
484 |
& + ar(ip1,j ,kp1) * frac0i * frac1j * frac0k
|
|
|
485 |
& + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
|
|
|
486 |
c print *,'int3dm 11: ',rid,rjd,rkd,int3dm
|
|
|
487 |
endif
|
|
|
488 |
endif
|
|
|
489 |
endif
|
|
|
490 |
end
|
|
|
491 |
real function int3dmlog(ar,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
492 |
c-----------------------------------------------------------------------
|
|
|
493 |
c Purpose:
|
|
|
494 |
c This subroutine interpolates a 3d-array to an arbitrary
|
|
|
495 |
c location within the grid. The interpolation includes the
|
|
|
496 |
c testing of the missing data flag 'misdat'.
|
|
|
497 |
c Prior to vertical interpolations the log is taken from the array.
|
|
|
498 |
c Arguments:
|
|
|
499 |
c ar real input surface pressure, define as ar(n1,n2,n3)
|
|
|
500 |
c n1,n2,n3 int input dimensions of ar
|
|
|
501 |
c ri,rj,rk real input grid location to be interpolated to
|
|
|
502 |
c misdat real input missing data flag (on if misdat<>0)
|
|
|
503 |
c Warning:
|
|
|
504 |
c This routine has not yet been seriously tested
|
|
|
505 |
c History:
|
|
|
506 |
c-----------------------------------------------------------------------
|
|
|
507 |
|
|
|
508 |
c argument declarations
|
|
|
509 |
integer n1,n2,n3
|
|
|
510 |
real ar(n1,n2,n3), rid,rjd,rkd, misdat
|
|
|
511 |
|
|
|
512 |
c local declarations
|
|
|
513 |
integer i,j,k,ip1,jp1,kp1,ih,jh,kh
|
|
|
514 |
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d
|
|
|
515 |
|
|
|
516 |
c print*,'hallo in SR int3dmlog'
|
|
|
517 |
c check if routine without missing data checking can be called instead
|
|
|
518 |
if (misdat.eq.0.) then
|
|
|
519 |
int3dmlog=int3d(ar,n1,n2,n3,rid,rjd,rkd)
|
|
|
520 |
return
|
|
|
521 |
endif
|
|
|
522 |
|
|
|
523 |
c do linear interpolation
|
|
|
524 |
ri=amax1(1.,amin1(float(n1),rid))
|
|
|
525 |
rj=amax1(1.,amin1(float(n2),rjd))
|
|
|
526 |
rk=amax1(1.,amin1(float(n3),rkd))
|
|
|
527 |
ih=nint(ri)
|
|
|
528 |
jh=nint(rj)
|
|
|
529 |
kh=nint(rk)
|
|
|
530 |
|
|
|
531 |
c Check for interpolation in i
|
|
|
532 |
* if (abs(float(ih)-ri).lt.1.e-3) then
|
|
|
533 |
* i =ih
|
|
|
534 |
* ip1=ih
|
|
|
535 |
* else
|
|
|
536 |
i =min0(int(ri),n1-1)
|
|
|
537 |
ip1=i+1
|
|
|
538 |
* endif
|
|
|
539 |
|
|
|
540 |
c Check for interpolation in j
|
|
|
541 |
* if (abs(float(jh)-rj).lt.1.e-3) then
|
|
|
542 |
* j =jh
|
|
|
543 |
* jp1=jh
|
|
|
544 |
* else
|
|
|
545 |
j =min0(int(rj),n2-1)
|
|
|
546 |
jp1=j+1
|
|
|
547 |
* endif
|
|
|
548 |
|
|
|
549 |
c Check for interpolation in k
|
|
|
550 |
* if (abs(float(kh)-rk).lt.1.e-3) then
|
|
|
551 |
* k =kh
|
|
|
552 |
* kp1=kh
|
|
|
553 |
* else
|
|
|
554 |
k =min0(int(rk),n3-1)
|
|
|
555 |
kp1=k+1
|
|
|
556 |
* endif
|
|
|
557 |
|
|
|
558 |
if (k.eq.kp1) then
|
|
|
559 |
c no interpolation in k
|
|
|
560 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
561 |
c no interpolation at all
|
|
|
562 |
if (misdat.eq.ar(i,j,k)) then
|
|
|
563 |
int3dmlog=misdat
|
|
|
564 |
else
|
|
|
565 |
int3dmlog=ar(i,j,k)
|
|
|
566 |
endif
|
|
|
567 |
c print *,'int3dmlog 00: ',rid,rjd,rkd,int3dmlog
|
|
|
568 |
else
|
|
|
569 |
c horizontal interpolation only
|
|
|
570 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
571 |
& (misdat.eq.ar(i ,jp1,k )).or.
|
|
|
572 |
& (misdat.eq.ar(ip1,j ,k )).or.
|
|
|
573 |
& (misdat.eq.ar(ip1,jp1,k ))) then
|
|
|
574 |
int3dmlog=misdat
|
|
|
575 |
else
|
|
|
576 |
frac0i=ri-float(i)
|
|
|
577 |
frac0j=rj-float(j)
|
|
|
578 |
frac1i=1.-frac0i
|
|
|
579 |
frac1j=1.-frac0j
|
|
|
580 |
int3dmlog = ar(i ,j ,k ) * frac1i * frac1j
|
|
|
581 |
& + ar(i ,jp1,k ) * frac1i * frac0j
|
|
|
582 |
& + ar(ip1,j ,k ) * frac0i * frac1j
|
|
|
583 |
& + ar(ip1,jp1,k ) * frac0i * frac0j
|
|
|
584 |
c print *,'int3dmlog 10: ',rid,rjd,rkd,int3dmlog
|
|
|
585 |
endif
|
|
|
586 |
endif
|
|
|
587 |
else
|
|
|
588 |
frac0k=rk-float(k)
|
|
|
589 |
frac1k=1.-frac0k
|
|
|
590 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
591 |
c vertical interpolation only
|
|
|
592 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
593 |
& (misdat.eq.ar(i ,j ,kp1))) then
|
|
|
594 |
int3dmlog=misdat
|
|
|
595 |
else
|
|
|
596 |
int3dmlog = log(ar(i ,j ,k )) * frac1k
|
|
|
597 |
& + log(ar(i ,j ,kp1)) * frac0k
|
|
|
598 |
int3dmlog = exp(int3dmlog)
|
|
|
599 |
c print *,'int3dmlog 01: ',rid,rjd,rkd,int3dmlog
|
|
|
600 |
endif
|
|
|
601 |
else
|
|
|
602 |
c full 3d interpolation
|
|
|
603 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
604 |
& (misdat.eq.ar(i ,jp1,k )).or.
|
|
|
605 |
& (misdat.eq.ar(ip1,j ,k )).or.
|
|
|
606 |
& (misdat.eq.ar(ip1,jp1,k )).or.
|
|
|
607 |
& (misdat.eq.ar(i ,j ,kp1)).or.
|
|
|
608 |
& (misdat.eq.ar(i ,jp1,kp1)).or.
|
|
|
609 |
& (misdat.eq.ar(ip1,j ,kp1)).or.
|
|
|
610 |
& (misdat.eq.ar(ip1,jp1,kp1))) then
|
|
|
611 |
int3dmlog=misdat
|
|
|
612 |
else
|
|
|
613 |
frac0i=ri-float(i)
|
|
|
614 |
frac0j=rj-float(j)
|
|
|
615 |
frac1i=1.-frac0i
|
|
|
616 |
frac1j=1.-frac0j
|
|
|
617 |
int3dmlog = log(ar(i ,j ,k )) * frac1i * frac1j * frac1k
|
|
|
618 |
& + log(ar(i ,jp1,k )) * frac1i * frac0j * frac1k
|
|
|
619 |
& + log(ar(ip1,j ,k )) * frac0i * frac1j * frac1k
|
|
|
620 |
& + log(ar(ip1,jp1,k )) * frac0i * frac0j * frac1k
|
|
|
621 |
& + log(ar(i ,j ,kp1)) * frac1i * frac1j * frac0k
|
|
|
622 |
& + log(ar(i ,jp1,kp1)) * frac1i * frac0j * frac0k
|
|
|
623 |
& + log(ar(ip1,j ,kp1)) * frac0i * frac1j * frac0k
|
|
|
624 |
& + log(ar(ip1,jp1,kp1)) * frac0i * frac0j * frac0k
|
|
|
625 |
int3dmlog = exp(int3dmlog)
|
|
|
626 |
c print *,'int3dmlog 11: ',rid,rjd,rkd,int3dmlog
|
|
|
627 |
endif
|
|
|
628 |
endif
|
|
|
629 |
endif
|
|
|
630 |
end
|
|
|
631 |
|
|
|
632 |
real function int3dmvc(ar,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
633 |
c-----------------------------------------------------------------------
|
|
|
634 |
c Purpose:
|
|
|
635 |
c This subroutine interpolates a 3d-array to an arbitrary
|
|
|
636 |
c location within the grid. The interpolation includes the
|
|
|
637 |
c testing of the missing data flag 'misdat'.
|
|
|
638 |
c In the vertical a Lagrangian cubic interpolation is
|
|
|
639 |
c performed.
|
|
|
640 |
c Arguments:
|
|
|
641 |
c ar real input surface pressure, define as ar(n1,n2,n3)
|
|
|
642 |
c n1,n2,n3 int input dimensions of ar
|
|
|
643 |
c ri,rj,rk real input grid location to be interpolated to
|
|
|
644 |
c misdat real input missing data flag (on if misdat<>0)
|
|
|
645 |
c Warning:
|
|
|
646 |
c This routine has not yet been seriously tested
|
|
|
647 |
c History:
|
|
|
648 |
c-----------------------------------------------------------------------
|
|
|
649 |
|
|
|
650 |
c argument declarations
|
|
|
651 |
integer n1,n2,n3
|
|
|
652 |
real ar(n1,n2,n3), rid,rjd,rkd, misdat
|
|
|
653 |
|
|
|
654 |
c local declarations
|
|
|
655 |
integer i,j,k,ip1,jp1,kp1,ih,jh,kh,klow,n
|
|
|
656 |
real frac0i,frac0j,frac1i,frac1j,ri,rj,rk
|
|
|
657 |
real int2d(4)
|
|
|
658 |
|
|
|
659 |
real int3dm
|
|
|
660 |
|
|
|
661 |
c if n3 < 4 then do call linear interpolation in the vertical
|
|
|
662 |
if (n3.lt.4) then
|
|
|
663 |
int3dmvc=int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
664 |
return
|
|
|
665 |
endif
|
|
|
666 |
|
|
|
667 |
c do linear interpolation in the horizontal
|
|
|
668 |
ri=amax1(1.,amin1(float(n1),rid))
|
|
|
669 |
rj=amax1(1.,amin1(float(n2),rjd))
|
|
|
670 |
rk=amax1(1.,amin1(float(n3),rkd))
|
|
|
671 |
ih=nint(ri)
|
|
|
672 |
jh=nint(rj)
|
|
|
673 |
kh=nint(rk)
|
|
|
674 |
|
|
|
675 |
c Check for interpolation in i
|
|
|
676 |
* if (abs(float(ih)-ri).lt.1.e-3) then
|
|
|
677 |
* i =ih
|
|
|
678 |
* ip1=ih
|
|
|
679 |
* else
|
|
|
680 |
i =min0(int(ri),n1-1)
|
|
|
681 |
ip1=i+1
|
|
|
682 |
* endif
|
|
|
683 |
|
|
|
684 |
c Check for interpolation in j
|
|
|
685 |
* if (abs(float(jh)-rj).lt.1.e-3) then
|
|
|
686 |
* j =jh
|
|
|
687 |
* jp1=jh
|
|
|
688 |
* else
|
|
|
689 |
j =min0(int(rj),n2-1)
|
|
|
690 |
jp1=j+1
|
|
|
691 |
* endif
|
|
|
692 |
|
|
|
693 |
c Check for interpolation in k
|
|
|
694 |
* if (abs(float(kh)-rk).lt.1.e-3) then
|
|
|
695 |
* k =kh
|
|
|
696 |
* kp1=kh
|
|
|
697 |
* else
|
|
|
698 |
k =min0(int(rk),n3-1)
|
|
|
699 |
kp1=k+1
|
|
|
700 |
* endif
|
|
|
701 |
|
|
|
702 |
if (k.eq.kp1) then
|
|
|
703 |
c no interpolation in k
|
|
|
704 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
705 |
c no interpolation at all
|
|
|
706 |
int3dmvc=ar(i,j,k)
|
|
|
707 |
c print *,'int3dmvc 00: ',rid,rjd,rkd,int3dmvc
|
|
|
708 |
else
|
|
|
709 |
c horizontal interpolation only
|
|
|
710 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
711 |
& (misdat.eq.ar(i ,jp1,k )).or.
|
|
|
712 |
& (misdat.eq.ar(ip1,j ,k )).or.
|
|
|
713 |
& (misdat.eq.ar(ip1,jp1,k ))) then
|
|
|
714 |
int3dmvc=misdat
|
|
|
715 |
else
|
|
|
716 |
frac0i=ri-float(i)
|
|
|
717 |
frac0j=rj-float(j)
|
|
|
718 |
frac1i=1.-frac0i
|
|
|
719 |
frac1j=1.-frac0j
|
|
|
720 |
int3dmvc = ar(i ,j ,k ) * frac1i * frac1j
|
|
|
721 |
& + ar(i ,jp1,k ) * frac1i * frac0j
|
|
|
722 |
& + ar(ip1,j ,k ) * frac0i * frac1j
|
|
|
723 |
& + ar(ip1,jp1,k ) * frac0i * frac0j
|
|
|
724 |
c print *,'int3dmvc 10: ',rid,rjd,rkd,int3dmvc
|
|
|
725 |
endif
|
|
|
726 |
endif
|
|
|
727 |
else
|
|
|
728 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
729 |
c vertical interpolation only
|
|
|
730 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
731 |
& (misdat.eq.ar(i ,j ,kp1))) then
|
|
|
732 |
int3dmvc=misdat
|
|
|
733 |
else
|
|
|
734 |
if (k-1.lt.1) then
|
|
|
735 |
klow=1
|
|
|
736 |
else if (k+2.gt.n3) then
|
|
|
737 |
klow=n3-3
|
|
|
738 |
else
|
|
|
739 |
klow=k-1
|
|
|
740 |
endif
|
|
|
741 |
call cubint(ar(i,j,klow),ar(i,j,klow+1),ar(i,j,klow+2),
|
|
|
742 |
& ar(i,j,klow+3),klow,rk,int3dmvc)
|
|
|
743 |
endif
|
|
|
744 |
else
|
|
|
745 |
c full 3d interpolation
|
|
|
746 |
if ((misdat.eq.ar(i ,j ,k )).or.
|
|
|
747 |
& (misdat.eq.ar(i ,jp1,k )).or.
|
|
|
748 |
& (misdat.eq.ar(ip1,j ,k )).or.
|
|
|
749 |
& (misdat.eq.ar(ip1,jp1,k )).or.
|
|
|
750 |
& (misdat.eq.ar(i ,j ,kp1)).or.
|
|
|
751 |
& (misdat.eq.ar(i ,jp1,kp1)).or.
|
|
|
752 |
& (misdat.eq.ar(ip1,j ,kp1)).or.
|
|
|
753 |
& (misdat.eq.ar(ip1,jp1,kp1))) then
|
|
|
754 |
int3dmvc=misdat
|
|
|
755 |
else
|
|
|
756 |
frac0i=ri-float(i)
|
|
|
757 |
frac0j=rj-float(j)
|
|
|
758 |
frac1i=1.-frac0i
|
|
|
759 |
frac1j=1.-frac0j
|
|
|
760 |
if (k-1.lt.1) then
|
|
|
761 |
klow=1
|
|
|
762 |
else if (k+2.gt.n3) then
|
|
|
763 |
klow=n3-3
|
|
|
764 |
else
|
|
|
765 |
klow=k-1
|
|
|
766 |
endif
|
|
|
767 |
do n=1,4
|
|
|
768 |
int2d(n) = ar(i ,j ,klow-1+n ) * frac1i * frac1j
|
|
|
769 |
& + ar(i ,jp1,klow-1+n ) * frac1i * frac0j
|
|
|
770 |
& + ar(ip1,j ,klow-1+n ) * frac0i * frac1j
|
|
|
771 |
& + ar(ip1,jp1,klow-1+n ) * frac0i * frac0j
|
|
|
772 |
enddo
|
|
|
773 |
call cubint(int2d(1),int2d(2),int2d(3),int2d(4),
|
|
|
774 |
& klow,rk,int3dmvc)
|
|
|
775 |
endif
|
|
|
776 |
endif
|
|
|
777 |
endif
|
|
|
778 |
end
|
|
|
779 |
real function int4d(ar,n1,n2,n3,n4,rid,rjd,rkd,rld)
|
|
|
780 |
c-----------------------------------------------------------------------
|
|
|
781 |
c Purpose:
|
|
|
782 |
c This subroutine interpolates a 4d-array to an arbitrary
|
|
|
783 |
c location within the grid.
|
|
|
784 |
c Arguments:
|
|
|
785 |
c ar real input surface pressure, define as ar(n1,n2,n3)
|
|
|
786 |
c n1,..,n4 int input dimensions of ar
|
|
|
787 |
c ri,..,rl real input grid location to be interpolated to
|
|
|
788 |
c History:
|
|
|
789 |
c-----------------------------------------------------------------------
|
|
|
790 |
|
|
|
791 |
c argument declarations
|
|
|
792 |
integer n1,n2,n3,n4
|
|
|
793 |
real ar(n1,n2,n3,n4), rid,rjd,rkd,rld
|
|
|
794 |
|
|
|
795 |
c local declarations
|
|
|
796 |
integer l,lp1,lh
|
|
|
797 |
real frac0l,frac1l,rl,int3d
|
|
|
798 |
|
|
|
799 |
c do linear interpolation in l-direction
|
|
|
800 |
rl=amax1(1.,amin1(float(n4),rld))
|
|
|
801 |
lh=nint(rl)
|
|
|
802 |
|
|
|
803 |
c Check for interpolation in l
|
|
|
804 |
* if (abs(float(lh)-rl).lt.1.e-3) then
|
|
|
805 |
* l =lh
|
|
|
806 |
* lp1=lh
|
|
|
807 |
* else
|
|
|
808 |
l =min0(int(rl),n4-1)
|
|
|
809 |
lp1=l+1
|
|
|
810 |
* endif
|
|
|
811 |
|
|
|
812 |
if (l.eq.lp1) then
|
|
|
813 |
c no interpolation in l
|
|
|
814 |
int4d=int3d(ar(1,1,1,l),n1,n2,n3,rid,rjd,rkd)
|
|
|
815 |
else
|
|
|
816 |
c interpolation in l
|
|
|
817 |
frac0l=rl-float(l)
|
|
|
818 |
frac1l=1.-frac0l
|
|
|
819 |
int4d = int3d(ar(1,1,1,l ),n1,n2,n3,rid,rjd,rkd) * frac1l
|
|
|
820 |
& + int3d(ar(1,1,1,lp1),n1,n2,n3,rid,rjd,rkd) * frac0l
|
|
|
821 |
endif
|
|
|
822 |
end
|
|
|
823 |
real function int4dm(ar,n1,n2,n3,n4,rid,rjd,rkd,rld,misdat)
|
|
|
824 |
c-----------------------------------------------------------------------
|
|
|
825 |
c Purpose:
|
|
|
826 |
c This subroutine interpolates a 4d-array to an arbitrary
|
|
|
827 |
c location within the grid. The interpolation includes the
|
|
|
828 |
c testing of the missing data flag 'misdat'.
|
|
|
829 |
c Arguments:
|
|
|
830 |
c ar real input surface pressure, define as ar(n1,n2,n3)
|
|
|
831 |
c n1,..,n4 int input dimensions of ar
|
|
|
832 |
c ri,..,rl real input grid location to be interpolated to
|
|
|
833 |
c misdat real input missing data flag (on if misdat<>0)
|
|
|
834 |
c Warning:
|
|
|
835 |
c This routine has not yet been seriously tested.
|
|
|
836 |
c History:
|
|
|
837 |
c-----------------------------------------------------------------------
|
|
|
838 |
|
|
|
839 |
c argument declarations
|
|
|
840 |
integer n1,n2,n3,n4
|
|
|
841 |
real ar(n1,n2,n3,n4), rid,rjd,rkd,rld, misdat
|
|
|
842 |
|
|
|
843 |
c local declarations
|
|
|
844 |
integer l,lp1,lh
|
|
|
845 |
real frac0l,frac1l,rl,rint0,rint1,int4d,int3dm
|
|
|
846 |
|
|
|
847 |
c check whether missing data checking is required
|
|
|
848 |
if (misdat.eq.0.) then
|
|
|
849 |
int4dm=int4d(ar,n1,n2,n3,n4,rid,rjd,rkd,rld)
|
|
|
850 |
return
|
|
|
851 |
endif
|
|
|
852 |
|
|
|
853 |
c do linear interpolation in l-direction
|
|
|
854 |
rl=amax1(1.,amin1(float(n4),rld))
|
|
|
855 |
lh=nint(rl)
|
|
|
856 |
|
|
|
857 |
c Check for interpolation in l
|
|
|
858 |
* if (abs(float(lh)-rl).lt.1.e-3) then
|
|
|
859 |
* l =lh
|
|
|
860 |
* lp1=lh
|
|
|
861 |
* else
|
|
|
862 |
l =min0(int(rl),n4-1)
|
|
|
863 |
lp1=l+1
|
|
|
864 |
* endif
|
|
|
865 |
|
|
|
866 |
if (l.eq.lp1) then
|
|
|
867 |
c no interpolation in l
|
|
|
868 |
int4dm = int3dm(ar(1,1,1,l),n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
869 |
else
|
|
|
870 |
c interpolation in l
|
|
|
871 |
frac0l=rl-float(l)
|
|
|
872 |
frac1l=1.-frac0l
|
|
|
873 |
rint0 = int3dm(ar(1,1,1,l ),n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
874 |
rint1 = int3dm(ar(1,1,1,lp1),n1,n2,n3,rid,rjd,rkd,misdat)
|
|
|
875 |
if ((rint0.eq.misdat).or.(rint1.eq.misdat)) then
|
|
|
876 |
int4dm = misdat
|
|
|
877 |
else
|
|
|
878 |
int4dm = rint0*frac1l + rint1*frac0l
|
|
|
879 |
endif
|
|
|
880 |
endif
|
|
|
881 |
end
|
|
|
882 |
real function int3dl(ar,n1,n2,n3,levels,rid,rjd,rkd)
|
|
|
883 |
c-----------------------------------------------------------------------
|
|
|
884 |
c Purpose:
|
|
|
885 |
c This subroutine interpolates a 3d-array to an arbitrary
|
|
|
886 |
c location within the grid. The vertical interpolation is linear
|
|
|
887 |
c in log(pressure).
|
|
|
888 |
c Arguments:
|
|
|
889 |
c ar real input surface pressure, define as ar(n1,n2,n3)
|
|
|
890 |
c n1,n2,n3 int input dimensions of ar
|
|
|
891 |
c levels real input array contains pressure levels for ar
|
|
|
892 |
c ri,rj,rk real input grid location to be interpolated to
|
|
|
893 |
c History:
|
|
|
894 |
c Based on int3d July 93
|
|
|
895 |
c-----------------------------------------------------------------------
|
|
|
896 |
|
|
|
897 |
c argument declarations
|
|
|
898 |
integer n1,n2,n3
|
|
|
899 |
real ar(n1,n2,n3), rid,rjd,rkd
|
|
|
900 |
real levels(n3)
|
|
|
901 |
|
|
|
902 |
c local declarations
|
|
|
903 |
real pval
|
|
|
904 |
integer i,j,k,ip1,jp1,kp1,ih,jh,kh
|
|
|
905 |
real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk
|
|
|
906 |
|
|
|
907 |
c do linear interpolation
|
|
|
908 |
ri=amax1(1.,amin1(float(n1),rid))
|
|
|
909 |
rj=amax1(1.,amin1(float(n2),rjd))
|
|
|
910 |
rk=amax1(1.,amin1(float(n3),rkd))
|
|
|
911 |
ih=nint(ri)
|
|
|
912 |
jh=nint(rj)
|
|
|
913 |
kh=nint(rk)
|
|
|
914 |
|
|
|
915 |
c Check for interpolation in i
|
|
|
916 |
* if (abs(float(ih)-ri).lt.1.e-3) then
|
|
|
917 |
* i =ih
|
|
|
918 |
* ip1=ih
|
|
|
919 |
* else
|
|
|
920 |
i =min0(int(ri),n1-1)
|
|
|
921 |
ip1=i+1
|
|
|
922 |
* endif
|
|
|
923 |
|
|
|
924 |
c Check for interpolation in j
|
|
|
925 |
* if (abs(float(jh)-rj).lt.1.e-3) then
|
|
|
926 |
* j =jh
|
|
|
927 |
* jp1=jh
|
|
|
928 |
* else
|
|
|
929 |
j =min0(int(rj),n2-1)
|
|
|
930 |
jp1=j+1
|
|
|
931 |
* endif
|
|
|
932 |
|
|
|
933 |
c Check for interpolation in k
|
|
|
934 |
* if (abs(float(kh)-rk).lt.1.e-3) then
|
|
|
935 |
* k =kh
|
|
|
936 |
* kp1=kh
|
|
|
937 |
* else
|
|
|
938 |
k =min0(int(rk),n3-1)
|
|
|
939 |
kp1=k+1
|
|
|
940 |
* endif
|
|
|
941 |
|
|
|
942 |
if (k.eq.kp1) then
|
|
|
943 |
c no interpolation in k
|
|
|
944 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
945 |
c no interpolation at all
|
|
|
946 |
int3dl=ar(i,j,k)
|
|
|
947 |
c print *,'int3dl 00: ',rid,rjd,rkd,int3dl
|
|
|
948 |
else
|
|
|
949 |
c horizontal interpolation only
|
|
|
950 |
frac0i=ri-float(i)
|
|
|
951 |
frac0j=rj-float(j)
|
|
|
952 |
frac1i=1.-frac0i
|
|
|
953 |
frac1j=1.-frac0j
|
|
|
954 |
int3dl = ar(i ,j ,k ) * frac1i * frac1j
|
|
|
955 |
& + ar(i ,jp1,k ) * frac1i * frac0j
|
|
|
956 |
& + ar(ip1,j ,k ) * frac0i * frac1j
|
|
|
957 |
& + ar(ip1,jp1,k ) * frac0i * frac0j
|
|
|
958 |
c print *,'int3dl 10: ',rid,rjd,rkd,int3dl
|
|
|
959 |
endif
|
|
|
960 |
else
|
|
|
961 |
* frac0k=rk-float(k)
|
|
|
962 |
c calculate the pressure value to be interpolated to
|
|
|
963 |
pval=levels(int(rk))
|
|
|
964 |
> -(rk-aint(rk))*(levels(int(rk))-levels(int(rk)+1))
|
|
|
965 |
frac0k=log(levels(int(rk))/pval)
|
|
|
966 |
& /log(levels(int(rk))/levels(int(rk)+1))
|
|
|
967 |
frac1k=1.-frac0k
|
|
|
968 |
if ((i.eq.ip1).and.(j.eq.jp1)) then
|
|
|
969 |
c vertical interpolation only
|
|
|
970 |
int3dl = ar(i ,j ,k ) * frac1k
|
|
|
971 |
& + ar(i ,j ,kp1) * frac0k
|
|
|
972 |
c print *,'int3dl 01: ',rid,rjd,rkd,int3dl
|
|
|
973 |
else
|
|
|
974 |
c full 3d interpolation
|
|
|
975 |
frac0i=ri-float(i)
|
|
|
976 |
frac0j=rj-float(j)
|
|
|
977 |
frac1i=1.-frac0i
|
|
|
978 |
frac1j=1.-frac0j
|
|
|
979 |
int3dl = ar(i ,j ,k ) * frac1i * frac1j * frac1k
|
|
|
980 |
& + ar(i ,jp1,k ) * frac1i * frac0j * frac1k
|
|
|
981 |
& + ar(ip1,j ,k ) * frac0i * frac1j * frac1k
|
|
|
982 |
& + ar(ip1,jp1,k ) * frac0i * frac0j * frac1k
|
|
|
983 |
& + ar(i ,j ,kp1) * frac1i * frac1j * frac0k
|
|
|
984 |
& + ar(i ,jp1,kp1) * frac1i * frac0j * frac0k
|
|
|
985 |
& + ar(ip1,j ,kp1) * frac0i * frac1j * frac0k
|
|
|
986 |
& + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
|
|
|
987 |
c print *,'int3dl 11: ',rid,rjd,rkd,int3dl
|
|
|
988 |
endif
|
|
|
989 |
endif
|
|
|
990 |
end
|
|
|
991 |
subroutine bcucof(y,y1,y2,y12,c)
|
|
|
992 |
c-----------------------------------------------------------------------
|
|
|
993 |
c Given arrays y,y1,y2 and y12, each of length 4, containing the
|
|
|
994 |
c function, gradients, and cross derivative at the four grid points
|
|
|
995 |
c of a rectangular grid cell (numbered counterclockwise from the
|
|
|
996 |
c lower left), and given d1 and d2, the length of the grid cell in
|
|
|
997 |
c the 1- and 2-directions, this routine returns the table c that is
|
|
|
998 |
c used by routine bcuint for biqubic interpolation.
|
|
|
999 |
c Source: Numerical Recipes, Fortran Version, p.99
|
|
|
1000 |
c-----------------------------------------------------------------------
|
|
|
1001 |
real c(4,4),y(4),y1(4),y2(4),y12(4),cl(16),x(16)
|
|
|
1002 |
integer wt(16,16)
|
|
|
1003 |
|
|
|
1004 |
data wt/1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,
|
|
|
1005 |
> 8*0,3,0,-9,6,-2,0,6,-4,
|
|
|
1006 |
> 10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,2*0,6,-4,
|
|
|
1007 |
> 4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,2,
|
|
|
1008 |
> 10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2,
|
|
|
1009 |
> 0,1,-2,1,5*0,-3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,
|
|
|
1010 |
> 10*0,-3,3,2*0,2,-2,2*0,-1,1,6*0,3,-3,2*0,-2,2,
|
|
|
1011 |
> 5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2,-1,0,1,-2,1,
|
|
|
1012 |
> 10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/
|
|
|
1013 |
|
|
|
1014 |
real xx
|
|
|
1015 |
integer i,j,k,l
|
|
|
1016 |
|
|
|
1017 |
do i=1,4 ! pack a temporary vector x
|
|
|
1018 |
x(i)=y(i)
|
|
|
1019 |
x(i+4)=y1(i)
|
|
|
1020 |
x(i+8)=y2(i)
|
|
|
1021 |
x(i+12)=y12(i)
|
|
|
1022 |
enddo
|
|
|
1023 |
|
|
|
1024 |
do i=1,16 ! matrix multiply by the stored table
|
|
|
1025 |
xx=0.
|
|
|
1026 |
do k=1,16
|
|
|
1027 |
xx=xx+wt(i,k)*x(k)
|
|
|
1028 |
enddo
|
|
|
1029 |
cl(i)=xx
|
|
|
1030 |
enddo
|
|
|
1031 |
|
|
|
1032 |
l=0
|
|
|
1033 |
do i=1,4 ! unpack the result into the output table
|
|
|
1034 |
do j=1,4
|
|
|
1035 |
l=l+1
|
|
|
1036 |
c(i,j)=cl(l)
|
|
|
1037 |
enddo
|
|
|
1038 |
enddo
|
|
|
1039 |
|
|
|
1040 |
return
|
|
|
1041 |
end
|
|
|
1042 |
subroutine bcuint(y,y1,y2,y12,x1l,x2l,x1,x2,ansy)
|
|
|
1043 |
c-----------------------------------------------------------------------
|
|
|
1044 |
c Bicubic interpolation within a grid square. Input quantities are
|
|
|
1045 |
c y,y1,y2,y12 (as described in bcucof); x1l and x1u, the lower and
|
|
|
1046 |
c upper coordinates of the grid square in the 1-direction; x2l and
|
|
|
1047 |
c x2u likewise for the 2-direction; and x1,x2, the coordinates of
|
|
|
1048 |
c the desired point for the interpolation. The interplated function
|
|
|
1049 |
c value is returned as ansy. This routine calls bcucof.
|
|
|
1050 |
c Source: Numerical Recipes, Fortran Version, p.99/100
|
|
|
1051 |
c !!! changed the proposed code !!!
|
|
|
1052 |
c-----------------------------------------------------------------------
|
|
|
1053 |
|
|
|
1054 |
real y(4),y1(4),y2(4),y12(4),c(4,4)
|
|
|
1055 |
real ansy,x1,x2,t,u
|
|
|
1056 |
integer i,x1l,x2l
|
|
|
1057 |
|
|
|
1058 |
call bcucof(y,y1,y2,y12,c)
|
|
|
1059 |
|
|
|
1060 |
t=x1-real(x1l)
|
|
|
1061 |
u=x2-real(x2l)
|
|
|
1062 |
|
|
|
1063 |
ansy=0.
|
|
|
1064 |
|
|
|
1065 |
do i=4,1,-1
|
|
|
1066 |
ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
|
|
|
1067 |
enddo
|
|
|
1068 |
|
|
|
1069 |
return
|
|
|
1070 |
end
|
|
|
1071 |
subroutine cubint(ya1,ya2,ya3,ya4,k,x,y)
|
|
|
1072 |
c-----------------------------------------------------------------------
|
|
|
1073 |
c Interface routine for SR polint for special case of cubic
|
|
|
1074 |
c interpolation in index space, with xa=k,k+1,k+2,k+3
|
|
|
1075 |
c-----------------------------------------------------------------------
|
|
|
1076 |
|
|
|
1077 |
integer k
|
|
|
1078 |
real ya1,ya2,ya3,ya4,x,y
|
|
|
1079 |
|
|
|
1080 |
integer n
|
|
|
1081 |
real xa(4),ya(4),dy
|
|
|
1082 |
|
|
|
1083 |
do n=1,4
|
|
|
1084 |
xa(1)=real(k)
|
|
|
1085 |
xa(2)=real(k+1)
|
|
|
1086 |
xa(3)=real(k+2)
|
|
|
1087 |
xa(4)=real(k+3)
|
|
|
1088 |
|
|
|
1089 |
ya(1)=ya1
|
|
|
1090 |
ya(2)=ya2
|
|
|
1091 |
ya(3)=ya3
|
|
|
1092 |
ya(4)=ya4
|
|
|
1093 |
enddo
|
|
|
1094 |
|
|
|
1095 |
call polint(xa,ya,4,x,y,dy)
|
|
|
1096 |
|
|
|
1097 |
return
|
|
|
1098 |
end
|
|
|
1099 |
subroutine polint(xa,ya,n,x,y,dy)
|
|
|
1100 |
c-----------------------------------------------------------------------
|
|
|
1101 |
c Given arrays xa and ya, each of length n, and given a value x, this
|
|
|
1102 |
c routine returns a value y, and an error estimate dy. If P(x) is the
|
|
|
1103 |
c polynomial of degree n-1 such that p(xa(i))=ya(i),i=1,...,n, then
|
|
|
1104 |
c the returned value y=p(x)
|
|
|
1105 |
c Source: Numerical Recipes, Fortran Version, p.82
|
|
|
1106 |
c-----------------------------------------------------------------------
|
|
|
1107 |
|
|
|
1108 |
integer nmax,n
|
|
|
1109 |
parameter(nmax=10)
|
|
|
1110 |
real xa(n),ya(n),x,y,dy
|
|
|
1111 |
integer i,m,ns
|
|
|
1112 |
real c(nmax),d(nmax),den,dif,dift,ho,hp,w
|
|
|
1113 |
|
|
|
1114 |
ns=1
|
|
|
1115 |
dif=abs(x-xa(1))
|
|
|
1116 |
|
|
|
1117 |
do i=1,n
|
|
|
1118 |
dift=abs(x-xa(i))
|
|
|
1119 |
if (dift.lt.dif) then
|
|
|
1120 |
ns=i
|
|
|
1121 |
dif=dift
|
|
|
1122 |
endif
|
|
|
1123 |
c(i)=ya(i)
|
|
|
1124 |
d(i)=ya(i)
|
|
|
1125 |
enddo
|
|
|
1126 |
|
|
|
1127 |
y=ya(ns)
|
|
|
1128 |
ns=ns-1
|
|
|
1129 |
do m=1,n-1
|
|
|
1130 |
do i=1,n-m
|
|
|
1131 |
ho=xa(i)-x
|
|
|
1132 |
hp=xa(i+m)-x
|
|
|
1133 |
w=c(i+1)-d(i)
|
|
|
1134 |
den=ho-hp
|
|
|
1135 |
den=w/den
|
|
|
1136 |
d(i)=hp*den
|
|
|
1137 |
c(i)=ho*den
|
|
|
1138 |
enddo
|
|
|
1139 |
if (2*ns.lt.n-m) then
|
|
|
1140 |
dy=c(ns+1)
|
|
|
1141 |
else
|
|
|
1142 |
dy=d(ns)
|
|
|
1143 |
ns=ns-1
|
|
|
1144 |
endif
|
|
|
1145 |
y=y+dy
|
|
|
1146 |
enddo
|
|
|
1147 |
|
|
|
1148 |
return
|
|
|
1149 |
end
|
|
|
1150 |
subroutine filt2d (a,af,f1,f2,nx,ny,fil,misdat,
|
|
|
1151 |
& iperx,ipery,ispol,inpol)
|
|
|
1152 |
c =============================================================
|
|
|
1153 |
c Apply a conservative diffusion operator onto the 2d field a,
|
|
|
1154 |
c with full missing data checking.
|
|
|
1155 |
c
|
|
|
1156 |
c a real inp array to be filtered, dimensioned (nx,ny)
|
|
|
1157 |
c af real out filtered array, dimensioned (nx,ny), can be
|
|
|
1158 |
c equivalenced with array a in the calling routine
|
|
|
1159 |
c f1 real workarray, dimensioned (nx+1,ny)
|
|
|
1160 |
c f2 real workarray, dimensioned (nx,ny+1)
|
|
|
1161 |
c fil real inp filter-coeff., 0<afil<=1. Maximum filtering with afil=1
|
|
|
1162 |
c corresponds to one application of the linear filter.
|
|
|
1163 |
c misdat real inp missing-data value, a(i,j)=misdat indicates that
|
|
|
1164 |
c the corresponding value is not available. The
|
|
|
1165 |
c misdat-checking can be switched off with with misdat=0.
|
|
|
1166 |
c iperx int inp periodic boundaries in the x-direction (1=yes,0=no)
|
|
|
1167 |
c ipery int inp periodic boundaries in the y-direction (1=yes,0=no)
|
|
|
1168 |
c inpol int inp northpole at j=ny (1=yes,0=no)
|
|
|
1169 |
c ispol int inp southpole at j=1 (1=yes,0=no)
|
|
|
1170 |
c
|
|
|
1171 |
c Christoph Schaer, 1993
|
|
|
1172 |
|
|
|
1173 |
c argument declaration
|
|
|
1174 |
integer nx,ny
|
|
|
1175 |
real a(nx,ny),af(nx,ny),f1(nx+1,ny),f2(nx,ny+1),fil,misdat
|
|
|
1176 |
integer iperx,ipery,inpol,ispol
|
|
|
1177 |
|
|
|
1178 |
c local variable declaration
|
|
|
1179 |
integer i,j,is
|
|
|
1180 |
real fh
|
|
|
1181 |
|
|
|
1182 |
c compute constant fh
|
|
|
1183 |
fh=0.125*fil
|
|
|
1184 |
|
|
|
1185 |
c compute fluxes in x-direction
|
|
|
1186 |
if (misdat.eq.0.) then
|
|
|
1187 |
do j=1,ny
|
|
|
1188 |
do i=2,nx
|
|
|
1189 |
f1(i,j)=a(i-1,j)-a(i,j)
|
|
|
1190 |
enddo
|
|
|
1191 |
enddo
|
|
|
1192 |
else
|
|
|
1193 |
do j=1,ny
|
|
|
1194 |
do i=2,nx
|
|
|
1195 |
if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
|
|
|
1196 |
f1(i,j)=0.
|
|
|
1197 |
else
|
|
|
1198 |
f1(i,j)=a(i-1,j)-a(i,j)
|
|
|
1199 |
endif
|
|
|
1200 |
enddo
|
|
|
1201 |
enddo
|
|
|
1202 |
endif
|
|
|
1203 |
if (iperx.eq.1) then
|
|
|
1204 |
c do periodic boundaries in the x-direction
|
|
|
1205 |
do j=1,ny
|
|
|
1206 |
f1(1,j)=f1(nx,j)
|
|
|
1207 |
f1(nx+1,j)=f1(2,j)
|
|
|
1208 |
enddo
|
|
|
1209 |
else
|
|
|
1210 |
c set boundary-fluxes to zero
|
|
|
1211 |
do j=1,ny
|
|
|
1212 |
f1(1,j)=0.
|
|
|
1213 |
f1(nx+1,j)=0.
|
|
|
1214 |
enddo
|
|
|
1215 |
endif
|
|
|
1216 |
|
|
|
1217 |
c compute fluxes in y-direction
|
|
|
1218 |
if (misdat.eq.0.) then
|
|
|
1219 |
do j=2,ny
|
|
|
1220 |
do i=1,nx
|
|
|
1221 |
f2(i,j)=a(i,j-1)-a(i,j)
|
|
|
1222 |
enddo
|
|
|
1223 |
enddo
|
|
|
1224 |
else
|
|
|
1225 |
do j=2,ny
|
|
|
1226 |
do i=1,nx
|
|
|
1227 |
if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
|
|
|
1228 |
f2(i,j)=0.
|
|
|
1229 |
else
|
|
|
1230 |
f2(i,j)=a(i,j-1)-a(i,j)
|
|
|
1231 |
endif
|
|
|
1232 |
enddo
|
|
|
1233 |
enddo
|
|
|
1234 |
endif
|
|
|
1235 |
c set boundary-fluxes to zero
|
|
|
1236 |
do i=1,nx
|
|
|
1237 |
f2(i,1)=0.
|
|
|
1238 |
f2(i,ny+1)=0.
|
|
|
1239 |
enddo
|
|
|
1240 |
if (ipery.eq.1) then
|
|
|
1241 |
c do periodic boundaries in the x-direction
|
|
|
1242 |
do i=1,nx
|
|
|
1243 |
f2(i,1)=f2(i,ny)
|
|
|
1244 |
f2(i,ny+1)=f2(i,2)
|
|
|
1245 |
enddo
|
|
|
1246 |
endif
|
|
|
1247 |
if (iperx.eq.1) then
|
|
|
1248 |
if (ispol.eq.1) then
|
|
|
1249 |
c do south-pole
|
|
|
1250 |
is=(nx-1)/2
|
|
|
1251 |
do i=1,nx
|
|
|
1252 |
f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
|
|
|
1253 |
enddo
|
|
|
1254 |
endif
|
|
|
1255 |
if (inpol.eq.1) then
|
|
|
1256 |
c do north-pole
|
|
|
1257 |
is=(nx-1)/2
|
|
|
1258 |
do i=1,nx
|
|
|
1259 |
f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
|
|
|
1260 |
enddo
|
|
|
1261 |
endif
|
|
|
1262 |
endif
|
|
|
1263 |
|
|
|
1264 |
c compute flux-convergence -> filter
|
|
|
1265 |
if (misdat.eq.0.) then
|
|
|
1266 |
do j=1,ny
|
|
|
1267 |
do i=1,nx
|
|
|
1268 |
af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
|
|
|
1269 |
enddo
|
|
|
1270 |
enddo
|
|
|
1271 |
else
|
|
|
1272 |
do j=1,ny
|
|
|
1273 |
do i=1,nx
|
|
|
1274 |
if (a(i,j).eq.misdat) then
|
|
|
1275 |
af(i,j)=misdat
|
|
|
1276 |
else
|
|
|
1277 |
af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
|
|
|
1278 |
endif
|
|
|
1279 |
enddo
|
|
|
1280 |
enddo
|
|
|
1281 |
endif
|
|
|
1282 |
end
|
|
|
1283 |
|
|
|
1284 |
subroutine pipo(var3d,p3d,lev,var,nx,ny,nz,mdv,mode)
|
|
|
1285 |
C ====================================================
|
|
|
1286 |
|
|
|
1287 |
C Interpolates the 3d variable var3d on the pressure surface
|
|
|
1288 |
C defined by lev. The interpolated field is returned as var.
|
|
|
1289 |
C p3d denotes the 3d pressure array.
|
|
|
1290 |
C mode determines the way of vertical interpolation:
|
|
|
1291 |
C mode=0 is for linear interpolation
|
|
|
1292 |
C mode=1 is for logarithmic interpolation
|
|
|
1293 |
|
|
|
1294 |
integer nx,ny,nz,mode
|
|
|
1295 |
real lev,mdv
|
|
|
1296 |
real var3d(nx,ny,nz),p3d(nx,ny,nz),var(nx,ny)
|
|
|
1297 |
|
|
|
1298 |
integer i,j,k
|
|
|
1299 |
real kind
|
|
|
1300 |
real int3dm
|
|
|
1301 |
|
|
|
1302 |
do i=1,nx
|
|
|
1303 |
do j=1,ny
|
|
|
1304 |
|
|
|
1305 |
kind=0.
|
|
|
1306 |
do k=1,nz-1
|
|
|
1307 |
if ((p3d(i,j,k).ge.lev).and.(p3d(i,j,k+1).le.lev)) then
|
|
|
1308 |
kind=float(k)+(p3d(i,j,k)-lev)/
|
|
|
1309 |
> (p3d(i,j,k)-p3d(i,j,k+1))
|
|
|
1310 |
goto 100
|
|
|
1311 |
endif
|
|
|
1312 |
enddo
|
|
|
1313 |
100 continue
|
|
|
1314 |
|
|
|
1315 |
if (kind.eq.0.) then
|
|
|
1316 |
var(i,j)=mdv
|
|
|
1317 |
else
|
|
|
1318 |
var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
|
|
|
1319 |
endif
|
|
|
1320 |
|
|
|
1321 |
enddo
|
|
|
1322 |
enddo
|
|
|
1323 |
|
|
|
1324 |
return
|
|
|
1325 |
end
|
|
|
1326 |
|
|
|
1327 |
subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
|
|
|
1328 |
C ======================================================
|
|
|
1329 |
|
|
|
1330 |
C Interpolates the 3d variable var3d on the isentropic surface
|
|
|
1331 |
C defined by lev. The interpolated field is returned as var.
|
|
|
1332 |
C th3d denotes the 3d theta array.
|
|
|
1333 |
C mode determines the way of vertical interpolation:
|
|
|
1334 |
C mode=0 is for linear interpolation
|
|
|
1335 |
C mode=1 is for logarithmic interpolation
|
|
|
1336 |
|
|
|
1337 |
integer nx,ny,nz,mode
|
|
|
1338 |
real lev,mdv
|
|
|
1339 |
real var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)
|
|
|
1340 |
|
|
|
1341 |
integer i,j,k
|
|
|
1342 |
real kind
|
|
|
1343 |
real int3dm
|
|
|
1344 |
|
|
|
1345 |
do i=1,nx
|
|
|
1346 |
do j=1,ny
|
|
|
1347 |
|
|
|
1348 |
kind=0
|
|
|
1349 |
do k=1,nz-1
|
|
|
1350 |
if ((th3d(i,j,k).le.lev).and.(th3d(i,j,k+1).ge.lev)) then
|
|
|
1351 |
kind=float(k)+(th3d(i,j,k)-lev)/
|
|
|
1352 |
> (th3d(i,j,k)-th3d(i,j,k+1))
|
|
|
1353 |
goto 100
|
|
|
1354 |
endif
|
|
|
1355 |
enddo
|
|
|
1356 |
100 continue
|
|
|
1357 |
|
|
|
1358 |
if (kind.eq.0) then
|
|
|
1359 |
var(i,j)=mdv
|
|
|
1360 |
else
|
|
|
1361 |
var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
|
|
|
1362 |
endif
|
|
|
1363 |
|
|
|
1364 |
enddo
|
|
|
1365 |
enddo
|
|
|
1366 |
|
|
|
1367 |
return
|
|
|
1368 |
end
|