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