NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
neighbor_interp_mod.F90
Go to the documentation of this file.
1
4
20module 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
37contains
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
559end 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