File: C:\NOAA\NEMS_11731\src\chem\gocart\src\GMAO_Shared\MAPL_Base\hinterp.F
1 subroutine hinterp ( qin,iin,jin,qout,iout,jout,mlev,undef )
2 implicit none
3 integer iin,jin, iout,jout, mlev
4 real qin(iin,jin,mlev), qout(iout,jout,mlev)
5 real undef,pi,dlin,dpin,dlout,dpout
6 real dlam(iin), lons(iout*jout), lon
7 real dphi(jin), lats(iout*jout), lat
8 integer i,j,loc
9
10 pi = 4.0*atan(1.0)
11 dlin = 2*pi/iin
12 dpin = pi/(jin-1)
13 dlam(:) = dlin
14 dphi(:) = dpin
15
16 dlout = 2*pi/iout
17 dpout = pi/(jout-1)
18
19 loc = 0
20 do j=1,jout
21 do i=1,iout
22 loc = loc + 1
23 lon = -pi + (i-1)*dlout
24 lons(loc) = lon
25 enddo
26 enddo
27
28 loc = 0
29 do j=1,jout
30 lat = -pi/2.0 + (j-1)*dpout
31 do i=1,iout
32 loc = loc + 1
33 lats(loc) = lat
34 enddo
35 enddo
36
37 call interp_hh ( qin,iin,jin,mlev,dlam,dphi,
38 . qout,iout*jout,lons,lats,undef, -pi )
39
40 return
41 end
42
43
44
45
46 subroutine hhinterp ( qin,iin,jin,qout,iout,jout,mlev,undef,
47 . lons_in, lats_in )
48 implicit none
49 integer iin,jin, iout,jout, mlev
50 real qin(iin,jin,mlev), qout(iout,jout,mlev)
51 real undef,pi,dlin,dpin,dlout,dpout
52 real :: lons_in(iin), lats_in(jin)
53 real dlam(iin), lons(iout*jout), lon
54 real dphi(jin), lats(iout*jout), lat
55 integer i,j,loc,rc
56 real lon_min
57
58 character(len=*), parameter :: Iam='hinterp_'
59 integer :: i1, i2, i3, i4, n1, n2, n3, n4, imh, k
60 real lons_save(iin)
61
62 pi = 4.0*atan(1.0)
63
64
65
66
67
68
69 = lons_in
70 imh = -1
71 if ( lons_in(1) >= 0 ) then
72 do i = 1, iin
73 if ( lons_in(i) >= pi ) then
74 lons_in(i) = lons_in(i) - 2*pi
75 if ( imh < 0 ) imh = i
76 end if
77 end do
78 i1 = 1; i2 = imh-1; i3=imh; i4=iin
79 n1 = 1; n2 = i4-i3+1; n3=n2+1; n4=iin
80 dlam = lons_in
81 lons_in(n1:n2) = dlam(i3:i4)
82 lons_in(n3:n4) = dlam(i1:i2)
83 do k = 1, mlev
84 do j = 1, jin
85 dlam = qin(:,j,k)
86 qin(n1:n2,j,k) = dlam(i3:i4)
87 qin(n3:n4,j,k) = dlam(i1:i2)
88 end do
89 end do
90 end if
91
92 do i = 1, iin-1
93 dlam(i) = lons_in(i+1) - lons_in(i)
94 enddo
95 dlam(iin) = lons_in(iin) - lons_in(iin-1)
96 = lons_in(1)
97
98 if ( abs(lats_in(1)+pi/2) >= 0.0001 ) then
99 qout = undef
100 return
101 end if
102 do j = 1, jin-1
103 dphi(j) = lats_in(j+1) - lats_in(j)
104 enddo
105 dphi(jin) = lats_in(jin) - lats_in(jin-1)
106
107 = 2*pi/iout
108 dpout = pi/(jout-1)
109
110 loc = 0
111 do j=1,jout
112 do i=1,iout
113 loc = loc + 1
114 lon = -pi + (i-1)*dlout
115 (loc) = lon
116 enddo
117 enddo
118
119 loc = 0
120 do j=1,jout
121 lat = -pi/2.0 + (j-1)*dpout
122 do i=1,iout
123 loc = loc + 1
124 lats(loc) = lat
125 enddo
126 enddo
127
128 call interp_hh ( qin,iin,jin,mlev,dlam,dphi,
129 . qout,iout*jout,lons,lats,undef, lon_min )
130
131
132
133 if ( imh > 0 ) then
134 lons_in = lons_save
135 do k = 1, mlev
136 do j = 1, jin
137 dlam = qin(:,j,k)
138 qin(i1:i2,j,k) = dlam(n3:n4)
139 qin(i3:i4,j,k) = dlam(n1:n2)
140 end do
141 end do
142 end if
143
144 return
145 end
146
147 subroutine interp_hh ( q_cmp,im,jm,lm,dlam,dphi,
148 . q_geo,irun,lon_geo,lat_geo, undef, lon_min )
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177 implicit none
178
179
180
181 integer im,jm,lm,irun
182
183 real q_geo(irun,lm)
184 real lon_geo(irun)
185 real lat_geo(irun)
186
187 real q_cmp(im,jm,lm)
188 real dlam(im)
189 real dphi(jm)
190
191 real :: lon_min
192
193
194
195 integer i,j,l,m,n
196 integer, allocatable :: ip1(:), ip0(:), im1(:), im2(:)
197 integer, allocatable :: jp1(:), jp0(:), jm1(:), jm2(:)
198
199 integer ip1_for_jp1, ip0_for_jp1, im1_for_jp1, im2_for_jp1
200 integer ip1_for_jm2, ip0_for_jm2, im1_for_jm2, im2_for_jm2
201 integer jm2_for_jm2, jp1_for_jp1
202
203
204
205 real, allocatable :: wl_ip0jp0 (:)
206 real, allocatable :: wl_im1jp0 (:)
207 real, allocatable :: wl_ip0jm1 (:)
208 real, allocatable :: wl_im1jm1 (:)
209
210
211
212 real, allocatable :: wc_ip1jp1 (:)
213 real, allocatable :: wc_ip0jp1 (:)
214 real, allocatable :: wc_im1jp1 (:)
215 real, allocatable :: wc_im2jp1 (:)
216 real, allocatable :: wc_ip1jp0 (:)
217 real, allocatable :: wc_ip0jp0 (:)
218 real, allocatable :: wc_im1jp0 (:)
219 real, allocatable :: wc_im2jp0 (:)
220 real, allocatable :: wc_ip1jm1 (:)
221 real, allocatable :: wc_ip0jm1 (:)
222 real, allocatable :: wc_im1jm1 (:)
223 real, allocatable :: wc_im2jm1 (:)
224 real, allocatable :: wc_ip1jm2 (:)
225 real, allocatable :: wc_ip0jm2 (:)
226 real, allocatable :: wc_im1jm2 (:)
227 real, allocatable :: wc_im2jm2 (:)
228
229 real ux, ap1, ap0, am1, am2
230 real uy, bp1, bp0, bm1, bm2
231
232 real lon_cmp(im)
233 real lat_cmp(jm)
234 real q_tmp(irun)
235
236 real pi,d
237 real lam,lam_ip1,lam_ip0,lam_im1,lam_im2
238 real phi,phi_jp1,phi_jp0,phi_jm1,phi_jm2
239 real dl,dp,phi_np,lam_0
240 real lam_geo, lam_cmp
241 real phi_geo, phi_cmp
242 real undef
243 integer im1_cmp,icmp
244 integer jm1_cmp,jcmp
245
246
247
248 = 4.*atan(1.)
249 dl = 2*pi/ im
250 = pi/(jm-1)
251
252
253
254 allocate ( wl_ip0jp0(irun) , wl_im1jp0(irun) )
255 allocate ( wl_ip0jm1(irun) , wl_im1jm1(irun) )
256 allocate ( wc_ip1jp1(irun) , wc_ip0jp1(irun) , wc_im1jp1(irun) , wc_im2
257 allocate ( wc_ip1jp0(irun) , wc_ip0jp0(irun) , wc_im1jp0(irun) , wc_im2
258 allocate ( wc_ip1jm1(irun) , wc_ip0jm1(irun) , wc_im1jm1(irun) , wc_im2
259 allocate ( wc_ip1jm2(irun) , wc_ip0jm2(irun) , wc_im1jm2(irun) , wc_im2
260 allocate ( ip1(irun) , ip0(irun) , im1(irun) ,
261 allocate ( jp1(irun) , jp0(irun) , jm1(irun) ,
262
263
264
265 (1) = lon_min
266 do i=2,im
267 lon_cmp(i) = lon_cmp(i-1) + dlam(i-1)
268 enddo
269 lat_cmp(1) = -pi*0.5
270 do j=2,jm-1
271 lat_cmp(j) = lat_cmp(j-1) + dphi(j-1)
272 enddo
273 lat_cmp(jm) = pi*0.5
274
275
276
277 do i=1,irun
278 lam_cmp = lon_geo(i)
279 phi_cmp = lat_geo(i)
280
281
282
283 = 1
284 do icmp = 2,im
285 if( lon_cmp(icmp).lt.lam_cmp ) im1_cmp = icmp
286 enddo
287 jm1_cmp = 1
288 do jcmp = 2,jm
289 if( lat_cmp(jcmp).lt.phi_cmp ) jm1_cmp = jcmp
290 enddo
291
292 im1(i) = im1_cmp
293 ip0(i) = im1(i) + 1
294 ip1(i) = ip0(i) + 1
295 im2(i) = im1(i) - 1
296
297 jm1(i) = jm1_cmp
298 jp0(i) = jm1(i) + 1
299 jp1(i) = jp0(i) + 1
300 jm2(i) = jm1(i) - 1
301
302
303
304 if(im1(i).eq.im) then
305 ip0(i) = 1
306 ip1(i) = 2
307 endif
308 if(im1(i).eq.1) then
309 im2(i) = im
310 endif
311 if(ip0(i).eq.im) then
312 ip1(i) = 1
313 endif
314
315
316
317
318 = lam_cmp
319 phi = phi_cmp
320
321
322
323 = lon_cmp(im2(i))
324 lam_im1 = lon_cmp(im1(i))
325 lam_ip0 = lon_cmp(ip0(i))
326 lam_ip1 = lon_cmp(ip1(i))
327
328 if( lam_im2.gt.lam_im1 ) lam_im2 = lam_im2 - 2*pi
329 if( lam_im1.gt.lam_ip0 ) lam_ip0 = lam_ip0 + 2*pi
330 if( lam_im1.gt.lam_ip1 ) lam_ip1 = lam_ip1 + 2*pi
331 if( lam_ip0.gt.lam_ip1 ) lam_ip1 = lam_ip1 + 2*pi
332
333
334
335
336
337 = lat_cmp(jm1(i))
338
339 if( jm2(i).eq.0 ) then
340 phi_jm2 = phi_jm1 - dphi(1)
341 else
342 phi_jm2 = lat_cmp(jm2(i))
343 endif
344
345 if( jm1(i).eq.jm ) then
346 phi_jp0 = phi_jm1 + dphi(jm-1)
347 phi_jp1 = phi_jp0 + dphi(jm-2)
348 else
349 phi_jp0 = lat_cmp(jp0(i))
350 if( jp1(i).eq.jm+1 ) then
351 phi_jp1 = phi_jp0 + dphi(jm-1)
352 else
353 phi_jp1 = lat_cmp(jp1(i))
354 endif
355 endif
356
357
358
359
360 = (lam_ip0-lam_im1)*(phi_jp0-phi_jm1)
361 wl_im1jm1(i) = (lam_ip0-lam )*(phi_jp0-phi )/d
362 wl_ip0jm1(i) = (lam -lam_im1)*(phi_jp0-phi )/d
363 wl_im1jp0(i) = (lam_ip0-lam )*(phi -phi_jm1)/d
364 wl_ip0jp0(i) = (lam -lam_im1)*(phi -phi_jm1)/d
365
366
367
368 = ( (lam -lam_ip0)*(lam -lam_im1)*(lam -lam_im2) )
369 . / ( (lam_ip1-lam_ip0)*(lam_ip1-lam_im1)*(lam_ip1-lam_im2) )
370 ap0 = ( (lam_ip1-lam )*(lam -lam_im1)*(lam -lam_im2) )
371 . / ( (lam_ip1-lam_ip0)*(lam_ip0-lam_im1)*(lam_ip0-lam_im2) )
372 am1 = ( (lam_ip1-lam )*(lam_ip0-lam )*(lam -lam_im2) )
373 . / ( (lam_ip1-lam_im1)*(lam_ip0-lam_im1)*(lam_im1-lam_im2) )
374 am2 = ( (lam_ip1-lam )*(lam_ip0-lam )*(lam_im1-lam ) )
375 . / ( (lam_ip1-lam_im2)*(lam_ip0-lam_im2)*(lam_im1-lam_im2) )
376
377 bp1 = ( (phi -phi_jp0)*(phi -phi_jm1)*(phi -phi_jm2) )
378 . / ( (phi_jp1-phi_jp0)*(phi_jp1-phi_jm1)*(phi_jp1-phi_jm2) )
379 bp0 = ( (phi_jp1-phi )*(phi -phi_jm1)*(phi -phi_jm2) )
380 . / ( (phi_jp1-phi_jp0)*(phi_jp0-phi_jm1)*(phi_jp0-phi_jm2) )
381 bm1 = ( (phi_jp1-phi )*(phi_jp0-phi )*(phi -phi_jm2) )
382 . / ( (phi_jp1-phi_jm1)*(phi_jp0-phi_jm1)*(phi_jm1-phi_jm2) )
383 bm2 = ( (phi_jp1-phi )*(phi_jp0-phi )*(phi_jm1-phi ) )
384 . / ( (phi_jp1-phi_jm2)*(phi_jp0-phi_jm2)*(phi_jm1-phi_jm2) )
385
386 wc_ip1jp1(i) = bp1*ap1
387 wc_ip0jp1(i) = bp1*ap0
388 wc_im1jp1(i) = bp1*am1
389 wc_im2jp1(i) = bp1*am2
390
391 wc_ip1jp0(i) = bp0*ap1
392 wc_ip0jp0(i) = bp0*ap0
393 wc_im1jp0(i) = bp0*am1
394 wc_im2jp0(i) = bp0*am2
395
396 wc_ip1jm1(i) = bm1*ap1
397 wc_ip0jm1(i) = bm1*ap0
398 wc_im1jm1(i) = bm1*am1
399 wc_im2jm1(i) = bm1*am2
400
401 wc_ip1jm2(i) = bm2*ap1
402 wc_ip0jm2(i) = bm2*ap0
403 wc_im1jm2(i) = bm2*am1
404 wc_im2jm2(i) = bm2*am2
405
406 enddo
407
408
409
410 do L=1,lm
411 do i=1,irun
412
413 if( lat_geo(i).le.lat_cmp(2) .or.
414 . lat_geo(i).ge.lat_cmp(jm-1) ) then
415
416
417
418 if( q_cmp( im1(i),jm1(i),L ).ne.undef .and.
419 . q_cmp( ip0(i),jm1(i),L ).ne.undef .and.
420 . q_cmp( im1(i),jp0(i),L ).ne.undef .and.
421 . q_cmp( ip0(i),jp0(i),L ).ne.undef ) then
422
423 q_tmp(i) = wl_im1jm1(i) * q_cmp( im1(i),jm1(i),L )
424 . + wl_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L )
425 . + wl_im1jp0(i) * q_cmp( im1(i),jp0(i),L )
426 . + wl_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L )
427
428 else
429 q_tmp(i) = undef
430 endif
431
432 else
433
434
435
436 if( q_cmp( ip1(i),jp0(i),L ).ne.undef .and.
437 . q_cmp( ip0(i),jp0(i),L ).ne.undef .and.
438 . q_cmp( im1(i),jp0(i),L ).ne.undef .and.
439 . q_cmp( im2(i),jp0(i),L ).ne.undef .and.
440
441 . q_cmp( ip1(i),jm1(i),L ).ne.undef .and.
442 . q_cmp( ip0(i),jm1(i),L ).ne.undef .and.
443 . q_cmp( im1(i),jm1(i),L ).ne.undef .and.
444 . q_cmp( im2(i),jm1(i),L ).ne.undef .and.
445
446 . q_cmp( ip1(i),jp1(i),L ).ne.undef .and.
447 . q_cmp( ip0(i),jp1(i),L ).ne.undef .and.
448 . q_cmp( im1(i),jp1(i),L ).ne.undef .and.
449 . q_cmp( im2(i),jp1(i),L ).ne.undef .and.
450
451 . q_cmp( ip1(i),jm2(i),L ).ne.undef .and.
452 . q_cmp( ip0(i),jm2(i),L ).ne.undef .and.
453 . q_cmp( im1(i),jm2(i),L ).ne.undef .and.
454 . q_cmp( im2(i),jm2(i),L ).ne.undef ) then
455
456 q_tmp(i) = wc_ip1jp1(i) * q_cmp( ip1(i),jp1(i),L )
457 . + wc_ip0jp1(i) * q_cmp( ip0(i),jp1(i),L )
458 . + wc_im1jp1(i) * q_cmp( im1(i),jp1(i),L )
459 . + wc_im2jp1(i) * q_cmp( im2(i),jp1(i),L )
460
461 . + wc_ip1jp0(i) * q_cmp( ip1(i),jp0(i),L )
462 . + wc_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L )
463 . + wc_im1jp0(i) * q_cmp( im1(i),jp0(i),L )
464 . + wc_im2jp0(i) * q_cmp( im2(i),jp0(i),L )
465
466 . + wc_ip1jm1(i) * q_cmp( ip1(i),jm1(i),L )
467 . + wc_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L )
468 . + wc_im1jm1(i) * q_cmp( im1(i),jm1(i),L )
469 . + wc_im2jm1(i) * q_cmp( im2(i),jm1(i),L )
470
471 . + wc_ip1jm2(i) * q_cmp( ip1(i),jm2(i),L )
472 . + wc_ip0jm2(i) * q_cmp( ip0(i),jm2(i),L )
473 . + wc_im1jm2(i) * q_cmp( im1(i),jm2(i),L )
474 . + wc_im2jm2(i) * q_cmp( im2(i),jm2(i),L )
475
476 elseif( q_cmp( im1(i),jm1(i),L ).ne.undef .and.
477 . q_cmp( ip0(i),jm1(i),L ).ne.undef .and.
478 . q_cmp( im1(i),jp0(i),L ).ne.undef .and.
479 . q_cmp( ip0(i),jp0(i),L ).ne.undef ) then
480
481 q_tmp(i) = wl_im1jm1(i) * q_cmp( im1(i),jm1(i),L )
482 . + wl_ip0jm1(i) * q_cmp( ip0(i),jm1(i),L )
483 . + wl_im1jp0(i) * q_cmp( im1(i),jp0(i),L )
484 . + wl_ip0jp0(i) * q_cmp( ip0(i),jp0(i),L )
485
486 else
487 q_tmp(i) = undef
488 endif
489
490 endif
491 enddo
492
493
494
495 do i=1,irun
496 q_geo(i,L) = q_tmp(i)
497 enddo
498 enddo
499
500 deallocate ( wl_ip0jp0 , wl_im1jp0 )
501 deallocate ( wl_ip0jm1 , wl_im1jm1 )
502 deallocate ( wc_ip1jp1 , wc_ip0jp1 , wc_im1jp1 , wc_im2jp1 )
503 deallocate ( wc_ip1jp0 , wc_ip0jp0 , wc_im1jp0 , wc_im2jp0 )
504 deallocate ( wc_ip1jm1 , wc_ip0jm1 , wc_im1jm1 , wc_im2jm1 )
505 deallocate ( wc_ip1jm2 , wc_ip0jm2 , wc_im1jm2 , wc_im2jm2 )
506 deallocate ( ip1 , ip0 , im1 , im2 )
507 deallocate ( jp1 , jp0 , jm1 , jm2 )
508
509 return
510
511 end
512