NCEPLIBS-ip 4.0.0
bicubic_interp_mod.f90
Go to the documentation of this file.
1
4
8module bicubic_interp_mod
9 use gdswzd_mod
10 use polfix_mod
11 use ip_grids_mod
12 implicit none
13
14 private
15 public :: interpolate_bicubic
16
17 interface interpolate_bicubic
18 module procedure interpolate_bicubic_scalar
19 module procedure interpolate_bicubic_vector
20 end interface interpolate_bicubic
21
22contains
23
75 SUBROUTINE interpolate_bicubic_scalar(IPOPT,grid_in,grid_out, &
76 MI,MO,KM,IBI,LI,GI, &
77 NO,RLAT,RLON,IBO,LO,GO,IRET)
78 class(ip_grid), intent(in) :: grid_in, grid_out
79 INTEGER, INTENT(IN ) :: IPOPT(20)
80 INTEGER, INTENT(IN ) :: MI,MO,KM
81 INTEGER, INTENT(IN ) :: IBI(KM)
82 INTEGER, INTENT(INOUT) :: NO
83 INTEGER, INTENT( OUT) :: IRET, IBO(KM)
84 !
85 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
86 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
87 !
88 REAL, INTENT(IN ) :: GI(MI,KM)
89 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
90 REAL, INTENT( OUT) :: GO(MO,KM)
91 !
92 REAL, PARAMETER :: FILL=-9999.
93 !
94 INTEGER :: IJX(4),IJY(4)
95 INTEGER :: MCON,MP,N,I,J,K
96 INTEGER :: NK,NV
97 LOGICAL :: SAME_GRIDI, SAME_GRIDO
98 !
99 REAL :: PMP,XIJ,YIJ,XF,YF
100 REAL :: G,W,GMIN,GMAX
101 REAL :: WX(4),WY(4)
102 REAL :: XPTS(MO),YPTS(MO)
103 logical :: to_station_points
104
105 ! Save coeffecients between calls and only compute if grids have changed
106 REAL, ALLOCATABLE,SAVE :: RLATX(:),RLONX(:)
107 REAL, ALLOCATABLE,SAVE :: WXY(:,:,:)
108 INTEGER, SAVE :: NOX=-1,iretx=-1
109 INTEGER, ALLOCATABLE,SAVE :: NXY(:,:,:),NC(:)
110 class(ip_grid), allocatable,save :: prev_grid_in, prev_grid_out
111
112 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113 ! SET PARAMETERS
114 iret=0
115 mcon=ipopt(1)
116 mp=ipopt(2)
117 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
118 IF(mp.LT.0.OR.mp.GT.100) iret=32
119 pmp=mp*0.01
120
121 if (.not. allocated(prev_grid_in) .or. .not. allocated(prev_grid_out)) then
122 allocate(prev_grid_in, source = grid_in)
123 allocate(prev_grid_out, source = grid_out)
124
125 same_gridi = .false.
126 same_grido = .false.
127 else
128 same_gridi = grid_in == prev_grid_in
129 same_grido = grid_out == prev_grid_out
130
131 if (.not. same_gridi .or. .not. same_grido) then
132 deallocate(prev_grid_in)
133 deallocate(prev_grid_out)
134
135 allocate(prev_grid_in, source = grid_in)
136 allocate(prev_grid_out, source = grid_out)
137 end if
138 end if
139
140 select type(grid_out)
141 type is(ip_station_points_grid)
142 to_station_points = .true.
143 class default
144 to_station_points = .false.
145 end select
146
147 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148 ! SAVE OR SKIP WEIGHT COMPUTATION
149 IF(iret.EQ.0.AND.(to_station_points.OR..NOT.same_gridi.OR..NOT.same_grido))THEN
150 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
152 IF(.not. to_station_points) THEN
153 CALL gdswzd(grid_out,0,mo,fill,xpts,ypts,rlon,rlat,no)
154 IF(no.EQ.0) iret=3
155 ENDIF
156 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157 ! LOCATE INPUT POINTS
158 CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
159 IF(iret.EQ.0.AND.nv.EQ.0) iret=2
160 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161 ! ALLOCATE AND SAVE GRID DATA
162 IF(nox.NE.no) THEN
163 IF(nox.GE.0) DEALLOCATE(rlatx,rlonx,nc,nxy,wxy)
164 ALLOCATE(rlatx(no),rlonx(no),nc(no),nxy(4,4,no),wxy(4,4,no))
165 nox=no
166 ENDIF
167 iretx=iret
168 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169 ! COMPUTE WEIGHTS
170 IF(iret.EQ.0) THEN
171 !$OMP PARALLEL DO PRIVATE(N,XIJ,YIJ,IJX,IJY,XF,YF,J,I,WX,WY) SCHEDULE(STATIC)
172 DO n=1,no
173 rlonx(n)=rlon(n)
174 rlatx(n)=rlat(n)
175 xij=xpts(n)
176 yij=ypts(n)
177 IF(xij.NE.fill.AND.yij.NE.fill) THEN
178 ijx(1:4)=floor(xij-1)+(/0,1,2,3/)
179 ijy(1:4)=floor(yij-1)+(/0,1,2,3/)
180 xf=xij-ijx(2)
181 yf=yij-ijy(2)
182 DO j=1,4
183 DO i=1,4
184 nxy(i,j,n) = grid_in%field_pos(ijx(i), ijy(j))
185 ENDDO
186 ENDDO
187 IF(minval(nxy(1:4,1:4,n)).GT.0) THEN
188 ! BICUBIC WHERE 16-POINT STENCIL IS AVAILABLE
189 nc(n)=1
190 wx(1)=xf*(1-xf)*(2-xf)/(-6.)
191 wx(2)=(xf+1)*(1-xf)*(2-xf)/2.
192 wx(3)=(xf+1)*xf*(2-xf)/2.
193 wx(4)=(xf+1)*xf*(1-xf)/(-6.)
194 wy(1)=yf*(1-yf)*(2-yf)/(-6.)
195 wy(2)=(yf+1)*(1-yf)*(2-yf)/2.
196 wy(3)=(yf+1)*yf*(2-yf)/2.
197 wy(4)=(yf+1)*yf*(1-yf)/(-6.)
198 ELSE
199 ! BILINEAR ELSEWHERE NEAR THE EDGE OF THE GRID
200 nc(n)=2
201 wx(1)=0
202 wx(2)=(1-xf)
203 wx(3)=xf
204 wx(4)=0
205 wy(1)=0
206 wy(2)=(1-yf)
207 wy(3)=yf
208 wy(4)=0
209 ENDIF
210 DO j=1,4
211 DO i=1,4
212 wxy(i,j,n)=wx(i)*wy(j)
213 ENDDO
214 ENDDO
215 ELSE
216 nc(n)=0
217 ENDIF
218 ENDDO
219 ENDIF
220 ENDIF
221 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222 ! INTERPOLATE OVER ALL FIELDS
223 IF(iret.EQ.0.AND.iretx.EQ.0) THEN
224 IF(.not. to_station_points) THEN
225 no=nox
226 DO n=1,no
227 rlon(n)=rlonx(n)
228 rlat(n)=rlatx(n)
229 ENDDO
230 ENDIF
231 !$OMP PARALLEL DO PRIVATE(NK,K,N,G,W,GMIN,GMAX,J,I) SCHEDULE(STATIC)
232 DO nk=1,no*km
233 k=(nk-1)/no+1
234 n=nk-no*(k-1)
235 IF(nc(n).GT.0) THEN
236 g=0
237 w=0
238 IF(mcon.GT.0) gmin=huge(gmin)
239 IF(mcon.GT.0) gmax=-huge(gmax)
240 DO j=nc(n),5-nc(n)
241 DO i=nc(n),5-nc(n)
242 IF(nxy(i,j,n).GT.0)THEN
243 IF(ibi(k).EQ.0.OR.li(nxy(i,j,n),k))THEN
244 g=g+wxy(i,j,n)*gi(nxy(i,j,n),k)
245 w=w+wxy(i,j,n)
246 IF(mcon.GT.0) gmin=min(gmin,gi(nxy(i,j,n),k))
247 IF(mcon.GT.0) gmax=max(gmax,gi(nxy(i,j,n),k))
248 ENDIF
249 ENDIF
250 ENDDO
251 ENDDO
252 lo(n,k)=w.GE.pmp
253 IF(lo(n,k)) THEN
254 go(n,k)=g/w
255 IF(mcon.GT.0) go(n,k)=min(max(go(n,k),gmin),gmax)
256 ELSE
257 go(n,k)=0.
258 ENDIF
259 ELSE
260 lo(n,k)=.false.
261 go(n,k)=0.
262 ENDIF
263 ENDDO
264 DO k=1,km
265 ibo(k)=ibi(k)
266 IF(.NOT.all(lo(1:no,k))) ibo(k)=1
267 ENDDO
268 select type(grid_out)
269 type is(ip_equid_cylind_grid)
270 CALL polfixs(no,mo,km,rlat,ibo,lo,go)
271 end select
272 ELSE
273 IF(iret.EQ.0) iret=iretx
274 IF(.not. to_station_points) no=0
275 ENDIF
276 end subroutine interpolate_bicubic_scalar
277
334 subroutine interpolate_bicubic_vector(ipopt, grid_in, grid_out, &
335 mi, mo, km, ibi, li, ui, vi, &
336 no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret)
337 class(ip_grid), intent(in) :: grid_in, grid_out
338 INTEGER, INTENT(IN ) :: IPOPT(20)
339 INTEGER, INTENT(IN ) :: IBI(KM),MI,MO,KM
340 INTEGER, INTENT(INOUT) :: NO
341 INTEGER, INTENT( OUT) :: IRET, IBO(KM)
342 !
343 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
344 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
345 !
346 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
347 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO),CROT(MO),SROT(MO)
348 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
349 !
350 REAL, PARAMETER :: FILL=-9999.
351 !
352 INTEGER :: IJX(4),IJY(4)
353 INTEGER :: MCON,MP,N,I,J,K,NK,NV
354 !
355 LOGICAL :: SAME_GRIDI,SAME_GRIDO
356 !
357 REAL :: CM,SM,UROT,VROT
358 REAL :: PMP,XIJ,YIJ,XF,YF
359 REAL :: U,V,W,UMIN,UMAX,VMIN,VMAX
360 REAL :: XPTS(MO),YPTS(MO)
361 REAL :: WX(4),WY(4)
362 REAL :: XPTI(MI),YPTI(MI),RLOI(MI),RLAI(MI)
363 REAL :: CROI(MI),SROI(MI)
364
365 logical :: to_station_points
366
367 ! Save coeffecients between calls and only compute if grids have changed
368 REAL, ALLOCATABLE, SAVE :: RLATX(:),RLONX(:),CROTX(:),SROTX(:)
369 REAL, ALLOCATABLE, SAVE :: WXY(:,:,:),CXY(:,:,:),SXY(:,:,:)
370 INTEGER, SAVE :: NOX=-1,iretx=-1
371 INTEGER, ALLOCATABLE, SAVE :: NXY(:,:,:),NC(:)
372 class(ip_grid), allocatable, save :: prev_grid_in, prev_grid_out
373 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
374 ! SET PARAMETERS
375 iret=0
376 mcon=ipopt(1)
377 mp=ipopt(2)
378 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
379 IF(mp.LT.0.OR.mp.GT.100) iret=32
380 pmp=mp*0.01
381
382
383 if (.not. allocated(prev_grid_in) .or. .not. allocated(prev_grid_out)) then
384 allocate(prev_grid_in, source = grid_in)
385 allocate(prev_grid_out, source = grid_out)
386
387 same_gridi = .false.
388 same_grido = .false.
389 else
390 same_gridi = grid_in == prev_grid_in
391 same_grido = grid_out == prev_grid_out
392
393 if (.not. same_gridi .or. .not. same_grido) then
394 deallocate(prev_grid_in)
395 deallocate(prev_grid_out)
396
397 allocate(prev_grid_in, source = grid_in)
398 allocate(prev_grid_out, source = grid_out)
399 end if
400 end if
401
402 select type(grid_out)
403 type is(ip_station_points_grid)
404 to_station_points = .true.
405 class default
406 to_station_points = .false.
407 end select
408
409 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
410 ! SAVE OR SKIP WEIGHT COMPUTATION
411 IF(iret.EQ.0.AND.(to_station_points.OR..NOT.same_gridi.OR..NOT.same_grido))THEN
412 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
413 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
414 IF(.not. to_station_points) then
415 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat, &
416 no,crot,srot)
417 IF(no.EQ.0) iret=3
418 ENDIF
419 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
420 ! LOCATE INPUT POINTS
421 CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
422 IF(iret.EQ.0.AND.nv.EQ.0) iret=2
423 CALL gdswzd(grid_in, 0,mi,fill,xpti,ypti,rloi,rlai, &
424 nv,croi,sroi)
425 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
426 ! ALLOCATE AND SAVE GRID DATA
427 IF(nox.NE.no) THEN
428 IF(nox.GE.0) DEALLOCATE(rlatx,rlonx,crotx,srotx,nc,nxy,wxy,cxy,sxy)
429 ALLOCATE(rlatx(no),rlonx(no),crotx(no),srotx(no),nc(no), &
430 nxy(4,4,no),wxy(4,4,no),cxy(4,4,no),sxy(4,4,no))
431 nox=no
432 ENDIF
433 iretx=iret
434 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
435 ! COMPUTE WEIGHTS
436 IF(iret.EQ.0) THEN
437 !$OMP PARALLEL DO PRIVATE(N,XIJ,YIJ,IJX,IJY,XF,YF,J,I,WX,WY,CM,SM) SCHEDULE(STATIC)
438 DO n=1,no
439 rlonx(n)=rlon(n)
440 rlatx(n)=rlat(n)
441 crotx(n)=crot(n)
442 srotx(n)=srot(n)
443 xij=xpts(n)
444 yij=ypts(n)
445 IF(xij.NE.fill.AND.yij.NE.fill) THEN
446 ijx(1:4)=floor(xij-1)+(/0,1,2,3/)
447 ijy(1:4)=floor(yij-1)+(/0,1,2,3/)
448 xf=xij-ijx(2)
449 yf=yij-ijy(2)
450 DO j=1,4
451 DO i=1,4
452 nxy(i,j,n) = grid_in%field_pos(ijx(i), ijy(j))
453 ENDDO
454 ENDDO
455 IF(minval(nxy(1:4,1:4,n)).GT.0) THEN
456 ! BICUBIC WHERE 16-POINT STENCIL IS AVAILABLE
457 nc(n)=1
458 wx(1)=xf*(1-xf)*(2-xf)/(-6.)
459 wx(2)=(xf+1)*(1-xf)*(2-xf)/2.
460 wx(3)=(xf+1)*xf*(2-xf)/2.
461 wx(4)=(xf+1)*xf*(1-xf)/(-6.)
462 wy(1)=yf*(1-yf)*(2-yf)/(-6.)
463 wy(2)=(yf+1)*(1-yf)*(2-yf)/2.
464 wy(3)=(yf+1)*yf*(2-yf)/2.
465 wy(4)=(yf+1)*yf*(1-yf)/(-6.)
466 ELSE
467 ! BILINEAR ELSEWHERE NEAR THE EDGE OF THE GRID
468 nc(n)=2
469 wx(1)=0
470 wx(2)=(1-xf)
471 wx(3)=xf
472 wx(4)=0
473 wy(1)=0
474 wy(2)=(1-yf)
475 wy(3)=yf
476 wy(4)=0
477 ENDIF
478 DO j=1,4
479 DO i=1,4
480 wxy(i,j,n)=wx(i)*wy(j)
481 IF(nxy(i,j,n).GT.0) THEN
482 CALL movect(rlai(nxy(i,j,n)),rloi(nxy(i,j,n)), &
483 rlat(n),rlon(n),cm,sm)
484 cxy(i,j,n)=cm*croi(nxy(i,j,n))+sm*sroi(nxy(i,j,n))
485 sxy(i,j,n)=sm*croi(nxy(i,j,n))-cm*sroi(nxy(i,j,n))
486 ENDIF
487 ENDDO
488 ENDDO
489 ELSE
490 nc(n)=0
491 ENDIF
492 ENDDO
493 ENDIF
494 ENDIF
495 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
496 ! INTERPOLATE OVER ALL FIELDS
497 IF(iret.EQ.0.AND.iretx.EQ.0) THEN
498 IF(.not. to_station_points) THEN
499 no=nox
500 DO n=1,no
501 rlon(n)=rlonx(n)
502 rlat(n)=rlatx(n)
503 crot(n)=crotx(n)
504 srot(n)=srotx(n)
505 ENDDO
506 ENDIF
507 !$OMP PARALLEL DO PRIVATE(NK,K,N,U,V,W,UMIN,UMAX,VMIN,VMAX,UROT,VROT,J,I) SCHEDULE(STATIC)
508 DO nk=1,no*km
509 k=(nk-1)/no+1
510 n=nk-no*(k-1)
511 IF(nc(n).GT.0) THEN
512 u=0
513 v=0
514 w=0
515 IF(mcon.GT.0) umin=huge(umin)
516 IF(mcon.GT.0) umax=-huge(umax)
517 IF(mcon.GT.0) vmin=huge(vmin)
518 IF(mcon.GT.0) vmax=-huge(vmax)
519 DO j=nc(n),5-nc(n)
520 DO i=nc(n),5-nc(n)
521 IF(nxy(i,j,n).GT.0) THEN
522 IF(ibi(k).EQ.0.OR.li(nxy(i,j,n),k)) THEN
523 urot=cxy(i,j,n)*ui(nxy(i,j,n),k)-sxy(i,j,n)*vi(nxy(i,j,n),k)
524 vrot=sxy(i,j,n)*ui(nxy(i,j,n),k)+cxy(i,j,n)*vi(nxy(i,j,n),k)
525 u=u+wxy(i,j,n)*urot
526 v=v+wxy(i,j,n)*vrot
527 w=w+wxy(i,j,n)
528 IF(mcon.GT.0) umin=min(umin,urot)
529 IF(mcon.GT.0) umax=max(umax,urot)
530 IF(mcon.GT.0) vmin=min(vmin,vrot)
531 IF(mcon.GT.0) vmax=max(vmax,vrot)
532 ENDIF
533 ENDIF
534 ENDDO
535 ENDDO
536 lo(n,k)=w.GE.pmp
537 IF(lo(n,k)) THEN
538 urot=crot(n)*u-srot(n)*v
539 vrot=srot(n)*u+crot(n)*v
540 uo(n,k)=urot/w
541 vo(n,k)=vrot/w
542 IF(mcon.GT.0) uo(n,k)=min(max(uo(n,k),umin),umax)
543 IF(mcon.GT.0) vo(n,k)=min(max(vo(n,k),vmin),vmax)
544 ELSE
545 uo(n,k)=0.
546 vo(n,k)=0.
547 ENDIF
548 ELSE
549 lo(n,k)=.false.
550 uo(n,k)=0.
551 vo(n,k)=0.
552 ENDIF
553 ENDDO
554 DO k=1,km
555 ibo(k)=ibi(k)
556 IF(.NOT.all(lo(1:no,k))) ibo(k)=1
557 ENDDO
558 select type(grid_out)
559 type is(ip_equid_cylind_grid)
560 CALL polfixv(no,mo,km,rlat,rlon,ibo,lo,uo,vo)
561 end select
562 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
563 ELSE
564 IF(iret.EQ.0) iret=iretx
565 IF(.not. to_station_points) no=0
566 ENDIF
567 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
568 end subroutine interpolate_bicubic_vector
569
570end module bicubic_interp_mod
bicubic_interp_mod::interpolate_bicubic
Definition: bicubic_interp_mod.f90:17
gdswzd_mod::gdswzd
Definition: gdswzd_mod.f90:37
bicubic_interp_mod
Bicubic interpolation routines for scalars and vectors.
Definition: bicubic_interp_mod.f90:8
bicubic_interp_mod::interpolate_bicubic_scalar
subroutine interpolate_bicubic_scalar(IPOPT, grid_in, grid_out, MI, MO, KM, IBI, LI, GI, NO, RLAT, RLON, IBO, LO, GO, IRET)
This subprogram performs bicubic interpolation from any grid to any grid for scalar fields.
Definition: bicubic_interp_mod.f90:78
bicubic_interp_mod::interpolate_bicubic_vector
subroutine interpolate_bicubic_vector(ipopt, grid_in, grid_out, mi, mo, km, ibi, li, ui, vi, no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret)
This subprogram performs bicubic interpolation from any grid to any grid for vector fields.
Definition: bicubic_interp_mod.f90:337
gdswzd_mod
Driver module for gdswzd routines.
Definition: gdswzd_mod.f90:25