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