NCEPLIBS-ip  4.3.0
neighbor_budget_interp_mod.F90
Go to the documentation of this file.
1 
4 
8 module neighbor_budget_interp_mod
9  use gdswzd_mod
10  use polfix_mod
11  use ip_grids_mod
12  implicit none
13 
14  private
15  public :: interpolate_neighbor_budget
16 
17  interface interpolate_neighbor_budget
18  module procedure interpolate_neighbor_budget_scalar
19  module procedure interpolate_neighbor_budget_vector
20  end interface interpolate_neighbor_budget
21 
22  ! Smallest positive real value (use for equality comparisons)
23  REAL :: TINYREAL=tiny(1.0)
24 
25 contains
26 
107  SUBROUTINE interpolate_neighbor_budget_scalar(IPOPT,grid_in,grid_out, &
108  MI,MO,KM,IBI,LI,GI, &
109  NO,RLAT,RLON,IBO,LO,GO,IRET)
110  class(ip_grid), intent(in) :: grid_in, grid_out
111 
112  INTEGER, INTENT(IN ) :: IBI(KM), IPOPT(20), KM, MI, MO
113  INTEGER, INTENT( OUT) :: IBO(KM), IRET, NO
114  !
115  LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
116  LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
117  !
118  REAL, INTENT(IN ) :: GI(MI,KM)
119  REAL, INTENT( OUT) :: GO(MO,KM), RLAT(MO), RLON(MO)
120  !
121  REAL, PARAMETER :: FILL=-9999.
122  !
123  INTEGER :: IB, I1
124  INTEGER :: JB, J1, K, LB, LSW, MP, N
125  INTEGER :: N11(MO), NB, NB1, NB2, NB3, NB4, NV
126  !
127  REAL :: PMP,RLOB(MO),RLAB(MO)
128  REAL :: WB, WO(MO,KM), XI, YI
129  REAL :: XPTB(MO),YPTB(MO),XPTS(MO),YPTS(MO)
130 
131  logical :: to_station_points
132 
133  select type(grid_out)
134  type is(ip_station_points_grid)
135  to_station_points = .true.
136  class default
137  to_station_points = .false.
138  end select
139 
140  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141  ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
142  iret=0
143  if(to_station_points) then
144  CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
145  IF(no.EQ.0) iret=3
146  CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
147  IF(nv.EQ.0) iret=2
148  else
149  CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
150  IF(no.EQ.0) iret=3
151  endif
152  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153  ! SET PARAMETERS
154  nb1=ipopt(1)
155  IF(nb1.EQ.-1) nb1=2
156  IF(iret.EQ.0.AND.nb1.LT.0) iret=32
157  lsw=1
158  IF(ipopt(1).EQ.-1.OR.ipopt(2).EQ.-1) lsw=0
159  IF(iret.EQ.0.AND.lsw.EQ.1.AND.nb1.GT.15) iret=32
160  mp=ipopt(3+ipopt(1))
161  IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
162  IF(mp.LT.0.OR.mp.GT.100) iret=32
163  pmp=mp*0.01
164  IF(iret.EQ.0) THEN
165  nb2=2*nb1+1
166  nb3=nb2*nb2
167  nb4=nb3
168  IF(lsw.EQ.1) THEN
169  nb4=ipopt(2)
170  DO ib=1,nb1
171  nb4=nb4+8*ib*ipopt(2+ib)
172  ENDDO
173  ENDIF
174  ELSE
175  nb2=0
176  nb3=0
177  nb4=0
178  ENDIF
179  DO k=1,km
180  DO n=1,no
181  go(n,k)=0.
182  wo(n,k)=0.
183  ENDDO
184  ENDDO
185  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186  ! LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX
187 
188  DO nb=1,nb3
189  ! LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS
190  jb=(nb-1)/nb2-nb1
191  ib=nb-(jb+nb1)*nb2-nb1-1
192  lb=max(abs(ib),abs(jb))
193  wb=1
194  IF(lsw.EQ.1) wb=ipopt(2+lb)
195  IF(abs(wb).GT.tinyreal) THEN
196  DO n=1,no
197  xptb(n)=xpts(n)+ib/real(nb2)
198  yptb(n)=ypts(n)+jb/real(nb2)
199  ENDDO
200  if(to_station_points)then
201  CALL gdswzd(grid_in, 1,no,fill,xptb,yptb,rlob,rlab,nv)
202  CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
203  else
204  CALL gdswzd(grid_out, 1,no,fill,xptb,yptb,rlob,rlab,nv)
205  CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
206  endif
207  IF(iret.EQ.0.AND.nv.EQ.0.AND.lb.EQ.0) iret=2
208  DO n=1,no
209  xi=xptb(n)
210  yi=yptb(n)
211  IF(abs(xi-fill).GT.tinyreal.AND.abs(yi-fill).GT.tinyreal) THEN
212  i1=nint(xi)
213  j1=nint(yi)
214  n11(n)=grid_in%field_pos(i1, j1)
215  ELSE
216  n11(n)=0
217  ENDIF
218  ENDDO
219  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220  ! INTERPOLATE WITH OR WITHOUT BITMAPS
221  DO k=1,km
222  DO n=1,no
223  IF(n11(n).GT.0) THEN
224  IF(ibi(k).EQ.0.OR.li(n11(n),k)) THEN
225  go(n,k)=go(n,k)+wb*gi(n11(n),k)
226  wo(n,k)=wo(n,k)+wb
227  ENDIF
228  ENDIF
229  ENDDO
230  ENDDO
231  ENDIF
232  ENDDO
233  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234  ! COMPUTE OUTPUT BITMAPS AND FIELDS
235  DO k=1,km
236  ibo(k)=ibi(k)
237  DO n=1,no
238  lo(n,k)=wo(n,k).GE.pmp*nb4
239  IF(lo(n,k)) THEN
240  go(n,k)=go(n,k)/wo(n,k)
241  ELSE
242  ibo(k)=1
243  go(n,k)=0.
244  ENDIF
245  ENDDO
246  ENDDO
247 
248  select type(grid_out)
249  type is(ip_equid_cylind_grid)
250  CALL polfixs(no,mo,km,rlat,ibo,lo,go)
251  end select
252 
253  END SUBROUTINE interpolate_neighbor_budget_scalar
254 
255 
351  SUBROUTINE interpolate_neighbor_budget_vector(IPOPT,grid_in,grid_out, &
352  MI,MO,KM,IBI,LI,UI,VI, &
353  NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
354  class(ip_grid), intent(in) :: grid_in, grid_out
355 
356  INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
357  INTEGER, INTENT(IN ) :: KM, MI, MO
358  INTEGER, INTENT( OUT) :: IRET, NO, IBO(KM)
359  !
360  LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
361  LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
362  !
363  REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
364  REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
365  REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
366  REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
367  !
368  REAL, PARAMETER :: FILL=-9999.
369  !
370  INTEGER :: N11(MO)
371  INTEGER :: IB, JB, I1, J1
372  INTEGER :: K, LB, LSW, MP, N, NV
373  INTEGER :: NB, NB1, NB2, NB3, NB4
374  !
375  LOGICAL :: SAME_GRID
376  !
377  REAL :: C11(MO),S11(MO)
378  REAL :: CM11, SM11, PMP
379  REAL :: U11, V11, UROT, VROT
380  REAL :: WB, WO(MO,KM), XI, YI
381  REAL :: RLOB(MO),RLAB(MO)
382  REAL :: XPTS(MO),YPTS(MO)
383  REAL :: XPTB(MO),YPTB(MO)
384 
385  logical :: to_station_points
386 
387  ! Save coeffecients between runs and only compute if grid has changed
388  INTEGER, SAVE :: MIX=-1
389  REAL, ALLOCATABLE,SAVE :: CROI(:),SROI(:)
390  REAL, ALLOCATABLE,SAVE :: XPTI(:),YPTI(:)
391  REAL, ALLOCATABLE,SAVE :: RLOI(:),RLAI(:)
392  class(ip_grid), allocatable, save :: prev_grid_in
393 
394  select type(grid_out)
395  type is(ip_station_points_grid)
396  to_station_points = .true.
397  class default
398  to_station_points = .false.
399  end select
400 
401  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
402  ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
403  iret=0
404  CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
405  IF(no.EQ.0) iret=3
406  if(to_station_points) then
407  CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv,crot,srot)
408  IF(nv.EQ.0) iret=2
409  endif
410 
411  if (.not. allocated(prev_grid_in)) then
412  allocate(prev_grid_in, source = grid_in)
413 
414  same_grid = .false.
415  else
416  same_grid = grid_in == prev_grid_in
417 
418  if (.not. same_grid) then
419  deallocate(prev_grid_in)
420  allocate(prev_grid_in, source = grid_in)
421  end if
422  end if
423 
424  IF(.NOT.same_grid) THEN
425  IF(mix.NE.mi) THEN
426  IF(mix.GE.0) DEALLOCATE(xpti,ypti,rloi,rlai,croi,sroi)
427  ALLOCATE(xpti(mi),ypti(mi),rloi(mi),rlai(mi),croi(mi),sroi(mi))
428  mix=mi
429  ENDIF
430  CALL gdswzd(grid_in,0,mi,fill,xpti,ypti, &
431  rloi,rlai,nv,croi,sroi)
432  ENDIF
433  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
434  ! SET PARAMETERS
435  nb1=ipopt(1)
436  IF(nb1.EQ.-1) nb1=2
437  IF(iret.EQ.0.AND.nb1.LT.0) iret=32
438  lsw=1
439  IF(ipopt(1).EQ.-1.OR.ipopt(2).EQ.-1) lsw=0
440  IF(iret.EQ.0.AND.lsw.EQ.1.AND.nb1.GT.15) iret=32
441  mp=ipopt(3+ipopt(1))
442  IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
443  IF(mp.LT.0.OR.mp.GT.100) iret=32
444  pmp=mp*0.01
445  IF(iret.EQ.0) THEN
446  nb2=2*nb1+1
447  nb3=nb2*nb2
448  nb4=nb3
449  IF(lsw.EQ.1) THEN
450  nb4=ipopt(2)
451  DO ib=1,nb1
452  nb4=nb4+8*ib*ipopt(2+ib)
453  ENDDO
454  ENDIF
455  ELSE
456  nb2=0
457  nb3=0
458  nb4=0
459  ENDIF
460  DO k=1,km
461  DO n=1,no
462  uo(n,k)=0
463  vo(n,k)=0
464  wo(n,k)=0.
465  ENDDO
466  ENDDO
467  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
468  ! LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX
469  DO nb=1,nb3
470  ! LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS AND ROTATIONS
471  jb=(nb-1)/nb2-nb1
472  ib=nb-(jb+nb1)*nb2-nb1-1
473  lb=max(abs(ib),abs(jb))
474  wb=1
475  IF(lsw.EQ.1) wb=ipopt(2+lb)
476  IF(abs(wb).GT.tinyreal) THEN
477  DO n=1,no
478  xptb(n)=xpts(n)+ib/real(nb2)
479  yptb(n)=ypts(n)+jb/real(nb2)
480  ENDDO
481  if(to_station_points)then
482  CALL gdswzd(grid_in, 1,no,fill,xptb,yptb,rlob,rlab,nv)
483  CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
484  else
485  CALL gdswzd(grid_out, 1,no,fill,xptb,yptb,rlob,rlab,nv)
486  CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
487  endif
488  IF(iret.EQ.0.AND.nv.EQ.0.AND.lb.EQ.0) iret=2
489  DO n=1,no
490  xi=xptb(n)
491  yi=yptb(n)
492  IF(abs(xi-fill).GT.tinyreal.AND.abs(yi-fill).GT.tinyreal) THEN
493  i1=nint(xi)
494  j1=nint(yi)
495  n11(n)=grid_in%field_pos(i1, j1)
496  IF(n11(n).GT.0) THEN
497  CALL movect(rlai(n11(n)),rloi(n11(n)),rlat(n),rlon(n),cm11,sm11)
498  c11(n)=cm11*croi(n11(n))+sm11*sroi(n11(n))
499  s11(n)=sm11*croi(n11(n))-cm11*sroi(n11(n))
500  ENDIF
501  ELSE
502  n11(n)=0
503  ENDIF
504  ENDDO
505  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
506  ! INTERPOLATE WITH OR WITHOUT BITMAPS
507  DO k=1,km
508  DO n=1,no
509  IF(n11(n).GT.0) THEN
510  IF(ibi(k).EQ.0.OR.li(n11(n),k)) THEN
511  u11=c11(n)*ui(n11(n),k)-s11(n)*vi(n11(n),k)
512  v11=s11(n)*ui(n11(n),k)+c11(n)*vi(n11(n),k)
513  uo(n,k)=uo(n,k)+wb*u11
514  vo(n,k)=vo(n,k)+wb*v11
515  wo(n,k)=wo(n,k)+wb
516  ENDIF
517  ENDIF
518  ENDDO
519  ENDDO
520  ENDIF
521  ENDDO ! NB LOOP
522  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
523  ! COMPUTE OUTPUT BITMAPS AND FIELDS
524  DO k=1,km
525  ibo(k)=ibi(k)
526  DO n=1,no
527  lo(n,k)=wo(n,k).GE.pmp*nb4
528  IF(lo(n,k)) THEN
529  uo(n,k)=uo(n,k)/wo(n,k)
530  vo(n,k)=vo(n,k)/wo(n,k)
531  urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
532  vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
533  uo(n,k)=urot
534  vo(n,k)=vrot
535  ELSE
536  ibo(k)=1
537  uo(n,k)=0.
538  vo(n,k)=0.
539  ENDIF
540  ENDDO
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  END SUBROUTINE interpolate_neighbor_budget_vector
550 
551 end module neighbor_budget_interp_mod
gdswzd_mod::gdswzd
Definition: gdswzd_mod.F90:37
neighbor_budget_interp_mod::interpolate_neighbor_budget
Definition: neighbor_budget_interp_mod.F90:17
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_budget_interp_mod
Interpolate scalar fields (neighbor).
Definition: neighbor_budget_interp_mod.F90:8
neighbor_budget_interp_mod::interpolate_neighbor_budget_vector
subroutine interpolate_neighbor_budget_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 (budget).
Definition: neighbor_budget_interp_mod.F90:354
neighbor_budget_interp_mod::interpolate_neighbor_budget_scalar
subroutine interpolate_neighbor_budget_scalar(IPOPT, grid_in, grid_out, MI, MO, KM, IBI, LI, GI, NO, RLAT, RLON, IBO, LO, GO, IRET)
Interpolate scalar fields (budget).
Definition: neighbor_budget_interp_mod.F90:110
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