Subversion Repositories tropofold.echam

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 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