NCEPLIBS-ip 4.0.0
budget_interp_mod.f90
Go to the documentation of this file.
1
5
10 use gdswzd_mod
11 use polfix_mod
12 use ip_grids_mod
15 implicit none
16
17 private
18 public :: interpolate_budget
19
21 module procedure interpolate_budget_scalar
22 module procedure interpolate_budget_vector
23 end interface interpolate_budget
24
25contains
26
89 SUBROUTINE interpolate_budget_scalar(IPOPT,grid_in,grid_out, &
90 MI,MO,KM,IBI,LI,GI, &
91 NO,RLAT,RLON,IBO,LO,GO,IRET)
92 class(ip_grid), intent(in) :: grid_in, grid_out
93 INTEGER, INTENT(IN ) :: IBI(KM), IPOPT(20)
94 INTEGER, INTENT(IN ) :: KM, MI, MO
95 INTEGER, INTENT( OUT) :: IBO(KM), IRET, NO
96 !
97 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
98 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
99 !
100 REAL, INTENT(IN ) :: GI(MI,KM)
101 REAL, INTENT(INOUT) :: RLAT(MO), RLON(MO)
102 REAL, INTENT( OUT) :: GO(MO,KM)
103 !
104 REAL, PARAMETER :: FILL=-9999.
105 !
106 INTEGER :: IJKGDS1, I1, J1, I2, J2, IB, JB
107 INTEGER :: IJKGDSA(20), IX, JX, IXS, JXS
108 INTEGER :: K, KXS, KXT, IGDTNUMO2
109 INTEGER :: LB, LSW, MP, MSPIRAL, MX
110 INTEGER :: N, NB, NB1, NB2, NB3, NB4, NV, NX
111 INTEGER :: N11(MO),N21(MO),N12(MO),N22(MO)
112 !
113 REAL :: GB, LAT(1), LON(1)
114 REAL :: PMP, RB2, RLOB(MO), RLAB(MO), WB
115 REAL :: W11(MO), W21(MO), W12(MO), W22(MO)
116 REAL :: WO(MO,KM), XF, YF, XI, YI, XX, YY
117 REAL :: XPTS(MO),YPTS(MO),XPTB(MO),YPTB(MO)
118 REAL :: XXX(1), YYY(1)
119 class(ip_grid), allocatable :: grid_out2
120 class(ip_grid_descriptor), allocatable :: grid_desc_out2
121
122 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
124 ! DO SUBSECTION OF GRID IF KGDSO(1) IS SUBTRACTED FROM 255.
125 iret=0
126
127 select type(grid_out)
128 type is(ip_station_points_grid)
129 allocate(grid_desc_out2, source = grid_out%descriptor)
130 grid_desc_out2%grid_num = 255 + grid_out%descriptor%grid_num
131 call init_grid(grid_out2, grid_desc_out2)
132
133 CALL gdswzd(grid_out2,-1,mo,fill,xpts,ypts,rlon,rlat,no)
134 IF(no.EQ.0) then
135 iret=3
136 end if
137 class default
138 allocate(grid_out2, source = grid_out)
139 CALL gdswzd(grid_out2, 0,mo,fill,xpts,ypts,rlon,rlat,no)
140 IF(no.EQ.0) iret=3
141 end select
142
143 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144 ! SET PARAMETERS
145 IF(ipopt(1).GT.16) iret=32
146 mspiral=max(ipopt(20),1)
147 nb1=ipopt(1)
148 IF(nb1.EQ.-1) nb1=2
149 IF(iret.EQ.0.AND.nb1.LT.0) iret=32
150 lsw=1
151 IF(ipopt(2).EQ.-2) lsw=2
152 IF(ipopt(1).EQ.-1.OR.ipopt(2).EQ.-1) lsw=0
153 IF(iret.EQ.0.AND.lsw.EQ.1.AND.nb1.GT.15) iret=32
154 mp=ipopt(3+ipopt(1))
155 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
156 IF(mp.LT.0.OR.mp.GT.100) iret=32
157 pmp=mp*0.01
158 IF(iret.EQ.0) THEN
159 nb2=2*nb1+1
160 rb2=1./nb2
161 nb3=nb2*nb2
162 nb4=nb3
163 IF(lsw.EQ.2) THEN
164 rb2=1./(nb1+1)
165 nb4=(nb1+1)**4
166 ELSEIF(lsw.EQ.1) THEN
167 nb4=ipopt(2)
168 DO ib=1,nb1
169 nb4=nb4+8*ib*ipopt(2+ib)
170 ENDDO
171 ENDIF
172 ELSE
173 nb3=0
174 nb4=1
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 DO nb=1,nb3
185 ! LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS
186 jb=(nb-1)/nb2-nb1
187 ib=nb-(jb+nb1)*nb2-nb1-1
188 lb=max(abs(ib),abs(jb))
189 wb=1
190 IF(lsw.EQ.2) THEN
191 wb=(nb1+1-abs(ib))*(nb1+1-abs(jb))
192 ELSEIF(lsw.EQ.1) THEN
193 wb=ipopt(2+lb)
194 ENDIF
195 IF(wb.NE.0) THEN
196 !$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC)
197 DO n=1,no
198 xptb(n)=xpts(n)+ib*rb2
199 yptb(n)=ypts(n)+jb*rb2
200 ENDDO
201 !$OMP END PARALLEL DO
202 CALL gdswzd(grid_out2, 1,no,fill,xptb,yptb,rlob,rlab,nv)
203 CALL gdswzd(grid_in,-1,no,fill,xptb,yptb,rlob,rlab,nv)
204 IF(iret.EQ.0.AND.nv.EQ.0.AND.lb.EQ.0) iret=2
205 !$OMP PARALLEL DO PRIVATE(N,XI,YI,I1,I2,J1,J2,XF,YF) SCHEDULE(STATIC)
206 DO n=1,no
207 xi=xptb(n)
208 yi=yptb(n)
209 IF(xi.NE.fill.AND.yi.NE.fill) THEN
210 i1=xi
211 i2=i1+1
212 j1=yi
213 j2=j1+1
214 xf=xi-i1
215 yf=yi-j1
216 n11(n)=grid_in%field_pos(i1, j1)
217 n21(n)=grid_in%field_pos(i2, j1)
218 n12(n)=grid_in%field_pos(i1, j2)
219 n22(n)=grid_in%field_pos(i2, j2)
220 IF(min(n11(n),n21(n),n12(n),n22(n)).GT.0) THEN
221 w11(n)=(1-xf)*(1-yf)
222 w21(n)=xf*(1-yf)
223 w12(n)=(1-xf)*yf
224 w22(n)=xf*yf
225 ELSE
226 n11(n)=0
227 n21(n)=0
228 n12(n)=0
229 n22(n)=0
230 ENDIF
231 ELSE
232 n11(n)=0
233 n21(n)=0
234 n12(n)=0
235 n22(n)=0
236 ENDIF
237 ENDDO
238 !$OMP END PARALLEL DO
239 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240 ! INTERPOLATE WITH OR WITHOUT BITMAPS
241 !$OMP PARALLEL DO PRIVATE(K,N,GB) SCHEDULE(STATIC)
242 DO k=1,km
243 DO n=1,no
244 IF(n11(n).GT.0) THEN
245 IF(ibi(k).EQ.0) THEN
246 gb=w11(n)*gi(n11(n),k)+w21(n)*gi(n21(n),k) &
247 +w12(n)*gi(n12(n),k)+w22(n)*gi(n22(n),k)
248 go(n,k)=go(n,k)+wb*gb
249 wo(n,k)=wo(n,k)+wb
250 ELSE
251 IF(li(n11(n),k)) THEN
252 go(n,k)=go(n,k)+wb*w11(n)*gi(n11(n),k)
253 wo(n,k)=wo(n,k)+wb*w11(n)
254 ENDIF
255 IF(li(n21(n),k)) THEN
256 go(n,k)=go(n,k)+wb*w21(n)*gi(n21(n),k)
257 wo(n,k)=wo(n,k)+wb*w21(n)
258 ENDIF
259 IF(li(n12(n),k)) THEN
260 go(n,k)=go(n,k)+wb*w12(n)*gi(n12(n),k)
261 wo(n,k)=wo(n,k)+wb*w12(n)
262 ENDIF
263 IF(li(n22(n),k)) THEN
264 go(n,k)=go(n,k)+wb*w22(n)*gi(n22(n),k)
265 wo(n,k)=wo(n,k)+wb*w22(n)
266 ENDIF
267 ENDIF
268 ENDIF
269 ENDDO
270 ENDDO
271 !$OMP END PARALLEL DO
272 ENDIF
273 ENDDO ! sub-grid points
274 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275 ! COMPUTE OUTPUT BITMAPS AND FIELDS
276 ! KM is often 1 .. do not do OMP PARALLEL DO here
277 km_loop : DO k=1,km
278 ibo(k)=ibi(k)
279 !$OMP PARALLEL DO PRIVATE(N,LAT,LON,XXX,YYY,NV,XX,YY,IXS,JXS,MX,KXS,KXT,IX,JX,NX) SCHEDULE(STATIC)
280 n_loop : DO n=1,no
281 lo(n,k)=wo(n,k).GE.pmp*nb4
282 IF(lo(n,k)) THEN
283 go(n,k)=go(n,k)/wo(n,k)
284 ELSEIF (mspiral.GT.1) THEN
285 lat(1)=rlat(n)
286 lon(1)=rlon(n)
287 CALL gdswzd(grid_in,-1,1,fill,xxx,yyy,lon,lat,nv)
288 xx=xxx(1)
289 yy=yyy(1)
290 IF(nv.EQ.1)THEN
291 i1=nint(xx)
292 j1=nint(yy)
293 ixs=sign(1.,xx-i1)
294 jxs=sign(1.,yy-j1)
295 spiral_loop : DO mx=2,mspiral**2
296 kxs=sqrt(4*mx-2.5)
297 kxt=mx-(kxs**2/4+1)
298 SELECT CASE(mod(kxs,4))
299 CASE(1)
300 ix=i1-ixs*(kxs/4-kxt)
301 jx=j1-jxs*kxs/4
302 CASE(2)
303 ix=i1+ixs*(1+kxs/4)
304 jx=j1-jxs*(kxs/4-kxt)
305 CASE(3)
306 ix=i1+ixs*(1+kxs/4-kxt)
307 jx=j1+jxs*(1+kxs/4)
308 CASE DEFAULT
309 ix=i1-ixs*kxs/4
310 jx=j1+jxs*(kxs/4-kxt)
311 END SELECT
312 nx=grid_in%field_pos(ix, jx)
313 IF(nx.GT.0.)THEN
314 IF(li(nx,k).OR.ibi(k).EQ.0) THEN
315 go(n,k)=gi(nx,k)
316 lo(n,k)=.true.
317 cycle n_loop
318 ENDIF
319 ENDIF
320 ENDDO spiral_loop
321 ibo(k)=1
322 go(n,k)=0.
323 ELSE
324 ibo(k)=1
325 go(n,k)=0.
326 ENDIF
327 ELSE ! no spiral search option
328 ibo(k)=1
329 go(n,k)=0.
330 ENDIF
331 ENDDO n_loop
332 !$OMP END PARALLEL DO
333 ENDDO km_loop
334
335 select type(grid_out2)
336 type is(ip_equid_cylind_grid)
337 CALL polfixs(no,mo,km,rlat,ibo,lo,go)
338 end select
339 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
340 END SUBROUTINE interpolate_budget_scalar
341
342
409 SUBROUTINE interpolate_budget_vector(IPOPT,grid_in,grid_out, &
410 MI,MO,KM,IBI,LI,UI,VI, &
411 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
412 class(ip_grid), intent(in) :: grid_in, grid_out
413 INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
414 INTEGER, INTENT(IN ) :: KM, MI, MO
415 INTEGER, INTENT( OUT) :: IRET, NO, IBO(KM)
416 !
417 LOGICAL*1, INTENT(IN ) :: LI(MI,KM)
418 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
419 !
420 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
421 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
422 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
423 REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
424 !
425 REAL, PARAMETER :: FILL=-9999.
426 !
427 INTEGER :: IGDTNUMO2
428 INTEGER :: I1,I2,J1,J2,IB,JB,LSW,MP
429 INTEGER :: K,LB,N,NB,NB1,NB2,NB3,NB4,NV
430 INTEGER :: N11(MO),N21(MO),N12(MO),N22(MO)
431 !
432 LOGICAL :: SAME_GRID
433 !
434 REAL :: CM11,SM11,CM12,SM12
435 REAL :: CM21,SM21,CM22,SM22
436 REAL :: PMP,RB2
437 REAL :: C11(MO),C21(MO),C12(MO),C22(MO)
438 REAL :: S11(MO),S21(MO),S12(MO),S22(MO)
439 REAL :: W11(MO),W21(MO),W12(MO),W22(MO)
440 REAL :: UB,VB,WB,UROT,VROT
441 REAL :: U11,V11,U21,V21,U12,V12,U22,V22
442 REAL :: WI1,WJ1,WI2,WJ2
443 REAL :: WO(MO,KM),XI,YI
444 REAL :: XPTS(MO),YPTS(MO)
445 REAL :: XPTB(MO),YPTB(MO),RLOB(MO),RLAB(MO)
446
447 class(ip_grid_descriptor), allocatable :: desc_out_subgrid
448 class(ip_grid), allocatable :: grid_out2
449
450 ! Save coeffecients between calls and only compute if grids have changed
451 INTEGER, SAVE :: MIX=-1
452 REAL, ALLOCATABLE, SAVE :: CROI(:),SROI(:)
453 REAL, ALLOCATABLE, SAVE :: XPTI(:),YPTI(:),RLOI(:),RLAI(:)
454
455 class(ip_grid), allocatable, save :: prev_grid_in
456
457 iret=0
458
459 ! Negative grid number means interpolate to subgrid
460 ! The type of the subgrid is calculated by 255 +
461 select type(grid_out)
462 type is(ip_station_points_grid)
463 allocate(desc_out_subgrid, source = grid_out%descriptor)
464 desc_out_subgrid%grid_num = 255 + grid_out%descriptor%grid_num
465
466 call init_grid(grid_out2, desc_out_subgrid)
467 CALL gdswzd(grid_out2,-1,mo,fill,xpts,ypts, &
468 rlon,rlat,no,crot,srot)
469 IF(no.EQ.0) iret=3
470 class default
471 allocate(grid_out2, source = grid_out)
472 CALL gdswzd(grid_out2, 0,mo,fill,xpts,ypts, &
473 rlon,rlat,no,crot,srot)
474 end select
475
476 if (.not. allocated(prev_grid_in)) then
477 allocate(prev_grid_in, source = grid_in)
478
479 same_grid = .false.
480 else
481 same_grid = grid_in == prev_grid_in
482
483 if (.not. same_grid) then
484 deallocate(prev_grid_in)
485 allocate(prev_grid_in, source = grid_in)
486 end if
487 end if
488
489 IF(.NOT.same_grid) THEN
490 IF(mix.NE.mi) THEN
491 IF(mix.GE.0) DEALLOCATE(xpti,ypti,rloi,rlai,croi,sroi)
492 ALLOCATE(xpti(mi),ypti(mi),rloi(mi),rlai(mi),croi(mi),sroi(mi))
493 mix=mi
494 ENDIF
495 CALL gdswzd(grid_in, 0,mi,fill,xpti,ypti, &
496 rloi,rlai,nv,croi,sroi)
497 ENDIF
498
499 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
500 ! SET PARAMETERS
501 nb1=ipopt(1)
502 IF(nb1.EQ.-1) nb1=2
503 IF(iret.EQ.0.AND.nb1.LT.0) iret=32
504 lsw=1
505 IF(ipopt(2).EQ.-2) lsw=2
506 IF(ipopt(1).EQ.-1.OR.ipopt(2).EQ.-1) lsw=0
507 IF(iret.EQ.0.AND.lsw.EQ.1.AND.nb1.GT.15) iret=32
508 mp=ipopt(3+ipopt(1))
509 IF(mp.EQ.-1.OR.mp.EQ.0) mp=50
510 IF(mp.LT.0.OR.mp.GT.100) iret=32
511 pmp=mp*0.01
512 IF(iret.EQ.0) THEN
513 nb2=2*nb1+1
514 rb2=1./nb2
515 nb3=nb2*nb2
516 nb4=nb3
517 IF(lsw.EQ.2) THEN
518 rb2=1./(nb1+1)
519 nb4=(nb1+1)**4
520 ELSEIF(lsw.EQ.1) THEN
521 nb4=ipopt(2)
522 DO ib=1,nb1
523 nb4=nb4+8*ib*ipopt(2+ib)
524 ENDDO
525 ENDIF
526 ELSE
527 nb3=0
528 nb4=1
529 ENDIF
530 DO k=1,km
531 DO n=1,no
532 uo(n,k)=0
533 vo(n,k)=0
534 wo(n,k)=0.
535 ENDDO
536 ENDDO
537 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
538 ! LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX
539 DO nb=1,nb3
540 ! LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS AND ROTATIONS
541 jb=(nb-1)/nb2-nb1
542 ib=nb-(jb+nb1)*nb2-nb1-1
543 lb=max(abs(ib),abs(jb))
544 wb=1
545 IF(ipopt(2).EQ.-2) THEN
546 wb=(nb1+1-abs(ib))*(nb1+1-abs(jb))
547 ELSEIF(ipopt(2).NE.-1) THEN
548 wb=ipopt(2+lb)
549 ENDIF
550 IF(wb.NE.0) THEN
551 !$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC)
552 DO n=1,no
553 xptb(n)=xpts(n)+ib*rb2
554 yptb(n)=ypts(n)+jb*rb2
555 ENDDO
556 !$OMP END PARALLEL DO
557 CALL gdswzd(grid_out2, 1,no,fill,xptb,yptb, &
558 rlob,rlab,nv)
559 CALL gdswzd(grid_in,-1,no,fill,xptb,yptb, &
560 rlob,rlab,nv)
561 IF(iret.EQ.0.AND.nv.EQ.0.AND.lb.EQ.0) iret=2
562 !$OMP PARALLEL DO PRIVATE(N,XI,YI,I1,I2,WI1,WI2,J1,J2,WJ1,WJ2,CM11,CM21,CM12,CM22,SM11,SM21,SM12,SM22) &
563 !$OMP SCHEDULE(STATIC)
564 DO n=1,no
565 xi=xptb(n)
566 yi=yptb(n)
567 IF(xi.NE.fill.AND.yi.NE.fill) THEN
568 i1=xi
569 i2=i1+1
570 wi2=xi-i1
571 wi1=1-wi2
572 j1=yi
573 j2=j1+1
574 wj2=yi-j1
575 wj1=1-wj2
576 n11(n) = grid_in%field_pos(i1,j1)
577 n21(n) = grid_in%field_pos(i2, j1)
578 n12(n) = grid_in%field_pos(i1, j2)
579 n22(n) = grid_in%field_pos(i2, j2)
580 IF(min(n11(n),n21(n),n12(n),n22(n)).GT.0) THEN
581 w11(n)=wi1*wj1
582 w21(n)=wi2*wj1
583 w12(n)=wi1*wj2
584 w22(n)=wi2*wj2
585 CALL movect(rlai(n11(n)),rloi(n11(n)),rlat(n),rlon(n),cm11,sm11)
586 CALL movect(rlai(n21(n)),rloi(n21(n)),rlat(n),rlon(n),cm21,sm21)
587 CALL movect(rlai(n12(n)),rloi(n12(n)),rlat(n),rlon(n),cm12,sm12)
588 CALL movect(rlai(n22(n)),rloi(n22(n)),rlat(n),rlon(n),cm22,sm22)
589 c11(n)=cm11*croi(n11(n))+sm11*sroi(n11(n))
590 s11(n)=sm11*croi(n11(n))-cm11*sroi(n11(n))
591 c21(n)=cm21*croi(n21(n))+sm21*sroi(n21(n))
592 s21(n)=sm21*croi(n21(n))-cm21*sroi(n21(n))
593 c12(n)=cm12*croi(n12(n))+sm12*sroi(n12(n))
594 s12(n)=sm12*croi(n12(n))-cm12*sroi(n12(n))
595 c22(n)=cm22*croi(n22(n))+sm22*sroi(n22(n))
596 s22(n)=sm22*croi(n22(n))-cm22*sroi(n22(n))
597 ELSE
598 n11(n)=0
599 n21(n)=0
600 n12(n)=0
601 n22(n)=0
602 ENDIF
603 ELSE
604 n11(n)=0
605 n21(n)=0
606 n12(n)=0
607 n22(n)=0
608 ENDIF
609 ENDDO
610 !$OMP END PARALLEL DO
611 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
612 ! INTERPOLATE WITH OR WITHOUT BITMAPS
613 ! KM IS OFTEN 1 .. DO NO PUT OMP PARALLEL DO HERE
614 DO k=1,km
615 !$OMP PARALLEL DO PRIVATE(N,U11,U12,U21,U22,UB,V11,V12,V21,V22,VB) SCHEDULE(STATIC)
616 DO n=1,no
617 IF(n11(n).GT.0) THEN
618 IF(ibi(k).EQ.0) THEN
619 u11=c11(n)*ui(n11(n),k)-s11(n)*vi(n11(n),k)
620 v11=s11(n)*ui(n11(n),k)+c11(n)*vi(n11(n),k)
621 u21=c21(n)*ui(n21(n),k)-s21(n)*vi(n21(n),k)
622 v21=s21(n)*ui(n21(n),k)+c21(n)*vi(n21(n),k)
623 u12=c12(n)*ui(n12(n),k)-s12(n)*vi(n12(n),k)
624 v12=s12(n)*ui(n12(n),k)+c12(n)*vi(n12(n),k)
625 u22=c22(n)*ui(n22(n),k)-s22(n)*vi(n22(n),k)
626 v22=s22(n)*ui(n22(n),k)+c22(n)*vi(n22(n),k)
627 ub=w11(n)*u11+w21(n)*u21+w12(n)*u12+w22(n)*u22
628 vb=w11(n)*v11+w21(n)*v21+w12(n)*v12+w22(n)*v22
629 uo(n,k)=uo(n,k)+wb*ub
630 vo(n,k)=vo(n,k)+wb*vb
631 wo(n,k)=wo(n,k)+wb
632 ELSE
633 IF(li(n11(n),k)) 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 uo(n,k)=uo(n,k)+wb*w11(n)*u11
637 vo(n,k)=vo(n,k)+wb*w11(n)*v11
638 wo(n,k)=wo(n,k)+wb*w11(n)
639 ENDIF
640 IF(li(n21(n),k)) THEN
641 u21=c21(n)*ui(n21(n),k)-s21(n)*vi(n21(n),k)
642 v21=s21(n)*ui(n21(n),k)+c21(n)*vi(n21(n),k)
643 uo(n,k)=uo(n,k)+wb*w21(n)*u21
644 vo(n,k)=vo(n,k)+wb*w21(n)*v21
645 wo(n,k)=wo(n,k)+wb*w21(n)
646 ENDIF
647 IF(li(n12(n),k)) THEN
648 u12=c12(n)*ui(n12(n),k)-s12(n)*vi(n12(n),k)
649 v12=s12(n)*ui(n12(n),k)+c12(n)*vi(n12(n),k)
650 uo(n,k)=uo(n,k)+wb*w12(n)*u12
651 vo(n,k)=vo(n,k)+wb*w12(n)*v12
652 wo(n,k)=wo(n,k)+wb*w12(n)
653 ENDIF
654 IF(li(n22(n),k)) THEN
655 u22=c22(n)*ui(n22(n),k)-s22(n)*vi(n22(n),k)
656 v22=s22(n)*ui(n22(n),k)+c22(n)*vi(n22(n),k)
657 uo(n,k)=uo(n,k)+wb*w22(n)*u22
658 vo(n,k)=vo(n,k)+wb*w22(n)*v22
659 wo(n,k)=wo(n,k)+wb*w22(n)
660 ENDIF
661 ENDIF
662 ENDIF
663 ENDDO
664 !$OMP END PARALLEL DO
665 ENDDO
666 ENDIF
667 ENDDO
668 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
669 ! COMPUTE OUTPUT BITMAPS AND FIELDS
670 ! KM is often 1, do not put OMP PARALLEL here
671 DO k=1,km
672 ibo(k)=ibi(k)
673 !$OMP PARALLEL DO PRIVATE(N,UROT,VROT) SCHEDULE(STATIC)
674 DO n=1,no
675 lo(n,k)=wo(n,k).GE.pmp*nb4
676 IF(lo(n,k)) THEN
677 uo(n,k)=uo(n,k)/wo(n,k)
678 vo(n,k)=vo(n,k)/wo(n,k)
679 urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
680 vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
681 uo(n,k)=urot
682 vo(n,k)=vrot
683 ELSE
684 ibo(k)=1
685 uo(n,k)=0.
686 vo(n,k)=0.
687 ENDIF
688 ENDDO
689 !$OMP END PARALLEL DO
690 ENDDO
691
692 select type(grid_out2)
693 type is(ip_equid_cylind_grid)
694 CALL polfixv(no,mo,km,rlat,rlon,ibo,lo,uo,vo)
695 end select
696 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
697 END SUBROUTINE interpolate_budget_vector
698
699end module budget_interp_mod
Budget interpolation routines for scalars and vectors.
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...
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)...
Driver module for gdswzd routines.
Definition: gdswzd_mod.f90:25
Uses derived type grid descriptor objects to abstract away the raw Grib-1 and Grib-2 grid definitions...
Routines for creating an ip_grid given a Grib descriptor.
Abstract descriptor object which represents a grib1 or grib2 descriptor.