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