NCEPLIBS-ip  5.0.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
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:52