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