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