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