NCEPLIBS-ip  5.1.0
spectral_interp_mod.F90
Go to the documentation of this file.
1 
4 
8 module 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 
35 contains
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
1450 end module spectral_interp_mod
gdswzd_mod::gdswzd
Definition: gdswzd_mod.F90:37
ip_grid_descriptor_mod::init_descriptor
Definition: ip_grid_descriptor_mod.F90:60
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::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
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::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::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_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
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:38
ip_grid_descriptor_mod::grib2_descriptor
Grib-2 descriptor containing a grib2 GDT represented by an integer array.
Definition: ip_grid_descriptor_mod.F90:47
ip_grid_mod::ip_grid
Abstract grid that holds fields and methods common to all grids.
Definition: ip_grid_mod.F90:58