NCEPLIBS-ip  4.2.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  ! Smallest positive real value (use for equality comparisons)
25  REAL :: TINYREAL=tiny(1.0)
26 
27 contains
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 
522 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: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