NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
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
22 ! Smallest positive real value (use for equality comparisons)
23 REAL :: TINYREAL=tiny(1.0)
24
25contains
26
80 SUBROUTINE interpolate_bicubic_scalar(IPOPT,grid_in,grid_out, &
81 MI,MO,KM,IBI,LI,GI, &
82 NO,RLAT,RLON,IBO,LO,GO,IRET)
83 class(ip_grid), intent(in) :: grid_in, grid_out
84 INTEGER, INTENT(IN ) :: IPOPT(20)
85 INTEGER, INTENT(IN ) :: MI,MO,KM
86 INTEGER, INTENT(IN ) :: IBI(KM)
87 INTEGER, INTENT(INOUT) :: NO
88 INTEGER, INTENT( OUT) :: IRET, IBO(KM)
89 !
90 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
91 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
92 !
93 REAL, INTENT(IN ) :: GI(MI,KM)
94 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
95 REAL, INTENT( OUT) :: GO(MO,KM)
96 !
97 REAL, PARAMETER :: FILL=-9999.
98 !
99 INTEGER :: IJX(4),IJY(4)
100 INTEGER :: MCON,MP,N,I,J,K
101 INTEGER :: NK,NV
102 LOGICAL :: SAME_GRIDI, SAME_GRIDO
103 !
104 REAL :: PMP,XIJ,YIJ,XF,YF
105 REAL :: G,W,GMIN,GMAX
106 REAL :: WX(4),WY(4)
107 REAL :: XPTS(MO),YPTS(MO)
108 logical :: to_station_points
109
110 ! Save coeffecients between calls and only compute if grids have changed
111 REAL, ALLOCATABLE,SAVE :: RLATX(:),RLONX(:)
112 REAL, ALLOCATABLE,SAVE :: WXY(:,:,:)
113 INTEGER, SAVE :: NOX=-1,iretx=-1
114 INTEGER, ALLOCATABLE,SAVE :: NXY(:,:,:),NC(:)
115 class(ip_grid), allocatable,save :: prev_grid_in, prev_grid_out
116
117 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118 ! SET PARAMETERS
119 iret=0
120 mcon=ipopt(1)
121 mp=ipopt(2)
122 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
123 IF(mp.LT.0.OR.mp.GT.100) iret=32
124 pmp=mp*0.01
125
126 if (.not. allocated(prev_grid_in) .or. .not. allocated(prev_grid_out)) then
127 allocate(prev_grid_in, source = grid_in)
128 allocate(prev_grid_out, source = grid_out)
129
130 same_gridi = .false.
131 same_grido = .false.
132 else
133 same_gridi = grid_in == prev_grid_in
134 same_grido = grid_out == prev_grid_out
135
136 if (.not. same_gridi .or. .not. same_grido) then
137 deallocate(prev_grid_in)
138 deallocate(prev_grid_out)
139
140 allocate(prev_grid_in, source = grid_in)
141 allocate(prev_grid_out, source = grid_out)
142 end if
143 end if
144
145 select type(grid_out)
146 type is(ip_station_points_grid)
147 to_station_points = .true.
148 class default
149 to_station_points = .false.
150 end select
151
152 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153 ! SAVE OR SKIP WEIGHT COMPUTATION
154 IF(iret.EQ.0.AND.(to_station_points.OR..NOT.same_gridi.OR..NOT.same_grido))THEN
155 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
157 CALL gdswzd(grid_out,0,mo,fill,xpts,ypts,rlon,rlat,no)
158 IF(no.EQ.0) iret=3
159 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160 ! LOCATE INPUT POINTS
161 CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
162 IF(iret.EQ.0.AND.nv.EQ.0) iret=2
163 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164 ! ALLOCATE AND SAVE GRID DATA
165 IF(nox.NE.no) THEN
166 IF(nox.GE.0) DEALLOCATE(rlatx,rlonx,nc,nxy,wxy)
167 ALLOCATE(rlatx(no),rlonx(no),nc(no),nxy(4,4,no),wxy(4,4,no))
168 nox=no
169 ENDIF
170 iretx=iret
171 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172 ! COMPUTE WEIGHTS
173 IF(iret.EQ.0) THEN
174 !$OMP PARALLEL DO PRIVATE(N,XIJ,YIJ,IJX,IJY,XF,YF,J,I,WX,WY) SCHEDULE(STATIC)
175 DO n=1,no
176 rlonx(n)=rlon(n)
177 rlatx(n)=rlat(n)
178 xij=xpts(n)
179 yij=ypts(n)
180 IF(abs(xij-fill).GT.tinyreal.AND.abs(yij-fill).GT.tinyreal) THEN
181 ijx(1:4)=floor(xij-1)+(/0,1,2,3/)
182 ijy(1:4)=floor(yij-1)+(/0,1,2,3/)
183 xf=xij-ijx(2)
184 yf=yij-ijy(2)
185 DO j=1,4
186 DO i=1,4
187 nxy(i,j,n) = grid_in%field_pos(ijx(i), ijy(j))
188 ENDDO
189 ENDDO
190 IF(minval(nxy(1:4,1:4,n)).GT.0) THEN
191 ! BICUBIC WHERE 16-POINT STENCIL IS AVAILABLE
192 nc(n)=1
193 wx(1)=xf*(1-xf)*(2-xf)/(-6.)
194 wx(2)=(xf+1)*(1-xf)*(2-xf)/2.
195 wx(3)=(xf+1)*xf*(2-xf)/2.
196 wx(4)=(xf+1)*xf*(1-xf)/(-6.)
197 wy(1)=yf*(1-yf)*(2-yf)/(-6.)
198 wy(2)=(yf+1)*(1-yf)*(2-yf)/2.
199 wy(3)=(yf+1)*yf*(2-yf)/2.
200 wy(4)=(yf+1)*yf*(1-yf)/(-6.)
201 ELSE
202 ! BILINEAR ELSEWHERE NEAR THE EDGE OF THE GRID
203 nc(n)=2
204 wx(1)=0
205 wx(2)=(1-xf)
206 wx(3)=xf
207 wx(4)=0
208 wy(1)=0
209 wy(2)=(1-yf)
210 wy(3)=yf
211 wy(4)=0
212 ENDIF
213 DO j=1,4
214 DO i=1,4
215 wxy(i,j,n)=wx(i)*wy(j)
216 ENDDO
217 ENDDO
218 ELSE
219 nc(n)=0
220 ENDIF
221 ENDDO
222 ENDIF
223 ENDIF
224 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225 ! INTERPOLATE OVER ALL FIELDS
226 IF(iret.EQ.0.AND.iretx.EQ.0) THEN
227 IF(.not. to_station_points) THEN
228 no=nox
229 DO n=1,no
230 rlon(n)=rlonx(n)
231 rlat(n)=rlatx(n)
232 ENDDO
233 ENDIF
234 !$OMP PARALLEL DO PRIVATE(NK,K,N,G,W,GMIN,GMAX,J,I) SCHEDULE(STATIC)
235 DO nk=1,no*km
236 k=(nk-1)/no+1
237 n=nk-no*(k-1)
238 IF(nc(n).GT.0) THEN
239 g=0
240 w=0
241 IF(mcon.GT.0) gmin=huge(gmin)
242 IF(mcon.GT.0) gmax=-huge(gmax)
243 DO j=nc(n),5-nc(n)
244 DO i=nc(n),5-nc(n)
245 IF(nxy(i,j,n).GT.0)THEN
246 IF(ibi(k).EQ.0.OR.li(nxy(i,j,n),k))THEN
247 g=g+wxy(i,j,n)*gi(nxy(i,j,n),k)
248 w=w+wxy(i,j,n)
249 IF(mcon.GT.0) gmin=min(gmin,gi(nxy(i,j,n),k))
250 IF(mcon.GT.0) gmax=max(gmax,gi(nxy(i,j,n),k))
251 ENDIF
252 ENDIF
253 ENDDO
254 ENDDO
255 lo(n,k)=w.GE.pmp
256 IF(lo(n,k)) THEN
257 go(n,k)=g/w
258 IF(mcon.GT.0) go(n,k)=min(max(go(n,k),gmin),gmax)
259 ELSE
260 go(n,k)=0.
261 ENDIF
262 ELSE
263 lo(n,k)=.false.
264 go(n,k)=0.
265 ENDIF
266 ENDDO
267 DO k=1,km
268 ibo(k)=ibi(k)
269 IF(.NOT.all(lo(1:no,k))) ibo(k)=1
270 ENDDO
271 select type(grid_out)
272 type is(ip_equid_cylind_grid)
273 CALL polfixs(no,mo,km,rlat,ibo,lo,go)
274 end select
275 ELSE
276 IF(iret.EQ.0) iret=iretx
277 IF(.not. to_station_points) no=0
278 ENDIF
279 end subroutine interpolate_bicubic_scalar
280
338 subroutine interpolate_bicubic_vector(ipopt, grid_in, grid_out, &
339 mi, mo, km, ibi, li, ui, vi, &
340 no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret)
341 class(ip_grid), intent(in) :: grid_in, grid_out
342 INTEGER, INTENT(IN ) :: IPOPT(20)
343 INTEGER, INTENT(IN ) :: IBI(KM),MI,MO,KM
344 INTEGER, INTENT(INOUT) :: NO
345 INTEGER, INTENT( OUT) :: IRET, IBO(KM)
346 !
347 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
348 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
349 !
350 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
351 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO),CROT(MO),SROT(MO)
352 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
353 !
354 REAL, PARAMETER :: FILL=-9999.
355 !
356 INTEGER :: IJX(4),IJY(4)
357 INTEGER :: MCON,MP,N,I,J,K,NK,NV
358 !
359 LOGICAL :: SAME_GRIDI,SAME_GRIDO
360 !
361 REAL :: CM,SM,UROT,VROT
362 REAL :: PMP,XIJ,YIJ,XF,YF
363 REAL :: U,V,W,UMIN,UMAX,VMIN,VMAX
364 REAL :: XPTS(MO),YPTS(MO)
365 REAL :: WX(4),WY(4)
366 REAL :: XPTI(MI),YPTI(MI),RLOI(MI),RLAI(MI)
367 REAL :: CROI(MI),SROI(MI)
368
369 logical :: to_station_points
370
371 ! Save coeffecients between calls and only compute if grids have changed
372 REAL, ALLOCATABLE, SAVE :: RLATX(:),RLONX(:),CROTX(:),SROTX(:)
373 REAL, ALLOCATABLE, SAVE :: WXY(:,:,:),CXY(:,:,:),SXY(:,:,:)
374 INTEGER, SAVE :: NOX=-1,iretx=-1
375 INTEGER, ALLOCATABLE, SAVE :: NXY(:,:,:),NC(:)
376 class(ip_grid), allocatable, save :: prev_grid_in, prev_grid_out
377 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
378 ! SET PARAMETERS
379 iret=0
380 mcon=ipopt(1)
381 mp=ipopt(2)
382 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
383 IF(mp.LT.0.OR.mp.GT.100) iret=32
384 pmp=mp*0.01
385
386
387 if (.not. allocated(prev_grid_in) .or. .not. allocated(prev_grid_out)) then
388 allocate(prev_grid_in, source = grid_in)
389 allocate(prev_grid_out, source = grid_out)
390
391 same_gridi = .false.
392 same_grido = .false.
393 else
394 same_gridi = grid_in == prev_grid_in
395 same_grido = grid_out == prev_grid_out
396
397 if (.not. same_gridi .or. .not. same_grido) then
398 deallocate(prev_grid_in)
399 deallocate(prev_grid_out)
400
401 allocate(prev_grid_in, source = grid_in)
402 allocate(prev_grid_out, source = grid_out)
403 end if
404 end if
405
406 select type(grid_out)
407 type is(ip_station_points_grid)
408 to_station_points = .true.
409 class default
410 to_station_points = .false.
411 end select
412
413 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414 ! SAVE OR SKIP WEIGHT COMPUTATION
415 IF(iret.EQ.0.AND.(to_station_points.OR..NOT.same_gridi.OR..NOT.same_grido))THEN
416 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
417 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
418 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
419 IF(no.EQ.0) iret=3
420 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
421 ! LOCATE INPUT POINTS
422 CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
423 IF(iret.EQ.0.AND.nv.EQ.0) iret=2
424 CALL gdswzd(grid_in, 0,mi,fill,xpti,ypti,rloi,rlai,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(abs(xij-fill).GT.tinyreal.AND.abs(yij-fill).GT.tinyreal) 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
movect
subroutine movect(flat, flon, tlat, tlon, crot, srot)
This subprogram provides the rotation parameters to move a vector along a great circle from one posit...
Definition movect.F90:26
bicubic_interp_mod
Bicubic interpolation routines for scalars and vectors.
Definition bicubic_interp_mod.F90:8
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:341
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:83
gdswzd_mod
Driver module for gdswzd routines.
Definition gdswzd_mod.F90:25
ip_grids_mod
Re-export the individual grids.
Definition ip_grids_mod.F90:7
polfix_mod
Make multiple pole scalar values consistent.
Definition polfix_mod.F90:7
polfix_mod::polfixs
subroutine, public polfixs(nm, nx, km, rlat, ib, lo, go)
Make multiple pole scalar values consistent.
Definition polfix_mod.F90:30
polfix_mod::polfixv
subroutine, public polfixv(nm, nx, km, rlat, rlon, ib, lo, uo, vo)
Make multiple pole vector values consistent,.
Definition polfix_mod.F90:125