NCEPLIBS-ip  5.0.0
neighbor_interp_mod.F90
Go to the documentation of this file.
1 
4 
20 module neighbor_interp_mod
21  use gdswzd_mod
22  use polfix_mod
23  use ip_grids_mod
24  implicit none
25 
26  private
27  public :: interpolate_neighbor
28 
29  interface interpolate_neighbor
30  module procedure interpolate_neighbor_scalar
31  module procedure interpolate_neighbor_vector
32  end interface interpolate_neighbor
33 
34  ! Smallest positive real value (use for equality comparisons)
35  REAL :: TINYREAL=tiny(1.0)
36 
37 contains
38 
100  SUBROUTINE interpolate_neighbor_scalar(IPOPT,grid_in,grid_out, &
101  MI,MO,KM,IBI,LI,GI, &
102  NO,RLAT,RLON,IBO,LO,GO,IRET)
103  class(ip_grid), intent(in) :: grid_in, grid_out
104  INTEGER, INTENT(IN ) :: IPOPT(20)
105  INTEGER, INTENT(IN ) :: MI,MO,KM
106  INTEGER, INTENT(IN ) :: IBI(KM)
107  INTEGER, INTENT(INOUT) :: NO
108  INTEGER, INTENT( OUT) :: IRET, IBO(KM)
109  !
110  LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
111  LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
112  !
113  REAL, INTENT(IN ) :: GI(MI,KM)
114  REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
115  REAL, INTENT( OUT) :: GO(MO,KM)
116  !
117  REAL, PARAMETER :: FILL=-9999.
118  !
119  INTEGER :: I1,J1,IXS,JXS
120  INTEGER :: MSPIRAL,N,K,NK
121  INTEGER :: NV
122  INTEGER :: MX,KXS,KXT,IX,JX,NX
123  !
124  LOGICAL :: SAME_GRIDI, SAME_GRIDO
125  !
126  REAL :: XPTS(MO),YPTS(MO)
127  logical :: to_station_points
128 
129  INTEGER, SAVE :: NOX=-1,iretx=-1
130  INTEGER, ALLOCATABLE, SAVE :: NXY(:)
131  REAL, ALLOCATABLE, SAVE :: RLATX(:),RLONX(:),XPTSX(:),YPTSX(:)
132  class(ip_grid), allocatable, save :: prev_grid_in, prev_grid_out
133  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134  ! SET PARAMETERS
135  iret=0
136  mspiral=max(ipopt(1),1)
137  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138  if (.not. allocated(prev_grid_in) .or. .not. allocated(prev_grid_out)) then
139  allocate(prev_grid_in, source = grid_in)
140  allocate(prev_grid_out, source = grid_out)
141 
142  same_gridi = .false.
143  same_grido = .false.
144  else
145  same_gridi = grid_in == prev_grid_in
146  same_grido = grid_out == prev_grid_out
147 
148  if (.not. same_gridi .or. .not. same_grido) then
149  deallocate(prev_grid_in)
150  deallocate(prev_grid_out)
151 
152  allocate(prev_grid_in, source = grid_in)
153  allocate(prev_grid_out, source = grid_out)
154  end if
155  end if
156 
157  select type(grid_out)
158  type is(ip_station_points_grid)
159  to_station_points = .true.
160  class default
161  to_station_points = .false.
162  end select
163  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164  ! SAVE OR SKIP WEIGHT COMPUTATION
165  IF(iret.EQ.0.AND.(to_station_points.OR..NOT.same_gridi.OR..NOT.same_grido))THEN
166  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167  ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
168  CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
169  IF(no.EQ.0) iret=3
170  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171  ! LOCATE INPUT POINTS
172  CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
173  IF(iret.EQ.0.AND.nv.EQ.0) iret=2
174  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175  ! ALLOCATE AND SAVE GRID DATA
176  IF(nox.NE.no) THEN
177  IF(nox.GE.0) DEALLOCATE(rlatx,rlonx,xptsx,yptsx,nxy)
178  ALLOCATE(rlatx(no),rlonx(no),xptsx(no),yptsx(no),nxy(no))
179  nox=no
180  ENDIF
181  iretx=iret
182  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183  ! COMPUTE WEIGHTS
184  IF(iret.EQ.0) THEN
185  !$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC)
186  DO n=1,no
187  rlonx(n)=rlon(n)
188  rlatx(n)=rlat(n)
189  xptsx(n)=xpts(n)
190  yptsx(n)=ypts(n)
191  IF(abs(xpts(n)-fill).GT.tinyreal.AND.abs(ypts(n)-fill).GT.tinyreal) THEN
192  nxy(n) = grid_in%field_pos(nint(xpts(n)), nint(ypts(n)))
193  ELSE
194  nxy(n)=0
195  ENDIF
196  ENDDO
197  ENDIF
198  ENDIF
199  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200  ! INTERPOLATE OVER ALL FIELDS
201  IF(iret.EQ.0.AND.iretx.EQ.0) THEN
202  IF(.not. to_station_points) THEN
203  no=nox
204  DO n=1,no
205  rlon(n)=rlonx(n)
206  rlat(n)=rlatx(n)
207  ENDDO
208  ENDIF
209  DO n=1,no
210  xpts(n)=xptsx(n)
211  ypts(n)=yptsx(n)
212  ENDDO
213  !$OMP PARALLEL DO PRIVATE(NK,K,N,I1,J1,IXS,JXS,MX,KXS,KXT,IX,JX,NX) SCHEDULE(STATIC)
214  DO nk=1,no*km
215  k=(nk-1)/no+1
216  n=nk-no*(k-1)
217  go(n,k)=0
218  lo(n,k)=.false.
219  IF(nxy(n).GT.0) THEN
220  IF(ibi(k).EQ.0.OR.li(nxy(n),k)) THEN
221  go(n,k)=gi(nxy(n),k)
222  lo(n,k)=.true.
223  ! SPIRAL AROUND UNTIL VALID DATA IS FOUND.
224  ELSEIF(mspiral.GT.1) THEN
225  i1=nint(xpts(n))
226  j1=nint(ypts(n))
227  ixs=int(sign(1.,xpts(n)-i1))
228  jxs=int(sign(1.,ypts(n)-j1))
229  DO mx=2,mspiral**2
230  kxs=int(sqrt(4*mx-2.5))
231  kxt=mx-(kxs**2/4+1)
232  SELECT CASE(mod(kxs,4))
233  CASE(1)
234  ix=i1-ixs*(kxs/4-kxt)
235  jx=j1-jxs*kxs/4
236  CASE(2)
237  ix=i1+ixs*(1+kxs/4)
238  jx=j1-jxs*(kxs/4-kxt)
239  CASE(3)
240  ix=i1+ixs*(1+kxs/4-kxt)
241  jx=j1+jxs*(1+kxs/4)
242  CASE DEFAULT
243  ix=i1-ixs*kxs/4
244  jx=j1+jxs*(kxs/4-kxt)
245  END SELECT
246  nx = grid_in%field_pos(ix, jx)
247  IF(nx.GT.0) THEN
248  IF(li(nx,k)) THEN
249  go(n,k)=gi(nx,k)
250  lo(n,k)=.true.
251  EXIT
252  ENDIF
253  ENDIF
254  ENDDO
255  ENDIF
256  ENDIF
257  ENDDO
258 
259  DO k=1,km
260  ibo(k)=ibi(k)
261  IF(.NOT.all(lo(1:no,k))) ibo(k)=1
262  ENDDO
263 
264  select type(grid_out)
265  type is(ip_equid_cylind_grid)
266  CALL polfixs(no,mo,km,rlat,ibo,lo,go)
267  end select
268  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269  ELSE
270  IF(iret.EQ.0) iret=iretx
271  IF(.not. to_station_points) no=0
272  ENDIF
273  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274  END SUBROUTINE interpolate_neighbor_scalar
275 
352  SUBROUTINE interpolate_neighbor_vector(IPOPT,grid_in,grid_out, &
353  MI,MO,KM,IBI,LI,UI,VI, &
354  NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
355 
356  class(ip_grid), intent(in) :: grid_in, grid_out
357  INTEGER, INTENT(IN ) :: IPOPT(20)
358  INTEGER, INTENT(IN ) :: IBI(KM),MI,MO,KM
359  INTEGER, INTENT(INOUT) :: NO
360  INTEGER, INTENT( OUT) :: IRET, IBO(KM)
361  !
362  LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
363  LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
364  !
365  REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
366  REAL, INTENT(INOUT) :: CROT(MO),SROT(MO)
367  REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
368  REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
369  !
370  REAL, PARAMETER :: FILL=-9999.
371  !
372  INTEGER :: I1,J1,IXS,JXS,MX
373  INTEGER :: KXS,KXT,IX,JX,NX
374  INTEGER :: MSPIRAL,N,K,NK,NV
375  !
376  LOGICAL :: SAME_GRIDI, SAME_GRIDO
377  !
378  REAL :: CX,SX,CM,SM,UROT,VROT
379  REAL :: XPTS(MO),YPTS(MO)
380  REAL :: CROI(MI),SROI(MI)
381  REAL :: XPTI(MI),YPTI(MI),RLOI(MI),RLAI(MI)
382 
383  logical :: to_station_points
384 
385  INTEGER, SAVE :: NOX=-1,iretx=-1
386  INTEGER, ALLOCATABLE, SAVE :: NXY(:)
387 
388  REAL, ALLOCATABLE, SAVE :: RLATX(:),RLONX(:),XPTSX(:),YPTSX(:)
389  REAL, ALLOCATABLE, SAVE :: CROTX(:),SROTX(:),CXY(:),SXY(:)
390  class(ip_grid), allocatable, save :: prev_grid_in, prev_grid_out
391  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
392  ! SET PARAMETERS
393  iret=0
394  mspiral=max(ipopt(1),1)
395  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
396 
397  if (.not. allocated(prev_grid_in) .or. .not. allocated(prev_grid_out)) then
398  allocate(prev_grid_in, source = grid_in)
399  allocate(prev_grid_out, source = grid_out)
400 
401  same_gridi = .false.
402  same_grido = .false.
403  else
404  same_gridi = grid_in == prev_grid_in
405  same_grido = grid_out == prev_grid_out
406 
407  if (.not. same_gridi .or. .not. same_grido) then
408  deallocate(prev_grid_in)
409  deallocate(prev_grid_out)
410 
411  allocate(prev_grid_in, source = grid_in)
412  allocate(prev_grid_out, source = grid_out)
413  end if
414  end if
415 
416  select type(grid_out)
417  type is(ip_station_points_grid)
418  to_station_points = .true.
419  class default
420  to_station_points = .false.
421  end select
422 
423  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
424  ! SAVE OR SKIP WEIGHT COMPUTATION
425  IF(iret.EQ.0.AND.(to_station_points.OR..NOT.same_gridi.OR..NOT.same_grido))THEN
426  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
427  ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
428  CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
429  IF(no.EQ.0) iret=3
430  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
431  ! LOCATE INPUT POINTS
432  CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
433  IF(iret.EQ.0.AND.nv.EQ.0) iret=2
434  CALL gdswzd(grid_in, 0,mi,fill,xpti,ypti,rloi,rlai,nv,croi,sroi)
435  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
436  ! ALLOCATE AND SAVE GRID DATA
437  IF(nox.NE.no) THEN
438  IF(nox.GE.0) DEALLOCATE(rlatx,rlonx,xptsx,yptsx,crotx,srotx,nxy,cxy,sxy)
439  ALLOCATE(rlatx(no),rlonx(no),xptsx(no),yptsx(no), &
440  crotx(no),srotx(no),nxy(no),cxy(no),sxy(no))
441  nox=no
442  ENDIF
443  iretx=iret
444  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
445  ! COMPUTE WEIGHTS
446  IF(iret.EQ.0) THEN
447  !$OMP PARALLEL DO PRIVATE(N,CM,SM) SCHEDULE(STATIC)
448  DO n=1,no
449  rlonx(n)=rlon(n)
450  rlatx(n)=rlat(n)
451  xptsx(n)=xpts(n)
452  yptsx(n)=ypts(n)
453  crotx(n)=crot(n)
454  srotx(n)=srot(n)
455  IF(abs(xpts(n)-fill).GT.tinyreal.AND.abs(ypts(n)-fill).GT.tinyreal) THEN
456  nxy(n) = grid_in%field_pos(nint(xpts(n)),nint(ypts(n)))
457  IF(nxy(n).GT.0) THEN
458  CALL movect(rlai(nxy(n)),rloi(nxy(n)),rlat(n),rlon(n),cm,sm)
459  cxy(n)=cm*croi(nxy(n))+sm*sroi(nxy(n))
460  sxy(n)=sm*croi(nxy(n))-cm*sroi(nxy(n))
461  ENDIF
462  ELSE
463  nxy(n)=0
464  ENDIF
465  ENDDO
466  ENDIF
467  ENDIF
468  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
469  ! INTERPOLATE OVER ALL FIELDS
470  IF(iret.EQ.0.AND.iretx.EQ.0) THEN
471  IF(.not. to_station_points) THEN
472  no=nox
473  DO n=1,no
474  rlon(n)=rlonx(n)
475  rlat(n)=rlatx(n)
476  crot(n)=crotx(n)
477  srot(n)=srotx(n)
478  ENDDO
479  ENDIF
480  DO n=1,no
481  xpts(n)=xptsx(n)
482  ypts(n)=yptsx(n)
483  ENDDO
484  !$OMP PARALLEL DO &
485  !$OMP PRIVATE(NK,K,N,I1,J1,IXS,JXS,MX,KXS,KXT,IX,JX,NX) &
486  !$OMP PRIVATE(CM,SM,CX,SX,UROT,VROT) SCHEDULE(STATIC)
487  DO nk=1,no*km
488  k=(nk-1)/no+1
489  n=nk-no*(k-1)
490  uo(n,k)=0
491  vo(n,k)=0
492  lo(n,k)=.false.
493  IF(nxy(n).GT.0) THEN
494  IF(ibi(k).EQ.0.OR.li(nxy(n),k)) THEN
495  urot=cxy(n)*ui(nxy(n),k)-sxy(n)*vi(nxy(n),k)
496  vrot=sxy(n)*ui(nxy(n),k)+cxy(n)*vi(nxy(n),k)
497  uo(n,k)=crot(n)*urot-srot(n)*vrot
498  vo(n,k)=srot(n)*urot+crot(n)*vrot
499  lo(n,k)=.true.
500  ! SPIRAL AROUND UNTIL VALID DATA IS FOUND.
501  ELSEIF(mspiral.GT.1) THEN
502  i1=nint(xpts(n))
503  j1=nint(ypts(n))
504  ixs=int(sign(1.,xpts(n)-i1))
505  jxs=int(sign(1.,ypts(n)-j1))
506  DO mx=2,mspiral**2
507  kxs=int(sqrt(4*mx-2.5))
508  kxt=mx-(kxs**2/4+1)
509  SELECT CASE(mod(kxs,4))
510  CASE(1)
511  ix=i1-ixs*(kxs/4-kxt)
512  jx=j1-jxs*kxs/4
513  CASE(2)
514  ix=i1+ixs*(1+kxs/4)
515  jx=j1-jxs*(kxs/4-kxt)
516  CASE(3)
517  ix=i1+ixs*(1+kxs/4-kxt)
518  jx=j1+jxs*(1+kxs/4)
519  CASE DEFAULT
520  ix=i1-ixs*kxs/4
521  jx=j1+jxs*(kxs/4-kxt)
522  END SELECT
523  nx = grid_in%field_pos(ix, jx)
524  IF(nx.GT.0) THEN
525  IF(li(nx,k)) THEN
526  CALL movect(rlai(nx),rloi(nx),rlat(n),rlon(n),cm,sm)
527  cx=cm*croi(nx)+sm*sroi(nx)
528  sx=sm*croi(nx)-cm*sroi(nx)
529  urot=cx*ui(nx,k)-sx*vi(nx,k)
530  vrot=sx*ui(nx,k)+cx*vi(nx,k)
531  uo(n,k)=crot(n)*urot-srot(n)*vrot
532  vo(n,k)=srot(n)*urot+crot(n)*vrot
533  lo(n,k)=.true.
534  EXIT
535  ENDIF
536  ENDIF
537  ENDDO
538  ENDIF
539  ENDIF
540  ENDDO
541  DO k=1,km
542  ibo(k)=ibi(k)
543  IF(.NOT.all(lo(1:no,k))) ibo(k)=1
544  ENDDO
545 
546  select type(grid_out)
547  type is(ip_equid_cylind_grid)
548  CALL polfixv(no,mo,km,rlat,rlon,ibo,lo,uo,vo)
549  end select
550 
551  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
552  ELSE
553  IF(iret.EQ.0) iret=iretx
554  IF(.not. to_station_points) no=0
555  ENDIF
556  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
557  END SUBROUTINE interpolate_neighbor_vector
558 
559 end module neighbor_interp_mod
gdswzd_mod::gdswzd
Definition: gdswzd_mod.F90:37
neighbor_interp_mod::interpolate_neighbor
Definition: neighbor_interp_mod.F90:29
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
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
neighbor_interp_mod
Interpolate scalar fields (neighbor).
Definition: neighbor_interp_mod.F90:20
neighbor_interp_mod::interpolate_neighbor_vector
subroutine interpolate_neighbor_vector(IPOPT, grid_in, grid_out, MI, MO, KM, IBI, LI, UI, VI, NO, RLAT, RLON, CROT, SROT, IBO, LO, UO, VO, IRET)
Interpolate vector fields (neighbor).
Definition: neighbor_interp_mod.F90:355
neighbor_interp_mod::interpolate_neighbor_scalar
subroutine interpolate_neighbor_scalar(IPOPT, grid_in, grid_out, MI, MO, KM, IBI, LI, GI, NO, RLAT, RLON, IBO, LO, GO, IRET)
Interpolate scalar fields (neighbor).
Definition: neighbor_interp_mod.F90:103
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