NCEPLIBS-ip 5.3.0
•All Data Structures Namespaces Files Functions Variables Pages
budget_interp_mod.F90
Go to the documentation of this file.
1!> @file
2!! @brief Budget interpolation routines for scalars and vectors
3!! @author Mark Iredell, Kyle Gerheiser, Eric Engle
4!! @date July 2021
5
6!> @brief Budget interpolation routines for scalars and vectors.
7!!
8!! @author George Gayno, Mark Iredell, Kyle Gerheiser, Eric Engle
9!! @date July 2021
10module budget_interp_mod
11 use gdswzd_mod
12 use polfix_mod
13 use ip_grids_mod
14 use ip_grid_descriptor_mod
15 use ip_grid_factory_mod
16 implicit none
17
18 private
19 public :: interpolate_budget
20
21 interface interpolate_budget
22 module procedure interpolate_budget_scalar
23 module procedure interpolate_budget_vector
24 end interface interpolate_budget
25
26 ! Smallest positive real value (use for equality comparisons)
27 REAL :: TINYREAL=tiny(1.0)
28
29contains
30
31 !> Performs budget interpolation from any grid to any grid (or to
32 !> random station points) for scalar fields.
33 !>
34 !> The algorithm simply computes (weighted) averages of bilinearly
35 !> interpolated points arranged in a square box centered around each
36 !> output grid point and stretching nearly halfway to each of the
37 !> neighboring grid points.
38 !>
39 !> Options allow choices of number of points in each radius from the
40 !> center point (ipopt(1)) which defaults to 2 (if ipopt(1)=-1)
41 !> meaning that 25 points will be averaged; further options are the
42 !> respective weights for the radius points starting at the center
43 !> point (ipopt(2:2+ipopt(1)) which defaults to all 1 (if
44 !> ipopt(1)=-1 or ipopt(2)=-1).
45 !>
46 !> A special interpolation is done if ipopt(2)=-2. in this case,
47 !> the boxes stretch nearly all the way to each of the neighboring
48 !> grid points and the weights are the adjoint of the bilinear
49 !> interpolation weights. This case gives quasi-second-order budget
50 !> interpolation.
51 !>
52 !> Another option is the minimum percentage for mask, i.e. percent
53 !> valid input data required to make output data, (ipopt(3+ipopt(1))
54 !> which defaults to 50 (if -1).
55 !>
56 !> In cases where there is no or insufficient valid input data, the
57 !> user may choose to search for the nearest valid data. this is
58 !> invoked by setting ipopt(20) to the width of the search
59 !> square. The default is 1 (no search). Squares are searched for
60 !> valid data in a spiral pattern starting from the center. No
61 !> searching is done where the output grid is outside the input
62 !> grid.
63 !>
64 !> Only horizontal interpolation is performed.
65 !>
66 !> @param[in] ipopt Interpolation options
67 !> - ipopt(1) is number of radius points (defaults to 2 if ipopt(1)=-1).
68 !> - ipopt(2:2+ipopt(1)) are respective weights (defaults to all 1 if ipopt(1)=-1 or ipopt(2)=-1).
69 !> - ipopt(3+ipopt(1)) is minimum percentage for mask (defaults to 50 if ipopt(3+ipopt(1)=-1).
70 !> @param[in] grid_in Input grid
71 !> @param[in] grid_out Output grid
72 !> @param[in] mi Skip number between input grid fields if km>1 or
73 !> dimension of input grid fields if km=1.
74 !> @param[out] mo Skip number between output grid fields if km>1 or
75 !> dimension of output grid fields if km=1.
76 !> @param[in] km Number of fields to interpolate.
77 !> @param[in] ibi Input bitmap flags.
78 !> @param[in] li Input bitmaps (if some ibi(k)=1).
79 !> @param[in] gi Input fields to interpolate.
80 !> @param[in,out] no Number of output points (only if igdtnumo<0).
81 !> @param[in,out] rlat Output latitudes in degrees (if igdtnumo<0).
82 !> @param[in,out] rlon Output longitudes in degrees (if igdtnumo<0).
83 !> @param[out] ibo Output bitmap flags.
84 !> @param[out] lo Output bitmaps (always output).
85 !> @param[out] go Output fields interpolated.
86 !> @param[out] iret Return code.
87 !> - 0 Successful interpolation.
88 !> - 2 Unrecognized input grid or no grid overlap.
89 !> - 3 Unrecognized output grid.
90 !> - 32 Invalid budget method parameters.
91 !>
92 !> @author Marke Iredell, George Gayno, Kyle Gerheiser, Eric Engle
93 !> @date July 2021
94 SUBROUTINE interpolate_budget_scalar(IPOPT,grid_in,grid_out, &
95 MI,MO,KM,IBI,LI,GI, &
96 NO,RLAT,RLON,IBO,LO,GO,IRET)
97 class(ip_grid), intent(in) :: grid_in, grid_out
98 INTEGER, INTENT(IN ) :: IBI(KM), IPOPT(20)
99 INTEGER, INTENT(IN ) :: KM, MI, MO
100 INTEGER, INTENT( OUT) :: IBO(KM), IRET, NO
101 !
102 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
103 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
104 !
105 REAL, INTENT(IN ) :: GI(MI,KM)
106 REAL, INTENT(INOUT) :: RLAT(MO), RLON(MO)
107 REAL, INTENT( OUT) :: GO(MO,KM)
108 !
109 REAL, PARAMETER :: FILL=-9999.
110 !
111 INTEGER :: I1, J1, I2, J2, IB, JB
112 INTEGER :: IX, JX, IXS, JXS
113 INTEGER :: K, KXS, KXT
114 INTEGER :: LB, LSW, MP, MSPIRAL, MX
115 INTEGER :: N, NB, NB1, NB2, NB3, NB4, NV, NX
116 INTEGER :: N11(MO),N21(MO),N12(MO),N22(MO)
117 !
118 REAL :: GB, LAT(1), LON(1)
119 REAL :: PMP, RB2, RLOB(MO), RLAB(MO), WB
120 REAL :: W11(MO), W21(MO), W12(MO), W22(MO)
121 REAL :: WO(MO,KM), XF, YF, XI, YI, XX, YY
122 REAL :: XPTS(MO),YPTS(MO),XPTB(MO),YPTB(MO)
123 REAL :: XXX(1), YYY(1)
124
125 logical :: to_station_points
126
127 class(ip_grid), allocatable :: grid_out2
128
129 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
131 ! DO SUBSECTION OF GRID IF KGDSO(1) IS SUBTRACTED FROM 255.
132 iret=0
133
134 select type(grid_out)
135 type is(ip_station_points_grid)
136 to_station_points = .true.
137 allocate(grid_out2, source = grid_out)
138 CALL gdswzd(grid_out2, 0,mo,fill,xpts,ypts,rlon,rlat,no)
139 IF(no.EQ.0) iret=3
140 CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv)
141 IF(nv.EQ.0) iret=2
142 class default
143 to_station_points = .false.
144 allocate(grid_out2, source = grid_out)
145 CALL gdswzd(grid_out2, 0,mo,fill,xpts,ypts,rlon,rlat,no)
146 IF(no.EQ.0) iret=3
147 end select
148
149 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150 ! SET PARAMETERS
151 IF(ipopt(1).GT.16) iret=32
152 mspiral=max(ipopt(20),1)
153 nb1=ipopt(1)
154 IF(nb1.EQ.-1) nb1=2
155 IF(iret.EQ.0.AND.nb1.LT.0) iret=32
156 lsw=1
157 IF(ipopt(2).EQ.-2) lsw=2
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 rb2=1./nb2
167 nb3=nb2*nb2
168 nb4=nb3
169 IF(lsw.EQ.2) THEN
170 rb2=1./(nb1+1)
171 nb4=(nb1+1)**4
172 ELSEIF(lsw.EQ.1) THEN
173 nb4=ipopt(2)
174 DO ib=1,nb1
175 nb4=nb4+8*ib*ipopt(2+ib)
176 ENDDO
177 ENDIF
178 ELSE
179 nb3=0
180 nb4=1
181 ENDIF
182 DO k=1,km
183 DO n=1,no
184 go(n,k)=0.
185 wo(n,k)=0.
186 ENDDO
187 ENDDO
188 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189 ! LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX
190 DO nb=1,nb3
191 ! LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS
192 jb=(nb-1)/nb2-nb1
193 ib=nb-(jb+nb1)*nb2-nb1-1
194 lb=max(abs(ib),abs(jb))
195 wb=1
196 IF(lsw.EQ.2) THEN
197 wb=(nb1+1-abs(ib))*(nb1+1-abs(jb))
198 ELSEIF(lsw.EQ.1) THEN
199 wb=ipopt(2+lb)
200 ENDIF
201 IF(abs(wb).GT.tinyreal) THEN
202 !$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC)
203 DO n=1,no
204 xptb(n)=xpts(n)+ib*rb2
205 yptb(n)=ypts(n)+jb*rb2
206 ENDDO
207 !$OMP END PARALLEL DO
208 if(to_station_points)then
209 CALL gdswzd(grid_in, 1,no,fill,xptb,yptb,rlob,rlab,nv)
210 CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
211 else
212 CALL gdswzd(grid_out2, 1,no,fill,xptb,yptb,rlob,rlab,nv)
213 CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
214 endif
215 IF(iret.EQ.0.AND.nv.EQ.0.AND.lb.EQ.0) iret=2
216 !$OMP PARALLEL DO PRIVATE(N,XI,YI,I1,I2,J1,J2,XF,YF) SCHEDULE(STATIC)
217 DO n=1,no
218 xi=xptb(n)
219 yi=yptb(n)
220 IF(abs(xi-fill).GT.tinyreal.AND.abs(yi-fill).GT.tinyreal) THEN
221 i1=int(xi)
222 i2=i1+1
223 j1=int(yi)
224 j2=j1+1
225 xf=xi-i1
226 yf=yi-j1
227 n11(n)=grid_in%field_pos(i1, j1)
228 n21(n)=grid_in%field_pos(i2, j1)
229 n12(n)=grid_in%field_pos(i1, j2)
230 n22(n)=grid_in%field_pos(i2, j2)
231 IF(min(n11(n),n21(n),n12(n),n22(n)).GT.0) THEN
232 w11(n)=(1-xf)*(1-yf)
233 w21(n)=xf*(1-yf)
234 w12(n)=(1-xf)*yf
235 w22(n)=xf*yf
236 ELSE
237 n11(n)=0
238 n21(n)=0
239 n12(n)=0
240 n22(n)=0
241 ENDIF
242 ELSE
243 n11(n)=0
244 n21(n)=0
245 n12(n)=0
246 n22(n)=0
247 ENDIF
248 ENDDO
249 !$OMP END PARALLEL DO
250 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251 ! INTERPOLATE WITH OR WITHOUT BITMAPS
252 !$OMP PARALLEL DO PRIVATE(K,N,GB) SCHEDULE(STATIC)
253 DO k=1,km
254 DO n=1,no
255 IF(n11(n).GT.0) THEN
256 IF(ibi(k).EQ.0) THEN
257 gb=w11(n)*gi(n11(n),k)+w21(n)*gi(n21(n),k) &
258 +w12(n)*gi(n12(n),k)+w22(n)*gi(n22(n),k)
259 go(n,k)=go(n,k)+wb*gb
260 wo(n,k)=wo(n,k)+wb
261 ELSE
262 IF(li(n11(n),k)) THEN
263 go(n,k)=go(n,k)+wb*w11(n)*gi(n11(n),k)
264 wo(n,k)=wo(n,k)+wb*w11(n)
265 ENDIF
266 IF(li(n21(n),k)) THEN
267 go(n,k)=go(n,k)+wb*w21(n)*gi(n21(n),k)
268 wo(n,k)=wo(n,k)+wb*w21(n)
269 ENDIF
270 IF(li(n12(n),k)) THEN
271 go(n,k)=go(n,k)+wb*w12(n)*gi(n12(n),k)
272 wo(n,k)=wo(n,k)+wb*w12(n)
273 ENDIF
274 IF(li(n22(n),k)) THEN
275 go(n,k)=go(n,k)+wb*w22(n)*gi(n22(n),k)
276 wo(n,k)=wo(n,k)+wb*w22(n)
277 ENDIF
278 ENDIF
279 ENDIF
280 ENDDO
281 ENDDO
282 !$OMP END PARALLEL DO
283 ENDIF
284 ENDDO ! sub-grid points
285 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286 ! COMPUTE OUTPUT BITMAPS AND FIELDS
287 ! KM is often 1 .. do not do OMP PARALLEL DO here
288 km_loop : DO k=1,km
289 ibo(k)=ibi(k)
290 !$OMP PARALLEL DO PRIVATE(N,LAT,LON,XXX,YYY,NV,XX,YY,IXS,JXS,MX,KXS,KXT,IX,JX,NX) SCHEDULE(STATIC)
291 n_loop : DO n=1,no
292 lo(n,k)=wo(n,k).GE.pmp*nb4
293 IF(lo(n,k)) THEN
294 go(n,k)=go(n,k)/wo(n,k)
295 ELSEIF (mspiral.GT.1) THEN
296 lat(1)=rlat(n)
297 lon(1)=rlon(n)
298 CALL gdswzd(grid_in,-1,1,fill,xxx,yyy,lon,lat,nv)
299 xx=xxx(1)
300 yy=yyy(1)
301 IF(nv.EQ.1)THEN
302 i1=nint(xx)
303 j1=nint(yy)
304 ixs=int(sign(1.,xx-i1))
305 jxs=int(sign(1.,yy-j1))
306 spiral_loop : DO mx=2,mspiral**2
307 kxs=int(sqrt(4*mx-2.5))
308 kxt=mx-(kxs**2/4+1)
309 SELECT CASE(mod(kxs,4))
310 CASE(1)
311 ix=i1-ixs*(kxs/4-kxt)
312 jx=j1-jxs*kxs/4
313 CASE(2)
314 ix=i1+ixs*(1+kxs/4)
315 jx=j1-jxs*(kxs/4-kxt)
316 CASE(3)
317 ix=i1+ixs*(1+kxs/4-kxt)
318 jx=j1+jxs*(1+kxs/4)
319 CASE DEFAULT
320 ix=i1-ixs*kxs/4
321 jx=j1+jxs*(kxs/4-kxt)
322 END SELECT
323 nx=grid_in%field_pos(ix, jx)
324 IF(nx.GT.0.)THEN
325 IF(li(nx,k).OR.ibi(k).EQ.0) THEN
326 go(n,k)=gi(nx,k)
327 lo(n,k)=.true.
328 cycle n_loop
329 ENDIF
330 ENDIF
331 ENDDO spiral_loop
332 ibo(k)=1
333 go(n,k)=0.
334 ELSE
335 ibo(k)=1
336 go(n,k)=0.
337 ENDIF
338 ELSE ! no spiral search option
339 ibo(k)=1
340 go(n,k)=0.
341 ENDIF
342 ENDDO n_loop
343 !$OMP END PARALLEL DO
344 ENDDO km_loop
345
346 select type(grid_out2)
347 type is(ip_equid_cylind_grid)
348 CALL polfixs(no,mo,km,rlat,ibo,lo,go)
349 end select
350 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351 END SUBROUTINE interpolate_budget_scalar
352
353 !> This subprogram performs budget interpolation from any grid to
354 !> any grid (or to random station points) for vector fields.
355 !>
356 !> The algorithm simply computes (weighted) averages of bilinearly
357 !> interpolated points arranged in a square box centered around each
358 !> output grid point and stretching nearly halfway to each of the
359 !> neighboring grid points.
360 !>
361 !> Options allow choices of number of points in each radius from the
362 !> center point (ipopt(1)) which defaults to 2 (if ipopt(1)=-1)
363 !> meaning that 25 points will be averaged; further options are the
364 !> respective weights for the radius points starting at the center
365 !> point (ipopt(2:2+ipopt(1)) which defaults to all 1 (if
366 !> ipopt(1)=-1 or ipopt(2)=-1).
367 !>
368 !> A special interpolation is done if ipopt(2)=-2. in this case,
369 !> the boxes stretch nearly all the way to each of the neighboring
370 !> grid points and the weights are the adjoint of the bilinear
371 !> interpolation weights. This case gives quasi-second-order budget
372 !> interpolation.
373 !>
374 !> Another option is the minimum percentage for mask, i.e. percent
375 !> valid input data required to make output data, (ipopt(3+ipopt(1))
376 !> which defaults to 50 (if -1).
377 !>
378 !> In cases where there is no or insufficient valid input data, the
379 !> user may choose to search for the nearest valid data. this is
380 !> invoked by setting ipopt(20) to the width of the search
381 !> square. The default is 1 (no search). Squares are searched for
382 !> valid data in a spiral pattern starting from the center. No
383 !> searching is done where the output grid is outside the input
384 !> grid.
385 !>
386 !> Only horizontal interpolation is performed.
387 !>
388 !> @param[in] ipopt interpolation options ipopt(1) Number of radius
389 !> points (defaults to 2 if ipopt(1)=-1); ipopt(2:2+ipopt(1))
390 !> Respective weights (defaults to all 1 if ipopt(1)=-1 or
391 !> ipopt(2)=-1). ipopt(3+ipopt(1)) Minimum percentage for mask
392 !> (defaults to 50 if ipopt(3+ipopt(1)=-1)
393 !> @param[in] grid_in Input grid.
394 !> @param[in] grid_out Output grid.
395 !> @param[in] mi skip Number between input grid fields if km>1 or
396 !> dimension of input grid fields if km=1.
397 !> @param[out] mo skip Number between output grid fields if km>1 or
398 !> dimension of output grid fields if km=1.
399 !> @param[in] km Number of fields to interpolate.
400 !> @param[in] ibi Input bitmap flags.
401 !> @param[in] li Input bitmaps (if some ibi(k)=1).
402 !> @param[in] ui Input u-component fields to interpolate.
403 !> @param[in] vi Input v-component fields to interpolate.
404 !> @param[in,out] no Number of output points (only if igdtnumo<0)
405 !> @param[in,out] rlat Output latitudes in degrees (if igdtnumo<0)
406 !> @param[in,out] rlon Output longitudes in degrees (if igdtnumo<0)
407 !> @param[in,out] crot Vector rotation cosines. If interpolating
408 !> subgrid ugrid=crot * uearth - srot * vearth.
409 !> @param[in,out] srot Vector rotation sines. If interpolating
410 !> subgrid vgrid = srot * uearth + crot * vearth.
411 !> @param[out] ibo Output bitmap flags.
412 !> @param[out] lo Output bitmaps (always output).
413 !> @param[out] uo Output u-component fields interpolated.
414 !> @param[out] vo Output v-component fields interpolated.
415 !> @param[out] iret Return code.
416 !> - 0 Successful interpolation.
417 !> - 2 Unrecognized input grid or no grid overlap.
418 !> - 3 Unrecognized output grid.
419 !> - 32 Invalid budget method parameters.
420 !>
421 !> @author Marke Iredell, George Gayno, Kyle Gerheiser, Eric Engle
422 !> @date July 2021
423 SUBROUTINE interpolate_budget_vector(IPOPT,grid_in,grid_out, &
424 MI,MO,KM,IBI,LI,UI,VI, &
425 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
426 class(ip_grid), intent(in) :: grid_in, grid_out
427 INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
428 INTEGER, INTENT(IN ) :: KM, MI, MO
429 INTEGER, INTENT( OUT) :: IRET, NO, IBO(KM)
430 !
431 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
432 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
433 !
434 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
435 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
436 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
437 REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
438 !
439 REAL, PARAMETER :: FILL=-9999.
440 !
441 INTEGER :: I1,I2,J1,J2,IB,JB,LSW,MP
442 INTEGER :: K,LB,N,NB,NB1,NB2,NB3,NB4,NV
443 INTEGER :: N11(MO),N21(MO),N12(MO),N22(MO)
444 !
445 LOGICAL :: SAME_GRID
446 !
447 REAL :: CM11,SM11,CM12,SM12
448 REAL :: CM21,SM21,CM22,SM22
449 REAL :: PMP,RB2
450 REAL :: C11(MO),C21(MO),C12(MO),C22(MO)
451 REAL :: S11(MO),S21(MO),S12(MO),S22(MO)
452 REAL :: W11(MO),W21(MO),W12(MO),W22(MO)
453 REAL :: UB,VB,WB,UROT,VROT
454 REAL :: U11,V11,U21,V21,U12,V12,U22,V22
455 REAL :: WI1,WJ1,WI2,WJ2
456 REAL :: WO(MO,KM),XI,YI
457 REAL :: XPTS(MO),YPTS(MO)
458 REAL :: XPTB(MO),YPTB(MO),RLOB(MO),RLAB(MO)
459
460 logical :: to_station_points
461
462 class(ip_grid), allocatable :: grid_out2
463
464 ! Save coeffecients between calls and only compute if grids have changed
465 INTEGER, SAVE :: MIX=-1
466 REAL, ALLOCATABLE, SAVE :: CROI(:),SROI(:)
467 REAL, ALLOCATABLE, SAVE :: XPTI(:),YPTI(:),RLOI(:),RLAI(:)
468
469 class(ip_grid), allocatable, save :: prev_grid_in
470
471 iret=0
472
473 ! Negative grid number means interpolate to subgrid
474 ! The type of the subgrid is calculated by 255 +
475 select type(grid_out)
476 type is(ip_station_points_grid)
477 to_station_points = .true.
478 allocate(grid_out2, source = grid_out)
479 CALL gdswzd(grid_out2, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
480 IF(no.EQ.0) iret=3
481 CALL gdswzd(grid_in,-1,no,fill,xpts,ypts,rlon,rlat,nv,crot,srot)
482 IF(nv.EQ.0) iret=2
483 class default
484 to_station_points = .false.
485 allocate(grid_out2, source = grid_out)
486 CALL gdswzd(grid_out2, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
487 end select
488
489 if (.not. allocated(prev_grid_in)) then
490 allocate(prev_grid_in, source = grid_in)
491
492 same_grid = .false.
493 else
494 same_grid = grid_in == prev_grid_in
495
496 if (.not. same_grid) then
497 deallocate(prev_grid_in)
498 allocate(prev_grid_in, source = grid_in)
499 end if
500 end if
501
502 IF(.NOT.same_grid) THEN
503 IF(mix.NE.mi) THEN
504 IF(mix.GE.0) DEALLOCATE(xpti,ypti,rloi,rlai,croi,sroi)
505 ALLOCATE(xpti(mi),ypti(mi),rloi(mi),rlai(mi),croi(mi),sroi(mi))
506 mix=mi
507 ENDIF
508 CALL gdswzd(grid_in, 0,mi,fill,xpti,ypti,rloi,rlai,nv,croi,sroi)
509 ENDIF
510
511 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
512 ! SET PARAMETERS
513 nb1=ipopt(1)
514 IF(nb1.EQ.-1) nb1=2
515 IF(iret.EQ.0.AND.nb1.LT.0) iret=32
516 lsw=1
517 IF(ipopt(2).EQ.-2) lsw=2
518 IF(ipopt(1).EQ.-1.OR.ipopt(2).EQ.-1) lsw=0
519 IF(iret.EQ.0.AND.lsw.EQ.1.AND.nb1.GT.15) iret=32
520 mp=ipopt(3+ipopt(1))
521 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
522 IF(mp.LT.0.OR.mp.GT.100) iret=32
523 pmp=mp*0.01
524 IF(iret.EQ.0) THEN
525 nb2=2*nb1+1
526 rb2=1./nb2
527 nb3=nb2*nb2
528 nb4=nb3
529 IF(lsw.EQ.2) THEN
530 rb2=1./(nb1+1)
531 nb4=(nb1+1)**4
532 ELSEIF(lsw.EQ.1) THEN
533 nb4=ipopt(2)
534 DO ib=1,nb1
535 nb4=nb4+8*ib*ipopt(2+ib)
536 ENDDO
537 ENDIF
538 ELSE
539 nb3=0
540 nb4=1
541 ENDIF
542 DO k=1,km
543 DO n=1,no
544 uo(n,k)=0
545 vo(n,k)=0
546 wo(n,k)=0.
547 ENDDO
548 ENDDO
549 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
550 ! LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX
551 DO nb=1,nb3
552 ! LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS AND ROTATIONS
553 jb=(nb-1)/nb2-nb1
554 ib=nb-(jb+nb1)*nb2-nb1-1
555 lb=max(abs(ib),abs(jb))
556 wb=1
557 IF(ipopt(2).EQ.-2) THEN
558 wb=(nb1+1-abs(ib))*(nb1+1-abs(jb))
559 ELSEIF(ipopt(2).NE.-1) THEN
560 wb=ipopt(2+lb)
561 ENDIF
562 IF(abs(wb).GT.tinyreal) THEN
563 !$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC)
564 DO n=1,no
565 xptb(n)=xpts(n)+ib*rb2
566 yptb(n)=ypts(n)+jb*rb2
567 ENDDO
568 !$OMP END PARALLEL DO
569 if(to_station_points)then
570 CALL gdswzd(grid_in, 1,no,fill,xptb,yptb,rlob,rlab,nv)
571 CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
572 else
573 CALL gdswzd(grid_out2, 1,no,fill,xptb,yptb,rlob,rlab,nv)
574 CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
575 endif
576 IF(iret.EQ.0.AND.nv.EQ.0.AND.lb.EQ.0) iret=2
577 !$OMP PARALLEL DO PRIVATE(N,XI,YI,I1,I2,WI1,WI2,J1,J2,WJ1,WJ2,CM11,CM21,CM12,CM22,SM11,SM21,SM12,SM22) &
578 !$OMP SCHEDULE(STATIC)
579 DO n=1,no
580 xi=xptb(n)
581 yi=yptb(n)
582 IF(abs(xi-fill).GT.tinyreal.AND.abs(yi-fill).GT.tinyreal) THEN
583 i1=int(xi)
584 i2=i1+1
585 wi2=xi-i1
586 wi1=1-wi2
587 j1=int(yi)
588 j2=j1+1
589 wj2=yi-j1
590 wj1=1-wj2
591 n11(n) = grid_in%field_pos(i1,j1)
592 n21(n) = grid_in%field_pos(i2, j1)
593 n12(n) = grid_in%field_pos(i1, j2)
594 n22(n) = grid_in%field_pos(i2, j2)
595 IF(min(n11(n),n21(n),n12(n),n22(n)).GT.0) THEN
596 w11(n)=wi1*wj1
597 w21(n)=wi2*wj1
598 w12(n)=wi1*wj2
599 w22(n)=wi2*wj2
600 CALL movect(rlai(n11(n)),rloi(n11(n)),rlat(n),rlon(n),cm11,sm11)
601 CALL movect(rlai(n21(n)),rloi(n21(n)),rlat(n),rlon(n),cm21,sm21)
602 CALL movect(rlai(n12(n)),rloi(n12(n)),rlat(n),rlon(n),cm12,sm12)
603 CALL movect(rlai(n22(n)),rloi(n22(n)),rlat(n),rlon(n),cm22,sm22)
604 c11(n)=cm11*croi(n11(n))+sm11*sroi(n11(n))
605 s11(n)=sm11*croi(n11(n))-cm11*sroi(n11(n))
606 c21(n)=cm21*croi(n21(n))+sm21*sroi(n21(n))
607 s21(n)=sm21*croi(n21(n))-cm21*sroi(n21(n))
608 c12(n)=cm12*croi(n12(n))+sm12*sroi(n12(n))
609 s12(n)=sm12*croi(n12(n))-cm12*sroi(n12(n))
610 c22(n)=cm22*croi(n22(n))+sm22*sroi(n22(n))
611 s22(n)=sm22*croi(n22(n))-cm22*sroi(n22(n))
612 ELSE
613 n11(n)=0
614 n21(n)=0
615 n12(n)=0
616 n22(n)=0
617 ENDIF
618 ELSE
619 n11(n)=0
620 n21(n)=0
621 n12(n)=0
622 n22(n)=0
623 ENDIF
624 ENDDO
625 !$OMP END PARALLEL DO
626 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
627 ! INTERPOLATE WITH OR WITHOUT BITMAPS
628 ! KM IS OFTEN 1 .. DO NO PUT OMP PARALLEL DO HERE
629 DO k=1,km
630 !$OMP PARALLEL DO PRIVATE(N,U11,U12,U21,U22,UB,V11,V12,V21,V22,VB) SCHEDULE(STATIC)
631 DO n=1,no
632 IF(n11(n).GT.0) THEN
633 IF(ibi(k).EQ.0) THEN
634 u11=c11(n)*ui(n11(n),k)-s11(n)*vi(n11(n),k)
635 v11=s11(n)*ui(n11(n),k)+c11(n)*vi(n11(n),k)
636 u21=c21(n)*ui(n21(n),k)-s21(n)*vi(n21(n),k)
637 v21=s21(n)*ui(n21(n),k)+c21(n)*vi(n21(n),k)
638 u12=c12(n)*ui(n12(n),k)-s12(n)*vi(n12(n),k)
639 v12=s12(n)*ui(n12(n),k)+c12(n)*vi(n12(n),k)
640 u22=c22(n)*ui(n22(n),k)-s22(n)*vi(n22(n),k)
641 v22=s22(n)*ui(n22(n),k)+c22(n)*vi(n22(n),k)
642 ub=w11(n)*u11+w21(n)*u21+w12(n)*u12+w22(n)*u22
643 vb=w11(n)*v11+w21(n)*v21+w12(n)*v12+w22(n)*v22
644 uo(n,k)=uo(n,k)+wb*ub
645 vo(n,k)=vo(n,k)+wb*vb
646 wo(n,k)=wo(n,k)+wb
647 ELSE
648 IF(li(n11(n),k)) THEN
649 u11=c11(n)*ui(n11(n),k)-s11(n)*vi(n11(n),k)
650 v11=s11(n)*ui(n11(n),k)+c11(n)*vi(n11(n),k)
651 uo(n,k)=uo(n,k)+wb*w11(n)*u11
652 vo(n,k)=vo(n,k)+wb*w11(n)*v11
653 wo(n,k)=wo(n,k)+wb*w11(n)
654 ENDIF
655 IF(li(n21(n),k)) THEN
656 u21=c21(n)*ui(n21(n),k)-s21(n)*vi(n21(n),k)
657 v21=s21(n)*ui(n21(n),k)+c21(n)*vi(n21(n),k)
658 uo(n,k)=uo(n,k)+wb*w21(n)*u21
659 vo(n,k)=vo(n,k)+wb*w21(n)*v21
660 wo(n,k)=wo(n,k)+wb*w21(n)
661 ENDIF
662 IF(li(n12(n),k)) THEN
663 u12=c12(n)*ui(n12(n),k)-s12(n)*vi(n12(n),k)
664 v12=s12(n)*ui(n12(n),k)+c12(n)*vi(n12(n),k)
665 uo(n,k)=uo(n,k)+wb*w12(n)*u12
666 vo(n,k)=vo(n,k)+wb*w12(n)*v12
667 wo(n,k)=wo(n,k)+wb*w12(n)
668 ENDIF
669 IF(li(n22(n),k)) THEN
670 u22=c22(n)*ui(n22(n),k)-s22(n)*vi(n22(n),k)
671 v22=s22(n)*ui(n22(n),k)+c22(n)*vi(n22(n),k)
672 uo(n,k)=uo(n,k)+wb*w22(n)*u22
673 vo(n,k)=vo(n,k)+wb*w22(n)*v22
674 wo(n,k)=wo(n,k)+wb*w22(n)
675 ENDIF
676 ENDIF
677 ENDIF
678 ENDDO
679 !$OMP END PARALLEL DO
680 ENDDO
681 ENDIF
682 ENDDO
683 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
684 ! COMPUTE OUTPUT BITMAPS AND FIELDS
685 ! KM is often 1, do not put OMP PARALLEL here
686 DO k=1,km
687 ibo(k)=ibi(k)
688 !$OMP PARALLEL DO PRIVATE(N,UROT,VROT) SCHEDULE(STATIC)
689 DO n=1,no
690 lo(n,k)=wo(n,k).GE.pmp*nb4
691 IF(lo(n,k)) THEN
692 uo(n,k)=uo(n,k)/wo(n,k)
693 vo(n,k)=vo(n,k)/wo(n,k)
694 urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
695 vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
696 uo(n,k)=urot
697 vo(n,k)=vrot
698 ELSE
699 ibo(k)=1
700 uo(n,k)=0.
701 vo(n,k)=0.
702 ENDIF
703 ENDDO
704 !$OMP END PARALLEL DO
705 ENDDO
706
707 select type(grid_out2)
708 type is(ip_equid_cylind_grid)
709 CALL polfixv(no,mo,km,rlat,rlon,ibo,lo,uo,vo)
710 end select
711 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
712 END SUBROUTINE interpolate_budget_vector
713
714end module budget_interp_mod
budget_interp_mod::interpolate_budget
Definition budget_interp_mod.F90:21
gdswzd_mod::gdswzd
Definition gdswzd_mod.F90:37
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
budget_interp_mod
Budget interpolation routines for scalars and vectors.
Definition budget_interp_mod.F90:10
budget_interp_mod::interpolate_budget_scalar
subroutine interpolate_budget_scalar(ipopt, grid_in, grid_out, mi, mo, km, ibi, li, gi, no, rlat, rlon, ibo, lo, go, iret)
Performs budget interpolation from any grid to any grid (or to random station points) for scalar fiel...
Definition budget_interp_mod.F90:97
budget_interp_mod::interpolate_budget_vector
subroutine interpolate_budget_vector(ipopt, grid_in, grid_out, mi, mo, km, ibi, li, ui, vi, no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret)
This subprogram performs budget interpolation from any grid to any grid (or to random station points)...
Definition budget_interp_mod.F90:426
gdswzd_mod
Driver module for gdswzd routines.
Definition gdswzd_mod.F90:25
ip_grid_descriptor_mod
Users derived type grid descriptor objects to abstract away the raw GRIB1 and GRIB2 grid definitions.
Definition ip_grid_descriptor_mod.F90:15
ip_grid_factory_mod
Routines for creating an ip_grid given a Grib descriptor.
Definition ip_grid_factory_mod.F90:11
ip_grids_mod
Re-export the individual grids.
Definition ip_grids_mod.F90:7
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