NCEPLIBS-ip  4.2.0
bicubic_interp_mod.F90
Go to the documentation of this file.
1 
4 
8 module 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 
25 contains
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 
570 end 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_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
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
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