NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
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
24 ! Smallest positive real value (use for equality comparisons)
25 REAL :: TINYREAL=tiny(1.0)
26
27contains
28
72 subroutine interpolate_bilinear_scalar(IPOPT,grid_in,grid_out,MI,MO,KM,IBI,LI,GI,NO,RLAT,RLON,IBO,LO,GO,IRET)
73 class(ip_grid), intent(in) :: grid_in, grid_out
74 INTEGER, INTENT(IN ) :: IPOPT(20)
75 INTEGER, INTENT(IN ) :: MI,MO,KM
76 INTEGER, INTENT(IN ) :: IBI(KM)
77 INTEGER, INTENT(INOUT) :: NO
78 INTEGER, INTENT( OUT) :: IRET, IBO(KM)
79 !
80 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
81 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
82 !
83 REAL, INTENT(IN ) :: GI(MI,KM)
84 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
85 REAL, INTENT( OUT) :: GO(MO,KM)
86 !
87 REAL, PARAMETER :: FILL=-9999.
88 !
89 INTEGER :: IJX(2),IJY(2)
90 INTEGER :: MP,N,I,J,K
91 INTEGER :: NK,NV
92 INTEGER :: MSPIRAL,I1,J1,IXS,JXS
93 INTEGER :: MX,KXS,KXT,IX,JX,NX
94 !
95 LOGICAL :: SAME_GRIDI, SAME_GRIDO
96 !
97 REAL :: WX(2),WY(2)
98 REAL :: XPTS(MO),YPTS(MO)
99 REAL :: PMP,XIJ,YIJ,XF,YF,G,W
100
101 logical :: to_station_points
102
103 ! Save coeffecients between calls and only compute if grids have changed
104 INTEGER, SAVE :: NOX=-1,iretx=-1
105 INTEGER, ALLOCATABLE,SAVE :: NXY(:,:,:)
106 REAL, ALLOCATABLE,SAVE :: RLATX(:),RLONX(:)
107 REAL, ALLOCATABLE,SAVE :: WXY(:,:,:)
108 class(ip_grid), allocatable,save :: prev_grid_in, prev_grid_out
109
110 iret=0
111 mp=ipopt(1)
112 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
113 IF(mp.LT.0.OR.mp.GT.100) iret=32
114 pmp=mp*0.01
115 mspiral=max(ipopt(2),0)
116
117 if (.not. allocated(prev_grid_in) .or. .not. allocated(prev_grid_out)) then
118 allocate(prev_grid_in, source = grid_in)
119 allocate(prev_grid_out, source = grid_out)
120
121 same_gridi = .false.
122 same_grido = .false.
123 else
124 same_gridi = grid_in == prev_grid_in
125 same_grido = grid_out == prev_grid_out
126
127 if (.not. same_gridi .or. .not. same_grido) then
128 deallocate(prev_grid_in)
129 deallocate(prev_grid_out)
130
131 allocate(prev_grid_in, source = grid_in)
132 allocate(prev_grid_out, source = grid_out)
133 end if
134 end if
135
136 select type(grid_out)
137 type is(ip_station_points_grid)
138 to_station_points = .true.
139 class default
140 to_station_points = .false.
141 end select
142
143 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144 ! SAVE OR SKIP WEIGHT COMPUTATION
145 IF(iret==0.AND.(to_station_points.OR..NOT.same_gridi.OR..NOT.same_grido))THEN
146 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
148 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
149 IF(no.EQ.0) iret=3
150 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151 ! LOCATE INPUT POINTS
152 CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
153 IF(iret.EQ.0.AND.nv.EQ.0) iret=2
154 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155 ! ALLOCATE AND SAVE GRID DATA
156 IF(nox.NE.no) THEN
157 IF(nox.GE.0) DEALLOCATE(rlatx,rlonx,nxy,wxy)
158 ALLOCATE(rlatx(no),rlonx(no),nxy(2,2,no),wxy(2,2,no))
159 nox=no
160 ENDIF
161 iretx=iret
162 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163 ! COMPUTE WEIGHTS
164 IF(iret.EQ.0) THEN
165 !$OMP PARALLEL DO PRIVATE(N,XIJ,YIJ,IJX,IJY,XF,YF,J,I,WX,WY) SCHEDULE(STATIC)
166 DO n=1,no
167 rlonx(n)=rlon(n)
168 rlatx(n)=rlat(n)
169 xij=xpts(n)
170 yij=ypts(n)
171 IF(abs(xij-fill).GT.tinyreal.AND.abs(yij-fill).GT.tinyreal) THEN
172 ijx(1:2)=floor(xij)+(/0,1/)
173 ijy(1:2)=floor(yij)+(/0,1/)
174 xf=xij-ijx(1)
175 yf=yij-ijy(1)
176 wx(1)=(1-xf)
177 wx(2)=xf
178 wy(1)=(1-yf)
179 wy(2)=yf
180 DO j=1,2
181 DO i=1,2
182 nxy(i,j,n)=grid_in%field_pos(ijx(i), ijy(j))
183 wxy(i,j,n)=wx(i)*wy(j)
184 ENDDO
185 ENDDO
186 ELSE
187 nxy(:,:,n)=0
188 ENDIF
189 ENDDO
190 ENDIF
191 ENDIF
192 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193 ! INTERPOLATE OVER ALL FIELDS
194 IF(iret.EQ.0.AND.iretx.EQ.0) THEN
195 IF(.not. to_station_points) THEN
196 no=nox
197 DO n=1,no
198 rlon(n)=rlonx(n)
199 rlat(n)=rlatx(n)
200 ENDDO
201 ENDIF
202 !$OMP PARALLEL DO &
203 !$OMP PRIVATE(NK,K,N,G,W,J,I) &
204 !$OMP PRIVATE(I1,J1,IXS,JXS,MX,KXS,KXT,IX,JX,NX) SCHEDULE(STATIC)
205 DO nk=1,no*km
206 k=(nk-1)/no+1
207 n=nk-no*(k-1)
208 g=0
209 w=0
210 DO j=1,2
211 DO i=1,2
212 IF(nxy(i,j,n).GT.0)THEN
213 IF(ibi(k).EQ.0.OR.li(nxy(i,j,n),k)) THEN
214 g=g+wxy(i,j,n)*gi(nxy(i,j,n),k)
215 w=w+wxy(i,j,n)
216 ENDIF
217 ENDIF
218 ENDDO
219 ENDDO
220 lo(n,k)=w.GE.pmp
221 IF(lo(n,k)) THEN
222 go(n,k)=g/w
223 ELSEIF(mspiral.GT.0.AND.abs(xpts(n)-fill).GT.tinyreal.AND.abs(ypts(n)-fill).GT.tinyreal) THEN
224 i1=nint(xpts(n))
225 j1=nint(ypts(n))
226 ixs=int(sign(1.,xpts(n)-i1))
227 jxs=int(sign(1.,ypts(n)-j1))
228 spiral : DO mx=1,mspiral**2
229 kxs=int(sqrt(4*mx-2.5))
230 kxt=mx-(kxs**2/4+1)
231 SELECT CASE(mod(kxs,4))
232 CASE(1)
233 ix=i1-ixs*(kxs/4-kxt)
234 jx=j1-jxs*kxs/4
235 CASE(2)
236 ix=i1+ixs*(1+kxs/4)
237 jx=j1-jxs*(kxs/4-kxt)
238 CASE(3)
239 ix=i1+ixs*(1+kxs/4-kxt)
240 jx=j1+jxs*(1+kxs/4)
241 CASE DEFAULT
242 ix=i1-ixs*kxs/4
243 jx=j1+jxs*(kxs/4-kxt)
244 END SELECT
245 nx=grid_in%field_pos(ix, jx)
246 IF(nx.GT.0.)THEN
247 IF(li(nx,k).OR.ibi(k).EQ.0)THEN
248 go(n,k)=gi(nx,k)
249 lo(n,k)=.true.
250 EXIT spiral
251 ENDIF
252 ENDIF
253 ENDDO spiral
254 IF(.NOT.lo(n,k))THEN
255 ibo(k)=1
256 go(n,k)=0.
257 ENDIF
258 ELSE
259 go(n,k)=0.
260 ENDIF
261 ENDDO
262 DO k=1,km
263 ibo(k)=ibi(k)
264 IF(.NOT.all(lo(1:no,k))) ibo(k)=1
265 ENDDO
266 select type(grid_out)
267 type is(ip_equid_cylind_grid)
268 CALL polfixs(no,mo,km,rlat,ibo,lo,go)
269 end select
270 ELSE
271 IF(iret.EQ.0) iret=iretx
272 IF(.not. to_station_points) no=0
273 ENDIF
274
275 end subroutine interpolate_bilinear_scalar
276
329 SUBROUTINE interpolate_bilinear_vector(ipopt,grid_in,grid_out, &
330 MI,MO,KM,IBI,LI,UI,VI, &
331 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
332 class(ip_grid), intent(in) :: grid_in, grid_out
333 INTEGER, INTENT(IN ) :: IPOPT(20),IBI(KM),MI,MO,KM
334 INTEGER, INTENT(INOUT) :: NO
335 INTEGER, INTENT( OUT) :: IRET, IBO(KM)
336 !
337 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
338 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
339 !
340 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
341 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO),CROT(MO),SROT(MO)
342 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
343 !
344 REAL, PARAMETER :: FILL=-9999.
345 !
346 INTEGER :: IJX(2),IJY(2)
347 INTEGER :: MP,N,I,J,K,NK,NV
348 !
349 LOGICAL :: SAME_GRIDI, SAME_GRIDO
350 !
351 REAL :: CM,SM,UROT,VROT
352 REAL :: PMP,XIJ,YIJ,XF,YF,U,V,W
353 REAL :: XPTS(MO),YPTS(MO)
354 REAL :: WX(2),WY(2)
355 REAL :: XPTI(MI),YPTI(MI)
356 REAL :: RLOI(MI),RLAI(MI)
357 REAL :: CROI(MI),SROI(MI)
358
359 logical :: to_station_points
360
361 ! Save coeffecients between calls and only compute if grids have changed
362 INTEGER, SAVE :: NOX=-1,iretx=-1
363 INTEGER, ALLOCATABLE,SAVE :: NXY(:,:,:)
364 REAL, ALLOCATABLE,SAVE :: RLATX(:),RLONX(:)
365 REAL, ALLOCATABLE,SAVE :: CROTX(:),SROTX(:)
366 REAL, ALLOCATABLE,SAVE :: WXY(:,:,:),CXY(:,:,:),SXY(:,:,:)
367 class(ip_grid), allocatable,save :: prev_grid_in, prev_grid_out
368
369 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
370 ! SET PARAMETERS
371 iret=0
372 mp=ipopt(1)
373 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
374 IF(mp.LT.0.OR.mp.GT.100) iret=32
375 pmp=mp*0.01
376
377 if (.not. allocated(prev_grid_in) .or. .not. allocated(prev_grid_out)) then
378 allocate(prev_grid_in, source = grid_in)
379 allocate(prev_grid_out, source = grid_out)
380
381 same_gridi = .false.
382 same_grido = .false.
383 else
384 same_gridi = grid_in == prev_grid_in
385 same_grido = grid_out == prev_grid_out
386
387 if (.not. same_gridi .or. .not. same_grido) then
388 deallocate(prev_grid_in)
389 deallocate(prev_grid_out)
390
391 allocate(prev_grid_in, source = grid_in)
392 allocate(prev_grid_out, source = grid_out)
393 end if
394 end if
395
396 select type(grid_out)
397 type is(ip_station_points_grid)
398 to_station_points = .true.
399 class default
400 to_station_points = .false.
401 end select
402
403 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
404 ! SAVE OR SKIP WEIGHT COMPUTATION
405 IF(iret.EQ.0.AND.(to_station_points.OR..NOT.same_gridi.OR..NOT.same_grido))THEN
406 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
407 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
408 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
409 IF(no.EQ.0) iret=3
410 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
411 ! LOCATE INPUT POINTS
412 CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
413 IF(iret.EQ.0.AND.nv.EQ.0) iret=2
414 CALL gdswzd(grid_in, 0,mi,fill,xpti,ypti,rloi,rlai,nv,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(abs(xij-fill).GT.tinyreal.AND.abs(yij-fill).GT.tinyreal) 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
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
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:73
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:332
gdswzd_mod
Driver module for gdswzd routines.
Definition gdswzd_mod.F90:25
ip_grid_descriptor_mod
Users derived type grid descriptor objects to abstract away the raw GRIB1 and GRIB2 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
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