NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
neighbor_budget_interp_mod.F90
Go to the documentation of this file.
1
4
8module 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
20 end interface interpolate_neighbor_budget
21
22 ! Smallest positive real value (use for equality comparisons)
23 REAL :: TINYREAL=tiny(1.0)
24
25contains
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
551end 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