NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
spectral_interp_mod.F90
Go to the documentation of this file.
1
4
8module spectral_interp_mod
9 use gdswzd_mod
10 use ip_grid_mod
11 use ip_grid_descriptor_mod
12 use ip_grid_factory_mod
13 use earth_radius_mod
14 use sp_mod
15 implicit none
16
17 private
18 public :: interpolate_spectral
19
20 interface interpolate_spectral
21 module procedure interpolate_spectral_scalar
22 module procedure interpolate_spectral_vector
23 end interface interpolate_spectral
24
25 interface polates4
26 module procedure polates4_grib1
27 module procedure polates4_grib2
28 end interface polates4
29
30 interface polatev4
31 module procedure polatev4_grib1
32 module procedure polatev4_grib2
33 end interface polatev4
34
35contains
36
61 subroutine interpolate_spectral_scalar(IPOPT,grid_in,grid_out, &
62 MI,MO,KM,IBI,GI, &
63 NO,RLAT,RLON,IBO,LO,GO,IRET)
64 INTEGER, INTENT(IN ) :: IPOPT(20)
65 class(ip_grid), intent(in) :: grid_in, grid_out
66 INTEGER, INTENT(IN ) :: MI, MO
67 INTEGER, INTENT(IN ) :: IBI(KM), KM
68 INTEGER, INTENT( OUT) :: IBO(KM), IRET, NO
69 !
70 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
71 !
72 REAL, INTENT(IN ) :: GI(MI,KM)
73 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
74 REAL, INTENT( OUT) :: GO(MO,KM)
75
76
77 select type(desc_in => grid_in%descriptor)
78 type is(grib1_descriptor)
79 select type(desc_out => grid_out%descriptor)
80 type is(grib1_descriptor)
81 CALL polates4(ipopt,desc_in%gds,desc_out%gds,mi,mo,km,ibi,gi,no,rlat,rlon,ibo,lo,go,iret)
82 end select
83
84 type is(grib2_descriptor)
85 select type(desc_out => grid_out%descriptor)
86 type is(grib2_descriptor)
87 CALL polates4(ipopt,desc_in%gdt_num,desc_in%gdt_tmpl,desc_in%gdt_len, &
88 desc_out%gdt_num,desc_out%gdt_tmpl,desc_out%gdt_len, &
89 mi,mo,km,ibi,gi,no,rlat,rlon,ibo,lo,go,iret)
90 end select
91 end select
92 end subroutine interpolate_spectral_scalar
93
122 subroutine interpolate_spectral_vector(IPOPT,grid_in,grid_out, &
123 MI,MO,KM,IBI,UI,VI, &
124 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
125 class(ip_grid), intent(in) :: grid_in, grid_out
126 INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
127 INTEGER, INTENT(IN ) :: KM, MI, MO
128 INTEGER, INTENT( OUT) :: IRET, IBO(KM), NO
129 !
130 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
131 !
132 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
133 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
134 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
135 REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
136
137
138 select type(desc_in => grid_in%descriptor)
139 type is(grib1_descriptor)
140 select type(desc_out => grid_out%descriptor)
141 type is(grib1_descriptor)
142 CALL polatev4_grib1(ipopt,desc_in%gds,desc_out%gds,mi,mo,km,ibi,ui,vi,&
143 no,rlat,rlon,crot,srot,ibo,lo,uo,vo,iret)
144 end select
145
146 type is(grib2_descriptor)
147 select type(desc_out => grid_out%descriptor)
148 type is(grib2_descriptor)
149 CALL polatev4(ipopt,desc_in%gdt_num,desc_in%gdt_tmpl,desc_in%gdt_len, &
150 desc_out%gdt_num,desc_out%gdt_tmpl,desc_out%gdt_len, &
151 mi,mo,km,ibi,ui,vi,&
152 no,rlat,rlon,crot,srot,ibo,lo,uo,vo,iret)
153 end select
154 end select
155
156
157 end subroutine interpolate_spectral_vector
158! @author Mark Iredell @date 96-04-10
255 SUBROUTINE polates4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
256 IGDTNUMO,IGDTMPLO,IGDTLENO, &
257 MI,MO,KM,IBI,GI, &
258 NO,RLAT,RLON,IBO,LO,GO,IRET)
259 INTEGER, INTENT(IN ) :: IGDTNUMI, IGDTLENI
260 INTEGER, INTENT(IN ) :: IGDTMPLI(IGDTLENI)
261 INTEGER, INTENT(IN ) :: IGDTNUMO, IGDTLENO
262 INTEGER, INTENT(IN ) :: IGDTMPLO(IGDTLENO)
263 INTEGER, INTENT(IN ) :: IPOPT(20)
264 INTEGER, INTENT(IN ) :: MI, MO
265 INTEGER, INTENT(IN ) :: IBI(KM), KM
266 INTEGER, INTENT( OUT) :: IBO(KM), IRET, NO
267 !
268 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
269 !
270 REAL, INTENT(IN ) :: GI(MI,KM)
271 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
272 REAL, INTENT( OUT) :: GO(MO,KM)
273 !
274 REAL, PARAMETER :: FILL=-9999.
275 REAL, PARAMETER :: PI=3.14159265358979
276 REAL, PARAMETER :: DPR=180./pi
277 !
278 INTEGER :: IDRTI, IDRTO, IG, JG, IM, JM
279 INTEGER :: IGO, JGO, IMO, JMO
280 INTEGER :: ISCAN, JSCAN, NSCAN
281 INTEGER :: ISCANO, JSCANO, NSCANO
282 INTEGER :: ISKIPI, JSKIPI, ISCALE
283 INTEGER :: IMAXI, JMAXI, ISPEC
284 INTEGER :: IP, IPRIME, IPROJ, IROMB, K
285 INTEGER :: MAXWV, N, NI, NJ, NPS
286 !
287 REAL :: DE, DR, DY
288 REAL :: DLAT, DLON, DLATO, DLONO
289 REAL :: GO2(MO,KM), H, HI, HJ
290 REAL :: ORIENT, SLAT, RERTH, E2
291 REAL :: RLAT1, RLON1, RLAT2, RLON2, RLATI
292 REAL :: XMESH, XP, YP
293 REAL :: XPTS(MO), YPTS(MO)
294
295 type(grib2_descriptor) :: desc_in, desc_out
296 class(ip_grid), allocatable :: grid_in, grid_out
297
298 desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
299 desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
300
301 call init_grid(grid_in, desc_in)
302 call init_grid(grid_out, desc_out)
303 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
304 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
305 iret=0
306 IF(igdtnumo.GE.0) THEN
307 !CALL GDSWZD(IGDTNUMO,IGDTMPLO,IGDTLENO, 0,MO,FILL,XPTS,YPTS,RLON,RLAT,NO)
308 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
309 IF(no.EQ.0) iret=3
310 ENDIF
311 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
312 ! AFFIRM APPROPRIATE INPUT GRID
313 ! LAT/LON OR GAUSSIAN
314 ! NO BITMAPS
315 ! FULL ZONAL COVERAGE
316 ! FULL MERIDIONAL COVERAGE
317 idrti=igdtnumi
318 IF(idrti==40) idrti=4
319 IF(idrti==0.OR.idrti==4)THEN
320 im=igdtmpli(8)
321 jm=igdtmpli(9)
322 iscale=igdtmpli(10)*igdtmpli(11)
323 IF(iscale==0) iscale=10**6
324 rlon1=float(igdtmpli(13))/float(iscale)
325 rlon2=float(igdtmpli(16))/float(iscale)
326 iscan=mod(igdtmpli(19)/128,2)
327 jscan=mod(igdtmpli(19)/64,2)
328 nscan=mod(igdtmpli(19)/32,2)
329 ELSE
330 iret=41
331 ENDIF
332 DO k=1,km
333 IF(ibi(k).NE.0) iret=41
334 ENDDO
335 IF(iret.EQ.0) THEN
336 IF(iscan.EQ.0) THEN
337 dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
338 ELSE
339 dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
340 ENDIF
341 ig=nint(360/abs(dlon))
342 iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
343 imaxi=ig
344 jmaxi=jm
345 IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
346 ENDIF
347 IF(iret.EQ.0.AND.idrti.EQ.0) THEN
348 iscale=igdtmpli(10)*igdtmpli(11)
349 IF(iscale==0) iscale=10**6
350 rlat1=float(igdtmpli(12))/float(iscale)
351 rlat2=float(igdtmpli(15))/float(iscale)
352 dlat=(rlat2-rlat1)/(jm-1)
353 jg=nint(180/abs(dlat))
354 IF(jm.EQ.jg) idrti=256
355 IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
356 ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
357 jg=igdtmpli(18)*2
358 IF(jm.NE.jg) iret=41
359 ENDIF
360 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
361 ! SET PARAMETERS
362 IF(iret.EQ.0) THEN
363 iromb=ipopt(1)
364 maxwv=ipopt(2)
365 IF(maxwv.EQ.-1) THEN
366 IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
367 IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
368 IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
369 IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
370 IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
371 IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
372 ENDIF
373 IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
374 ENDIF
375 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
376 ! INTERPOLATE
377 IF(iret.EQ.0) THEN
378 IF(nscan.EQ.0) THEN
379 iskipi=1
380 jskipi=im
381 ELSE
382 iskipi=jm
383 jskipi=1
384 ENDIF
385 IF(iscan.EQ.1) iskipi=-iskipi
386 IF(jscan.EQ.0) jskipi=-jskipi
387 ispec=0
388 ! SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
389 IF((igdtnumo.EQ.0.OR.igdtnumo.EQ.40).AND. &
390 mod(igdtmplo(8),2).EQ.0.AND.igdtmplo(13).EQ.0.AND.igdtmplo(19).EQ.0) THEN
391 idrto=igdtnumo
392 IF(idrto==40)idrto=4
393 imo=igdtmplo(8)
394 jmo=igdtmplo(9)
395 iscale=igdtmplo(10)*igdtmplo(11)
396 IF(iscale==0) iscale=10**6
397 rlon2=float(igdtmplo(16))/float(iscale)
398 dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
399 igo=nint(360/abs(dlono))
400 IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
401 rlat1=float(igdtmplo(12))/float(iscale)
402 rlat2=float(igdtmplo(15))/float(iscale)
403 dlat=(rlat2-rlat1)/(jmo-1)
404 jgo=nint(180/abs(dlat))
405 IF(jmo.EQ.jgo) idrto=256
406 IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
407 ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
408 jgo=igdtmplo(18)*2
409 IF(jmo.EQ.jgo) ispec=1
410 ENDIF
411 IF(ispec.EQ.1) THEN
412 CALL sptrun(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
413 km,iprime,iskipi,jskipi,mi,0,0,mo,0,gi,go)
414 ENDIF
415 ! SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
416 ELSEIF(igdtnumo.EQ.20.AND. &
417 igdtmplo(8).EQ.igdtmplo(9).AND.mod(igdtmplo(8),2).EQ.1.AND. &
418 igdtmplo(15).EQ.igdtmplo(16).AND.igdtmplo(18).EQ.64) THEN
419 nps=igdtmplo(8)
420 rlat1=float(igdtmplo(10))*1.e-6
421 rlon1=float(igdtmplo(11))*1.e-6
422 orient=float(igdtmplo(14))*1.e-6
423 xmesh=float(igdtmplo(15))*1.e-3
424 iproj=mod(igdtmplo(17)/128,2)
425 ip=(nps+1)/2
426 h=(-1.)**iproj
427 slat=float(abs(igdtmplo(13)))*1.e-6
428 CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
429 de=(1.+sin(slat/dpr))*rerth
430 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
431 xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
432 yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
433 IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
434 IF(iproj.EQ.0) THEN
435 CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
436 iprime,iskipi,jskipi,mi,mo,0,0,0, &
437 slat,xmesh,orient,gi,go,go2)
438 ELSE
439 CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
440 iprime,iskipi,jskipi,mi,mo,0,0,0, &
441 slat,xmesh,orient,gi,go2,go)
442 ENDIF
443 ispec=1
444 ENDIF
445 ! SPECIAL CASE OF MERCATOR GRID
446 ELSEIF(igdtnumo.EQ.10) THEN
447 ni=igdtmplo(8)
448 nj=igdtmplo(9)
449 rlat1=float(igdtmplo(10))*1.0e-6
450 rlon1=float(igdtmplo(11))*1.0e-6
451 rlon2=float(igdtmplo(15))*1.0e-6
452 rlati=float(igdtmplo(13))*1.0e-6
453 iscano=mod(igdtmplo(16)/128,2)
454 jscano=mod(igdtmplo(16)/64,2)
455 nscano=mod(igdtmplo(16)/32,2)
456 dy=float(igdtmplo(19))*1.0e-3
457 hi=(-1.)**iscano
458 hj=(-1.)**(1-jscano)
459 CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
460 dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
461 dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
462 IF(nscano.EQ.0) THEN
463 CALL sptrunm(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
464 iprime,iskipi,jskipi,mi,mo,0,0,0, &
465 rlat1,rlon1,dlato,dlono,gi,go)
466 ispec=1
467 ENDIF
468 ENDIF
469 ! GENERAL SLOW CASE
470 IF(ispec.EQ.0) THEN
471 CALL sptrung(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
472 iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon,gi,go)
473 ENDIF
474 DO k=1,km
475 ibo(k)=0
476 DO n=1,no
477 lo(n,k)=.true.
478 ENDDO
479 ENDDO
480 ELSE
481 DO k=1,km
482 ibo(k)=1
483 DO n=1,no
484 lo(n,k)=.false.
485 go(n,k)=0.
486 ENDDO
487 ENDDO
488 ENDIF
489 END SUBROUTINE polates4_grib2
490
560 suBROUTINE polates4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,GI, &
561 NO,RLAT,RLON,IBO,LO,GO,IRET)
562 INTEGER, INTENT(IN ) :: IPOPT(20), KGDSI(200)
563 INTEGER, INTENT(IN ) :: KGDSO(200), MI, MO
564 INTEGER, INTENT(IN ) :: IBI(KM), KM
565 INTEGER, INTENT( OUT) :: IBO(KM), IRET
566 !
567 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
568 !
569 REAL, INTENT(IN ) :: GI(MI,KM)
570 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
571 REAL, INTENT( OUT) :: GO(MO,KM)
572 !
573 REAL, PARAMETER :: FILL=-9999.
574 REAL, PARAMETER :: RERTH=6.3712e6
575 REAL, PARAMETER :: PI=3.14159265358979
576 REAL, PARAMETER :: DPR=180./pi
577 !
578 INTEGER :: IDRTI, IDRTO, IG, JG, IM, JM
579 INTEGER :: IGO, JGO, IMO, JMO
580 INTEGER :: ISCAN, JSCAN, NSCAN
581 INTEGER :: ISCANO, JSCANO, NSCANO
582 INTEGER :: ISKIPI, JSKIPI
583 INTEGER :: IMAXI, JMAXI, ISPEC
584 INTEGER :: IP, IPRIME, IPROJ, IROMB, K
585 INTEGER :: MAXWV, N, NI, NJ, NPS, NO
586 !
587 REAL :: DE, DR, DY
588 REAL :: DLAT, DLON, DLATO, DLONO
589 REAL :: GO2(MO,KM), H, HI, HJ
590 REAL :: ORIENT
591 REAL :: RLAT1, RLON1, RLAT2, RLON2, RLATI
592 REAL :: XMESH, XP, YP
593 REAL :: XPTS(MO), YPTS(MO)
594
595 type(grib1_descriptor) :: desc_in, desc_out
596 class(ip_grid), allocatable :: grid_in, grid_out
597
598 desc_in = init_descriptor(kgdsi)
599 desc_out = init_descriptor(kgdso)
600
601 call init_grid(grid_in, desc_in)
602 call init_grid(grid_out, desc_out)
603
604 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
605 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
606 iret=0
607 IF(kgdso(1).GE.0) THEN
608 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
609 IF(no.EQ.0) iret=3
610 ENDIF
611 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
612 ! AFFIRM APPROPRIATE INPUT GRID
613 ! LAT/LON OR GAUSSIAN
614 ! NO BITMAPS
615 ! FULL ZONAL COVERAGE
616 ! FULL MERIDIONAL COVERAGE
617 idrti=kgdsi(1)
618 im=kgdsi(2)
619 jm=kgdsi(3)
620 rlon1=kgdsi(5)*1.e-3
621 rlon2=kgdsi(8)*1.e-3
622 iscan=mod(kgdsi(11)/128,2)
623 jscan=mod(kgdsi(11)/64,2)
624 nscan=mod(kgdsi(11)/32,2)
625 IF(idrti.NE.0.AND.idrti.NE.4) iret=41
626 DO k=1,km
627 IF(ibi(k).NE.0) iret=41
628 ENDDO
629 IF(iret.EQ.0) THEN
630 IF(iscan.EQ.0) THEN
631 dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
632 ELSE
633 dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
634 ENDIF
635 ig=nint(360/abs(dlon))
636 iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
637 imaxi=ig
638 jmaxi=jm
639 IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
640 ENDIF
641 IF(iret.EQ.0.AND.idrti.EQ.0) THEN
642 rlat1=kgdsi(4)*1.e-3
643 rlat2=kgdsi(7)*1.e-3
644 dlat=(rlat2-rlat1)/(jm-1)
645 jg=nint(180/abs(dlat))
646 IF(jm.EQ.jg) idrti=256
647 IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
648 ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
649 jg=kgdsi(10)*2
650 IF(jm.NE.jg) iret=41
651 ENDIF
652 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
653 ! SET PARAMETERS
654 IF(iret.EQ.0) THEN
655 iromb=ipopt(1)
656 maxwv=ipopt(2)
657 IF(maxwv.EQ.-1) THEN
658 IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
659 IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
660 IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
661 IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
662 IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
663 IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
664 ENDIF
665 IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
666 ENDIF
667 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
668 ! INTERPOLATE
669 IF(iret.EQ.0) THEN
670 IF(nscan.EQ.0) THEN
671 iskipi=1
672 jskipi=im
673 ELSE
674 iskipi=jm
675 jskipi=1
676 ENDIF
677 IF(iscan.EQ.1) iskipi=-iskipi
678 IF(jscan.EQ.0) jskipi=-jskipi
679 ispec=0
680 ! SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
681 IF((kgdso(1).EQ.0.OR.kgdso(1).EQ.4).AND. &
682 mod(kgdso(2),2).EQ.0.AND.kgdso(5).EQ.0.AND.kgdso(11).EQ.0) THEN
683 idrto=kgdso(1)
684 imo=kgdso(2)
685 jmo=kgdso(3)
686 rlon2=kgdso(8)*1.e-3
687 dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
688 igo=nint(360/abs(dlono))
689 IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
690 rlat1=kgdso(4)*1.e-3
691 rlat2=kgdso(7)*1.e-3
692 dlat=(rlat2-rlat1)/(jmo-1)
693 jgo=nint(180/abs(dlat))
694 IF(jmo.EQ.jgo) idrto=256
695 IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
696 ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
697 jgo=kgdso(10)*2
698 IF(jmo.EQ.jgo) ispec=1
699 ENDIF
700 IF(ispec.EQ.1) THEN
701 CALL sptrun(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
702 km,iprime,iskipi,jskipi,mi,0,0,mo,0,gi,go)
703 ENDIF
704 ! SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
705 ELSEIF(kgdso(1).EQ.5.AND. &
706 kgdso(2).EQ.kgdso(3).AND.mod(kgdso(2),2).EQ.1.AND. &
707 kgdso(8).EQ.kgdso(9).AND.kgdso(11).EQ.64) THEN
708 nps=kgdso(2)
709 rlat1=kgdso(4)*1.e-3
710 rlon1=kgdso(5)*1.e-3
711 orient=kgdso(7)*1.e-3
712 xmesh=kgdso(8)
713 iproj=mod(kgdso(10)/128,2)
714 ip=(nps+1)/2
715 h=(-1.)**iproj
716 de=(1.+sin(60./dpr))*rerth
717 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
718 xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
719 yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
720 IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
721 IF(iproj.EQ.0) THEN
722 CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
723 iprime,iskipi,jskipi,mi,mo,0,0,0, &
724 60.,xmesh,orient,gi,go,go2)
725 ELSE
726 CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
727 iprime,iskipi,jskipi,mi,mo,0,0,0, &
728 60.,xmesh,orient,gi,go2,go)
729 ENDIF
730 ispec=1
731 ENDIF
732 ! SPECIAL CASE OF MERCATOR GRID
733 ELSEIF(kgdso(1).EQ.1) THEN
734 ni=kgdso(2)
735 nj=kgdso(3)
736 rlat1=kgdso(4)*1.e-3
737 rlon1=kgdso(5)*1.e-3
738 rlon2=kgdso(8)*1.e-3
739 rlati=kgdso(9)*1.e-3
740 iscano=mod(kgdso(11)/128,2)
741 jscano=mod(kgdso(11)/64,2)
742 nscano=mod(kgdso(11)/32,2)
743 dy=kgdso(13)
744 hi=(-1.)**iscano
745 hj=(-1.)**(1-jscano)
746 dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
747 dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
748 IF(nscano.EQ.0) THEN
749 CALL sptrunm(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
750 iprime,iskipi,jskipi,mi,mo,0,0,0, &
751 rlat1,rlon1,dlato,dlono,gi,go)
752 ispec=1
753 ENDIF
754 ENDIF
755 ! GENERAL SLOW CASE
756 IF(ispec.EQ.0) THEN
757 CALL sptrung(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
758 iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon,gi,go)
759 ENDIF
760 DO k=1,km
761 ibo(k)=0
762 DO n=1,no
763 lo(n,k)=.true.
764 ENDDO
765 ENDDO
766 ELSE
767 DO k=1,km
768 ibo(k)=1
769 DO n=1,no
770 lo(n,k)=.false.
771 go(n,k)=0.
772 ENDDO
773 ENDDO
774 ENDIF
775 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
776 END SUBROUTINE polates4_grib1
777
778
888 SUBROUTINE polatev4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
889 IGDTNUMO,IGDTMPLO,IGDTLENO, &
890 MI,MO,KM,IBI,UI,VI, &
891 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
892 INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
893 INTEGER, INTENT(IN ) :: KM, MI, MO
894 INTEGER, INTENT( OUT) :: IRET, IBO(KM), NO
895 INTEGER, INTENT(IN ) :: IGDTNUMI, IGDTLENI
896 INTEGER, INTENT(IN ) :: IGDTMPLI(IGDTLENI)
897 INTEGER, INTENT(IN ) :: IGDTNUMO, IGDTLENO
898 INTEGER, INTENT(IN ) :: IGDTMPLO(IGDTLENO)
899 !
900 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
901 !
902 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
903 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
904 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
905 REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
906 !
907 REAL, PARAMETER :: FILL=-9999.
908 REAL, PARAMETER :: PI=3.14159265358979
909 REAL, PARAMETER :: DPR=180./pi
910 !
911 INTEGER :: IDRTO, IROMB, ISKIPI, ISPEC
912 INTEGER :: IDRTI, IMAXI, JMAXI, IM, JM
913 INTEGER :: IPRIME, IG, IMO, JMO, IGO, JGO
914 INTEGER :: ISCAN, JSCAN, NSCAN
915 INTEGER :: ISCANO, JSCANO, NSCANO
916 INTEGER :: ISCALE, IP, IPROJ, JSKIPI, JG
917 INTEGER :: K, MAXWV, N, NI, NJ, NPS
918 !
919 REAL :: DLAT, DLON, DLATO, DLONO, DE, DR, DY
920 REAL :: E2, H, HI, HJ, DUMM(1)
921 REAL :: ORIENT, RERTH, SLAT
922 REAL :: RLAT1, RLON1, RLAT2, RLON2, RLATI
923 REAL :: UROT, VROT, UO2(MO,KM),VO2(MO,KM)
924 REAL :: XMESH, XP, YP, XPTS(MO),YPTS(MO)
925
926 type(grib2_descriptor) :: desc_in, desc_out
927 class(ip_grid), allocatable :: grid_in, grid_out
928
929 desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
930 desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
931
932 call init_grid(grid_in, desc_in)
933 call init_grid(grid_out, desc_out)
934
935 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
936 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
937 iret=0
938 IF(igdtnumo.GE.0) THEN
939 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts, &
940 rlon,rlat,no,crot,srot)
941 IF(no.EQ.0) iret=3
942 ENDIF
943 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
944 ! AFFIRM APPROPRIATE INPUT GRID
945 ! LAT/LON OR GAUSSIAN
946 ! NO BITMAPS
947 ! FULL ZONAL COVERAGE
948 ! FULL MERIDIONAL COVERAGE
949 idrti=igdtnumi
950 IF(idrti==40) idrti=4
951 IF(idrti==0.OR.idrti==4)THEN
952 im=igdtmpli(8)
953 jm=igdtmpli(9)
954 iscale=igdtmpli(10)*igdtmpli(11)
955 IF(iscale==0) iscale=10**6
956 rlon1=float(igdtmpli(13))/float(iscale)
957 rlon2=float(igdtmpli(16))/float(iscale)
958 iscan=mod(igdtmpli(19)/128,2)
959 jscan=mod(igdtmpli(19)/64,2)
960 nscan=mod(igdtmpli(19)/32,2)
961 ELSE
962 iret=41
963 ENDIF
964 DO k=1,km
965 IF(ibi(k).NE.0) iret=41
966 ENDDO
967 IF(iret.EQ.0) THEN
968 IF(iscan.EQ.0) THEN
969 dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
970 ELSE
971 dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
972 ENDIF
973 ig=nint(360/abs(dlon))
974 iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
975 imaxi=ig
976 jmaxi=jm
977 IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
978 ENDIF
979 IF(iret.EQ.0.AND.idrti.EQ.0) THEN
980 iscale=igdtmpli(10)*igdtmpli(11)
981 IF(iscale==0) iscale=10**6
982 rlat1=float(igdtmpli(12))/float(iscale)
983 rlat2=float(igdtmpli(15))/float(iscale)
984 dlat=(rlat2-rlat1)/(jm-1)
985 jg=nint(180/abs(dlat))
986 IF(jm.EQ.jg) idrti=256
987 IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
988 ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
989 jg=igdtmpli(18)*2
990 IF(jm.NE.jg) iret=41
991 ENDIF
992 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
993 ! SET PARAMETERS
994 IF(iret.EQ.0) THEN
995 iromb=ipopt(1)
996 maxwv=ipopt(2)
997 IF(maxwv.EQ.-1) THEN
998 IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
999 IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
1000 IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
1001 IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
1002 IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
1003 IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
1004 ENDIF
1005 IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
1006 ENDIF
1007 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1008 ! INTERPOLATE
1009 IF(iret.EQ.0) THEN
1010 IF(nscan.EQ.0) THEN
1011 iskipi=1
1012 jskipi=im
1013 ELSE
1014 iskipi=jm
1015 jskipi=1
1016 ENDIF
1017 IF(iscan.EQ.1) iskipi=-iskipi
1018 IF(jscan.EQ.0) jskipi=-jskipi
1019 ispec=0
1020 ! SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
1021 IF((igdtnumo.EQ.0.OR.igdtnumo.EQ.40).AND. &
1022 mod(igdtmplo(8),2).EQ.0.AND.igdtmplo(13).EQ.0.AND. &
1023 igdtmplo(19).EQ.0) THEN
1024 idrto=igdtnumo
1025 IF(idrto==40)idrto=4
1026 imo=igdtmplo(8)
1027 jmo=igdtmplo(9)
1028 iscale=igdtmplo(10)*igdtmplo(11)
1029 IF(iscale==0) iscale=10**6
1030 rlon2=float(igdtmplo(16))/float(iscale)
1031 dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
1032 igo=nint(360/abs(dlono))
1033 IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
1034 rlat1=float(igdtmplo(12))/float(iscale)
1035 rlat2=float(igdtmplo(15))/float(iscale)
1036 dlat=(rlat2-rlat1)/(jmo-1)
1037 jgo=nint(180/abs(dlat))
1038 IF(jmo.EQ.jgo) idrto=256
1039 IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
1040 ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
1041 jgo=igdtmplo(18)*2
1042 IF(jmo.EQ.jgo) ispec=1
1043 ENDIF
1044 IF(ispec.EQ.1) THEN
1045 CALL sptrunv(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
1046 km,iprime,iskipi,jskipi,mi,0,0,mo,0,ui,vi, &
1047 .true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
1048 ENDIF
1049 ! SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
1050 ELSEIF(igdtnumo.EQ.20.AND. &
1051 igdtmplo(8).EQ.igdtmplo(9).AND.mod(igdtmplo(8),2).EQ.1.AND. &
1052 igdtmplo(15).EQ.igdtmplo(16).AND.igdtmplo(18).EQ.64.AND. &
1053 mod(igdtmplo(12)/8,2).EQ.1) THEN
1054 nps=igdtmplo(8)
1055 rlat1=float(igdtmplo(10))*1.e-6
1056 rlon1=float(igdtmplo(11))*1.e-6
1057 orient=float(igdtmplo(14))*1.e-6
1058 xmesh=float(igdtmplo(15))*1.e-3
1059 iproj=mod(igdtmplo(17)/128,2)
1060 ip=(nps+1)/2
1061 h=(-1.)**iproj
1062 slat=float(abs(igdtmplo(13)))*1.e-6
1063 CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
1064 de=(1.+sin(slat/dpr))*rerth
1065 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
1066 xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
1067 yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
1068 IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
1069 IF(iproj.EQ.0) THEN
1070 CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
1071 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1072 slat,xmesh,orient,ui,vi,.true.,uo,vo,uo2,vo2, &
1073 .false.,dumm,dumm,dumm,dumm, &
1074 .false.,dumm,dumm,dumm,dumm)
1075 ELSE
1076 CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
1077 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1078 slat,xmesh,orient,ui,vi,.true.,uo2,vo2,uo,vo, &
1079 .false.,dumm,dumm,dumm,dumm, &
1080 .false.,dumm,dumm,dumm,dumm)
1081 ENDIF
1082 ispec=1
1083 ENDIF
1084 ! SPECIAL CASE OF MERCATOR GRID
1085 ELSEIF(igdtnumo.EQ.10) THEN
1086 ni=igdtmplo(8)
1087 nj=igdtmplo(9)
1088 rlat1=float(igdtmplo(10))*1.0e-6
1089 rlon1=float(igdtmplo(11))*1.0e-6
1090 rlon2=float(igdtmplo(15))*1.0e-6
1091 rlati=float(igdtmplo(13))*1.0e-6
1092 iscano=mod(igdtmplo(16)/128,2)
1093 jscano=mod(igdtmplo(16)/64,2)
1094 nscano=mod(igdtmplo(16)/32,2)
1095 dy=float(igdtmplo(19))*1.0e-3
1096 hi=(-1.)**iscano
1097 hj=(-1.)**(1-jscano)
1098 CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
1099 dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
1100 dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
1101 IF(nscano.EQ.0) THEN
1102 CALL sptrunmv(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
1103 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1104 rlat1,rlon1,dlato,dlono,ui,vi, &
1105 .true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
1106 ispec=1
1107 ENDIF
1108 ENDIF
1109 ! GENERAL SLOW CASE
1110 IF(ispec.EQ.0) THEN
1111 CALL sptrungv(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
1112 iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon, &
1113 ui,vi,.true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
1114 DO k=1,km
1115 ibo(k)=0
1116 DO n=1,no
1117 lo(n,k)=.true.
1118 urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
1119 vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
1120 uo(n,k)=urot
1121 vo(n,k)=vrot
1122 ENDDO
1123 ENDDO
1124 ENDIF
1125 ELSE
1126 DO k=1,km
1127 ibo(k)=1
1128 DO n=1,no
1129 lo(n,k)=.false.
1130 uo(n,k)=0.
1131 vo(n,k)=0.
1132 ENDDO
1133 ENDDO
1134 ENDIF
1135 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1136 END SUBROUTINE polatev4_grib2
1137
1222 SUBROUTINE polatev4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,UI,VI, &
1223 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
1224 INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
1225 INTEGER, INTENT(IN ) :: KM, MI, MO
1226 INTEGER, INTENT( OUT) :: IRET, IBO(KM)
1227 INTEGER, INTENT(IN) :: KGDSI(200),KGDSO(200)
1228 !
1229 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
1230 !
1231 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
1232 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
1233 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
1234 REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
1235 !
1236 REAL, PARAMETER :: FILL=-9999.
1237 REAL, PARAMETER :: RERTH=6.3712e6
1238 REAL, PARAMETER :: PI=3.14159265358979
1239 REAL, PARAMETER :: DPR=180./pi
1240 !
1241 INTEGER :: IDRTO, IROMB, ISKIPI, ISPEC
1242 INTEGER :: IDRTI, IMAXI, JMAXI, IM, JM
1243 INTEGER :: IPRIME, IG, IMO, JMO, IGO, JGO
1244 INTEGER :: ISCAN, JSCAN, NSCAN
1245 INTEGER :: ISCANO, JSCANO, NSCANO
1246 INTEGER :: IP, IPROJ, JSKIPI, JG
1247 INTEGER :: K, MAXWV, N, NI, NJ, NO, NPS
1248 !
1249 REAL :: DLAT, DLON, DLATO, DLONO, DE, DR, DY
1250 REAL :: H, HI, HJ, DUMM(1)
1251 REAL :: ORIENT
1252 REAL :: RLAT1, RLON1, RLAT2, RLON2, RLATI
1253 REAL :: UROT, VROT, UO2(MO,KM),VO2(MO,KM)
1254 REAL :: XMESH, XP, YP, XPTS(MO),YPTS(MO)
1255
1256 type(grib1_descriptor) :: desc_in, desc_out
1257 class(ip_grid), allocatable :: grid_in, grid_out
1258
1259 desc_in = init_descriptor(kgdsi)
1260 desc_out = init_descriptor(kgdso)
1261
1262 call init_grid(grid_in, desc_in)
1263 call init_grid(grid_out, desc_out)
1264 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1265 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
1266 iret=0
1267 IF(kgdso(1).GE.0) THEN
1268 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
1269 IF(no.EQ.0) iret=3
1270 ENDIF
1271 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1272 ! AFFIRM APPROPRIATE INPUT GRID
1273 ! LAT/LON OR GAUSSIAN
1274 ! NO BITMAPS
1275 ! FULL ZONAL COVERAGE
1276 ! FULL MERIDIONAL COVERAGE
1277 idrti=kgdsi(1)
1278 im=kgdsi(2)
1279 jm=kgdsi(3)
1280 rlon1=kgdsi(5)*1.e-3
1281 rlon2=kgdsi(8)*1.e-3
1282 iscan=mod(kgdsi(11)/128,2)
1283 jscan=mod(kgdsi(11)/64,2)
1284 nscan=mod(kgdsi(11)/32,2)
1285 IF(idrti.NE.0.AND.idrti.NE.4) iret=41
1286 DO k=1,km
1287 IF(ibi(k).NE.0) iret=41
1288 ENDDO
1289 IF(iret.EQ.0) THEN
1290 IF(iscan.EQ.0) THEN
1291 dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
1292 ELSE
1293 dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
1294 ENDIF
1295 ig=nint(360/abs(dlon))
1296 iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
1297 imaxi=ig
1298 jmaxi=jm
1299 IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
1300 ENDIF
1301 IF(iret.EQ.0.AND.idrti.EQ.0) THEN
1302 rlat1=kgdsi(4)*1.e-3
1303 rlat2=kgdsi(7)*1.e-3
1304 dlat=(rlat2-rlat1)/(jm-1)
1305 jg=nint(180/abs(dlat))
1306 IF(jm.EQ.jg) idrti=256
1307 IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
1308 ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
1309 jg=kgdsi(10)*2
1310 IF(jm.NE.jg) iret=41
1311 ENDIF
1312 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1313 ! SET PARAMETERS
1314 IF(iret.EQ.0) THEN
1315 iromb=ipopt(1)
1316 maxwv=ipopt(2)
1317 IF(maxwv.EQ.-1) THEN
1318 IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
1319 IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
1320 IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
1321 IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
1322 IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
1323 IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
1324 ENDIF
1325 IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
1326 ENDIF
1327 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1328 ! INTERPOLATE
1329 IF(iret.EQ.0) THEN
1330 IF(nscan.EQ.0) THEN
1331 iskipi=1
1332 jskipi=im
1333 ELSE
1334 iskipi=jm
1335 jskipi=1
1336 ENDIF
1337 IF(iscan.EQ.1) iskipi=-iskipi
1338 IF(jscan.EQ.0) jskipi=-jskipi
1339 ispec=0
1340 ! SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
1341 IF((kgdso(1).EQ.0.OR.kgdso(1).EQ.4).AND. &
1342 mod(kgdso(2),2).EQ.0.AND.kgdso(5).EQ.0.AND. &
1343 kgdso(11).EQ.0) THEN
1344 idrto=kgdso(1)
1345 imo=kgdso(2)
1346 jmo=kgdso(3)
1347 rlon2=kgdso(8)*1.e-3
1348 dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
1349 igo=nint(360/abs(dlono))
1350 IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
1351 rlat1=kgdso(4)*1.e-3
1352 rlat2=kgdso(7)*1.e-3
1353 dlat=(rlat2-rlat1)/(jmo-1)
1354 jgo=nint(180/abs(dlat))
1355 IF(jmo.EQ.jgo) idrto=256
1356 IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
1357 ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
1358 jgo=kgdso(10)*2
1359 IF(jmo.EQ.jgo) ispec=1
1360 ENDIF
1361 IF(ispec.EQ.1) THEN
1362 CALL sptrunv(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
1363 km,iprime,iskipi,jskipi,mi,0,0,mo,0,ui,vi, &
1364 .true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
1365 ENDIF
1366 ! SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
1367 ELSEIF(kgdso(1).EQ.5.AND. &
1368 kgdso(2).EQ.kgdso(3).AND.mod(kgdso(2),2).EQ.1.AND. &
1369 kgdso(8).EQ.kgdso(9).AND.kgdso(11).EQ.64.AND. &
1370 mod(kgdso(6)/8,2).EQ.1) THEN
1371 nps=kgdso(2)
1372 rlat1=kgdso(4)*1.e-3
1373 rlon1=kgdso(5)*1.e-3
1374 orient=kgdso(7)*1.e-3
1375 xmesh=kgdso(8)
1376 iproj=mod(kgdso(10)/128,2)
1377 ip=(nps+1)/2
1378 h=(-1.)**iproj
1379 de=(1.+sin(60./dpr))*rerth
1380 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
1381 xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
1382 yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
1383 IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
1384 IF(iproj.EQ.0) THEN
1385 CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
1386 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1387 60.,xmesh,orient,ui,vi,.true.,uo,vo,uo2,vo2, &
1388 .false.,dumm,dumm,dumm,dumm, &
1389 .false.,dumm,dumm,dumm,dumm)
1390 ELSE
1391 CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
1392 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1393 60.,xmesh,orient,ui,vi,.true.,uo2,vo2,uo,vo, &
1394 .false.,dumm,dumm,dumm,dumm, &
1395 .false.,dumm,dumm,dumm,dumm)
1396 ENDIF
1397 ispec=1
1398 ENDIF
1399 ! SPECIAL CASE OF MERCATOR GRID
1400 ELSEIF(kgdso(1).EQ.1) THEN
1401 ni=kgdso(2)
1402 nj=kgdso(3)
1403 rlat1=kgdso(4)*1.e-3
1404 rlon1=kgdso(5)*1.e-3
1405 rlon2=kgdso(8)*1.e-3
1406 rlati=kgdso(9)*1.e-3
1407 iscano=mod(kgdso(11)/128,2)
1408 jscano=mod(kgdso(11)/64,2)
1409 nscano=mod(kgdso(11)/32,2)
1410 dy=kgdso(13)
1411 hi=(-1.)**iscano
1412 hj=(-1.)**(1-jscano)
1413 dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
1414 dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
1415 IF(nscano.EQ.0) THEN
1416 CALL sptrunmv(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
1417 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1418 rlat1,rlon1,dlato,dlono,ui,vi, &
1419 .true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
1420 ispec=1
1421 ENDIF
1422 ENDIF
1423 ! GENERAL SLOW CASE
1424 IF(ispec.EQ.0) THEN
1425 CALL sptrungv(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
1426 iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon, &
1427 ui,vi,.true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
1428 DO k=1,km
1429 ibo(k)=0
1430 DO n=1,no
1431 lo(n,k)=.true.
1432 urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
1433 vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
1434 uo(n,k)=urot
1435 vo(n,k)=vrot
1436 ENDDO
1437 ENDDO
1438 ENDIF
1439 ELSE
1440 DO k=1,km
1441 ibo(k)=1
1442 DO n=1,no
1443 lo(n,k)=.false.
1444 uo(n,k)=0.
1445 vo(n,k)=0.
1446 ENDDO
1447 ENDDO
1448 ENDIF
1449 END SUBROUTINE polatev4_grib1
1450end module spectral_interp_mod
gdswzd_mod::gdswzd
Definition gdswzd_mod.F90:37
ip_grid_descriptor_mod::init_descriptor
Definition ip_grid_descriptor_mod.F90:63
ip_grid_factory_mod::init_grid
Definition ip_grid_factory_mod.F90:20
spectral_interp_mod::interpolate_spectral
Definition spectral_interp_mod.F90:20
spectral_interp_mod::polates4
Definition spectral_interp_mod.F90:25
spectral_interp_mod::polatev4
Definition spectral_interp_mod.F90:30
earth_radius_mod
Determine earth radius and shape.
Definition earth_radius_mod.F90:7
earth_radius_mod::earth_radius
subroutine, public earth_radius(igdtmpl, igdtlen, radius, eccen_squared)
Determine earth radius and shape.
Definition earth_radius_mod.F90:41
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_grid_mod
Abstract ip_grid type.
Definition ip_grid_mod.F90:10
spectral_interp_mod
Interpolate spectral.
Definition spectral_interp_mod.F90:8
spectral_interp_mod::polates4_grib1
subroutine polates4_grib1(ipopt, kgdsi, kgdso, mi, mo, km, ibi, gi, no, rlat, rlon, ibo, lo, go, iret)
Interpolate scalar fields (spectral).
Definition spectral_interp_mod.F90:562
spectral_interp_mod::polatev4_grib1
subroutine polatev4_grib1(ipopt, kgdsi, kgdso, mi, mo, km, ibi, ui, vi, no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret)
Interpolate vector fields (spectral).
Definition spectral_interp_mod.F90:1224
spectral_interp_mod::interpolate_spectral_scalar
subroutine interpolate_spectral_scalar(ipopt, grid_in, grid_out, mi, mo, km, ibi, gi, no, rlat, rlon, ibo, lo, go, iret)
Interpolate spectral scalar.
Definition spectral_interp_mod.F90:64
spectral_interp_mod::polates4_grib2
subroutine polates4_grib2(ipopt, igdtnumi, igdtmpli, igdtleni, igdtnumo, igdtmplo, igdtleno, mi, mo, km, ibi, gi, no, rlat, rlon, ibo, lo, go, iret)
Interpolate scalar fields (spectral).
Definition spectral_interp_mod.F90:259
spectral_interp_mod::interpolate_spectral_vector
subroutine interpolate_spectral_vector(ipopt, grid_in, grid_out, mi, mo, km, ibi, ui, vi, no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret)
Interpolate spectral vector.
Definition spectral_interp_mod.F90:125
spectral_interp_mod::polatev4_grib2
subroutine polatev4_grib2(ipopt, igdtnumi, igdtmpli, igdtleni, igdtnumo, igdtmplo, igdtleno, mi, mo, km, ibi, ui, vi, no, rlat, rlon, crot, srot, ibo, lo, uo, vo, iret)
Interpolate vector fields (spectral).
Definition spectral_interp_mod.F90:892
sptrun
subroutine sptrun(iromb, maxwv, idrti, imaxi, jmaxi, idrto, imaxo, jmaxo, kmax, iprime, iskipi, jskipi, kskipi, iskipo, jskipo, kskipo, jcpu, gridi, grido)
This subprogram spectrally truncates scalar fields on a global cylindrical grid, returning the fields...
Definition sptrun.f:58
sptrung
subroutine sptrung(iromb, maxwv, idrti, imaxi, jmaxi, kmax, nmax, iprime, iskipi, jskipi, kskipi, kgskip, nrskip, ngskip, jcpu, rlat, rlon, gridi, gp)
This subprogram spectrally truncates scalar fields on a global cylindrical grid, returning the fields...
Definition sptrung.f:68
sptrungv
subroutine sptrungv(iromb, maxwv, idrti, imaxi, jmaxi, kmax, nmax, iprime, iskipi, jskipi, kskipi, kgskip, nrskip, ngskip, jcpu, rlat, rlon, gridui, gridvi, luv, up, vp, ldz, dp, zp, lps, pp, sp)
THIS SUBPROGRAM SPECTRALLY TRUNCATES VECTORS FIELDS ON A GLOBAL CYLINDRICAL GRID, RETURNING THE FIELD...
Definition sptrungv.f:85
sptrunm
subroutine sptrunm(iromb, maxwv, idrti, imaxi, jmaxi, kmax, mi, mj, iprime, iskipi, jskipi, kskipi, kgskip, niskip, njskip, jcpu, rlat1, rlon1, dlat, dlon, gridi, gm)
This subprogram spectrally truncates scalar fields on a global cylindrical grid, returning the fields...
Definition sptrunm.f:80
sptrunmv
subroutine sptrunmv(iromb, maxwv, idrti, imaxi, jmaxi, kmax, mi, mj, iprime, iskipi, jskipi, kskipi, kgskip, niskip, njskip, jcpu, rlat1, rlon1, dlat, dlon, gridui, gridvi, luv, um, vm, ldz, dm, zm, lps, pm, sm)
THIS SUBPROGRAM SPECTRALLY TRUNCATES VECTOR FIELDS ON A GLOBAL CYLINDRICAL GRID, RETURNING THE FIELDS...
Definition sptrunmv.f:96
sptruns
subroutine sptruns(iromb, maxwv, idrti, imaxi, jmaxi, kmax, nps, iprime, iskipi, jskipi, kskipi, kgskip, niskip, njskip, jcpu, true, xmesh, orient, gridi, gn, gs)
This subprogram spectrally truncates scalar fields on a global cylindrical grid, returning the fields...
Definition sptruns.f:75
sptrunsv
subroutine sptrunsv(iromb, maxwv, idrti, imaxi, jmaxi, kmax, nps, iprime, iskipi, jskipi, kskipi, kgskip, niskip, njskip, jcpu, true, xmesh, orient, gridui, gridvi, luv, un, vn, us, vs, ldz, dn, zn, ds, zs, lps, pn, sn, ps, ss)
This subprogram spectrally truncates vector fields on a global cylindrical grid, returning the fields...
Definition sptrunsv.f:94
sptrunv
subroutine sptrunv(iromb, maxwv, idrti, imaxi, jmaxi, idrto, imaxo, jmaxo, kmax, iprime, iskipi, jskipi, kskipi, iskipo, jskipo, kskipo, jcpu, gridui, gridvi, luv, griduo, gridvo, ldz, griddo, gridzo, lps, gridpo, gridso)
This subprogram spectrally truncates vector fields on a global cylindrical grid, returning the fields...
Definition sptrunv.f:96
ip_grid_descriptor_mod::grib1_descriptor
Descriptor representing a grib1 grib descriptor section (GDS) with an integer array.
Definition ip_grid_descriptor_mod.F90:39
ip_grid_descriptor_mod::grib2_descriptor
Grib-2 descriptor containing a grib2 GDT represented by an integer array.
Definition ip_grid_descriptor_mod.F90:49
ip_grid_mod::ip_grid
Abstract grid that holds fields and methods common to all grids.
Definition ip_grid_mod.F90:58