spectral_interp_mod.F90

Go to the documentation of this file.

!> @file
!! @brief Interpolate spectral.
!! @author Mark Iredell @date 96-04-10
 
!> @brief Interpolate spectral.
!!
!! @author Mark Iredell @date 96-04-10
module spectral_interp_mod
  use gdswzd_mod
  use ip_grid_mod
  use ip_grid_descriptor_mod
  use ip_grid_factory_mod
  use earth_radius_mod
  use sp_mod
  implicit none
 
  private
  public :: interpolate_spectral
 
  interface interpolate_spectral
     module procedure interpolate_spectral_scalar
     module procedure interpolate_spectral_vector
  interface interpolate_spectral …
  end interface interpolate_spectral
 
  interface polates4
     module procedure polates4_grib1
     module procedure polates4_grib2
  interface polates4 …
  end interface polates4
 
  interface polatev4
     module procedure polatev4_grib1
     module procedure polatev4_grib2
  interface polatev4 …
  end interface polatev4
 
contains
 
  !> Interpolate spectral scalar.
  !>
  !> @param[in] ipopt interpolation options; ipopt(1)=0 for triangular;
  !> ipopt(1)=1 for rhomboidal; ipopt(2) is truncation number
  !> (defaults to a sensible truncation if ipopt(2)=-1).
  !> @param[in] grid_in input grid descriptor.
  !> @param[in] grid_out output grid descriptor.
  !> @param[in] MI skip number between input grid fields if km>1 or
  !> dimension of input grid fields if km=1.
  !> @param[in] MO skip number between output grid fields if km>1 or
  !> dimension of output grid fields if km=1.
  !> @param[in] KM number of fields to interpolate.
  !> @param[in] IBI input bitmap flags (Must be all 0. Routine does
  !> not do bitmapped interpolation.)
  !> @param[in] GI input fields to interpolate.
  !> @param[out] NO number of output points.
  !> @param[inout] RLAT output latitudes in degrees.
  !> @param[inout] RLON output longitudes in degrees.
  !> @param[out] IBO output bitmap flags.
  !> @param[out] LO output bitmaps.
  !> @param[out] GO output fields interpolated.
  !> @param[out] IRET return code. 0/non-0 - successful/not successful.
  !>
  !! @author Mark Iredell @date 96-04-10
  subroutine interpolate_spectral_scalar(IPOPT,grid_in,grid_out, &
       MI,MO,KM,IBI,GI, &
       NO,RLAT,RLON,IBO,LO,GO,IRET)
    INTEGER,          INTENT(IN   ) :: IPOPT(20)
    class(ip_grid), intent(in) :: grid_in, grid_out
    INTEGER,          INTENT(IN   ) :: MI, MO
    INTEGER,          INTENT(IN   ) :: IBI(KM), KM
    INTEGER,          INTENT(  OUT) :: IBO(KM), IRET, NO
    !
    LOGICAL*1,        INTENT(  OUT) :: LO(MO,KM)
    !
    REAL,             INTENT(IN   ) :: GI(MI,KM)
    REAL,             INTENT(INOUT) :: RLAT(MO),RLON(MO)
    REAL,             INTENT(  OUT) :: GO(MO,KM)
 
 
    select type(desc_in => grid_in%descriptor)
    type is(grib1_descriptor)
       select type(desc_out => grid_out%descriptor)
       type is(grib1_descriptor)
          CALL polates4(ipopt,desc_in%gds,desc_out%gds,mi,mo,km,ibi,gi,no,rlat,rlon,ibo,lo,go,iret)
       end select
 
    type is(grib2_descriptor)
       select type(desc_out => grid_out%descriptor)
       type is(grib2_descriptor)
          CALL polates4(ipopt,desc_in%gdt_num,desc_in%gdt_tmpl,desc_in%gdt_len, &
               desc_out%gdt_num,desc_out%gdt_tmpl,desc_out%gdt_len, &
               mi,mo,km,ibi,gi,no,rlat,rlon,ibo,lo,go,iret)
       end select
    end select
  subroutine interpolate_spectral_scalar(IPOPT,grid_in,grid_out, & …
  end subroutine interpolate_spectral_scalar
 
  !> Interpolate spectral vector.
  !>
  !> @param ipopt interpolation options; ipopt(1)=0 for triangular;
  !> ipopt(1)=1 for rhomboidal; ipopt(2) is truncation number
  !> (defaults to a sensible truncation if ipopt(2)=-1).
  !> @param grid_in input grid descriptor.
  !> @param grid_out output grid descriptor.
  !> @param MI skip number between input grid fields if km>1 or
  !> dimension of input grid fields if km=1.
  !> @param MO skip number between output grid fields if km>1 or
  !> dimension of output grid fields if km=1.
  !> @param KM number of fields to interpolate.
  !> @param IBI input bitmap flags (Must be all 0. Routine does not do
  !> bitmapped interpolation.)
  !> @param UI input u-component fields to interpolate.
  !> @param VI input v-component fields to interpolate.
  !> @param NO number of output points.
  !> @param RLAT output latitudes in degrees.
  !> @param RLON output longitudes in degrees.
  !> @param CROT vector rotation cosines.
  !> @param SROT vector rotation sines.
  !> @param IBO output bitmap flags.
  !> @param LO output bitmaps.
  !> @param UO output u-component fields interpolated.
  !> @param VO output v-component fields interpolated.
  !> @param IRET return code. 0/non-0 - successful/not successful.
  !>
  !! @author Mark Iredell @date 96-04-10
  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)
    class(ip_grid), intent(in) :: grid_in, grid_out
    INTEGER,          INTENT(IN   ) :: IPOPT(20), IBI(KM)
    INTEGER,          INTENT(IN   ) :: KM, MI, MO
    INTEGER,          INTENT(  OUT) :: IRET, IBO(KM), NO
    !
    LOGICAL*1,        INTENT(  OUT) :: LO(MO,KM)
    !
    REAL,             INTENT(IN   ) :: UI(MI,KM),VI(MI,KM)
    REAL,             INTENT(  OUT) :: UO(MO,KM),VO(MO,KM)
    REAL,             INTENT(INOUT) :: RLAT(MO),RLON(MO)
    REAL,             INTENT(  OUT) :: CROT(MO),SROT(MO)
 
 
    select type(desc_in => grid_in%descriptor)
    type is(grib1_descriptor)
       select type(desc_out => grid_out%descriptor)
       type is(grib1_descriptor)
         CALL polatev4_grib1(ipopt,desc_in%gds,desc_out%gds,mi,mo,km,ibi,ui,vi,&
            no,rlat,rlon,crot,srot,ibo,lo,uo,vo,iret)
       end select
 
    type is(grib2_descriptor)
       select type(desc_out => grid_out%descriptor)
       type is(grib2_descriptor)
           CALL polatev4(ipopt,desc_in%gdt_num,desc_in%gdt_tmpl,desc_in%gdt_len, &
            desc_out%gdt_num,desc_out%gdt_tmpl,desc_out%gdt_len, &
            mi,mo,km,ibi,ui,vi,&
            no,rlat,rlon,crot,srot,ibo,lo,uo,vo,iret)
       end select
    end select
 
 
  subroutine interpolate_spectral_vector(IPOPT,grid_in,grid_out, & …
  end subroutine interpolate_spectral_vector
! @author Mark Iredell @date 96-04-10
  !> Interpolate scalar fields (spectral).
  !>
  !> This subprogram performs spectral interpolation from any grid to
  !> any grid for scalar fields. It requires that the input fields be
  !> uniformly global. Options allow choices between triangular shape
  !> (ipopt(1)=0) and rhomboidal shape (ipopt(1)=1) which has no
  !> default; a second option is the truncation (ipopt(2)) which
  !> defaults to a sensible truncation for the input grid (if
  !> opt(2)=-1).
  !>
  !> @note If the output grid is not found in a special list, then the
  !> transform back to grid is not very fast. This special list
  !> contains global cylindrical grids, polar stereographic grids
  !> centered at the pole and mercator grids.
  !>
  !> Only horizontal interpolation is performed.
  !>
  !> The code recognizes the following projections, where "igdtnumi/o"
  !> is the GRIB2 grid defintion template number for the input and
  !> onutput grids, respectively:
  !> - igdtnumi/o = 00 equidistant cylindrical
  !> - igdtnumo = 01 rotated equidistant cylindrical. "e" and non-"e" staggered
  !> - igdtnumo = 10 mercator cylindrical
  !> - igdtnumo = 20 polar stereographic azimuthal
  !> - igdtnumo = 30 lambert conformal conical
  !> - igdtnumi/o = 40 gaussian cylindrical
  !>
  !> As an added bonus the number of output grid points and their
  !> latitudes and longitudes are also returned. On the other hand,
  !> the output can be a set of station points if igdtnumo < 0, in which
  !> case the number of points and their latitudes and longitudes must
  !> be input. Output bitmaps will not be created.
  !>        
  !> ### Program History Log
  !> Date | Programmer | Comments
  !> -----|------------|---------
  !>   96-04-10 | Iredell | initial
  !> 2001-06-18 | Iredell | improve detection of special fast transform
  !> 2015-01-27 | Gayno | replace calls to gdswiz with new merged version of gdswzd.
  !> 2015-07-13 | Gayno | convert to grib 2. replace grib 1 kgds arrays with grib 2 grid definition template arrays.
  !>
  !> @param[in] ipopt (20) interpolation options; ipopt(1)=0 for
  !> triangular, ipopt(1)=1 for rhomboidal; ipopt(2) is truncation
  !> number (defaults to sensible if ipopt(2)=-1).
  !> @param[in] igdtnumi grid definition template number - input
  !> grid. Corresponds to the gfld%igdtnum component of the
  !> [NCEPLIBS-g2](https://github.com/NOAA-EMC/NCEPLIBS-g2) library
  !> gridmod data structure.
  !> - 00 - equidistant cylindrical
  !> - 01 - rotated equidistant cylindrical.  "e" and non-"e" staggered
  !> - 10 - mercator cyclindrical
  !> - 20 - polar stereographic azimuthal
  !> - 30 - lambert conformal conical
  !> - 40 - gaussian equidistant cyclindrical
  !> @param[in] igdtmpli (igdtleni) grid definition template array -
  !> input grid. corresponds to the gfld%igdtmpl component of the ncep
  !> g2 library gridmod data structure: (section 3 info). See comments
  !> in routine ipolates() for complete definition.
  !> @param[in] igdtleni number of elements of the grid definition
  !> template array - input grid.  corresponds to the gfld%igdtlen
  !> component of the ncep g2 library gridmod data structure.
  !> @param[in] igdtnumo grid definition template number - output
  !> grid. Corresponds to the gfld%igdtnum component of the
  !> [NCEPLIBS-g2](https://github.com/NOAA-EMC/NCEPLIBS-g2) library
  !> gridmod data structure. igdtnumo<0 means interpolate to random
  !> station points. Otherwise, same definition as igdtnumi.
  !> @param[in] igdtmplo (igdtleno) grid definition template array -
  !> output grid. Corresponds to the gfld%igdtmpl component of the
  !> ncep g2 library gridmod data structure (section 3 info). See
  !> comments in routine ipolates() for complete definition.
  !> @param[in] igdtleno number of elements of the grid definition
  !> template array - output grid. Corresponds to the gfld%igdtlen
  !> component of the
  !> [NCEPLIBS-g2](https://github.com/NOAA-EMC/NCEPLIBS-g2) library
  !> gridmod data structure.
  !> @param[in] mi skip number between input grid fields if km>1 or
  !> dimension of input grid fields if km=1
  !> @param[in] mo skip number between output grid fields if km>1 or
  !> dimension of output grid fields if km=1
  !> @param[out] km number of fields to interpolate
  !> @param[out] ibi (km) input bitmap flags (must be all 0)
  !> @param[out] gi (mi,km) input fields to interpolate
  !> @param[out] no number of output points (only if igdtnumo>=0)
  !> @param[out] rlat (mo) output latitudes in degrees (if igdtnumo<0)
  !> @param[out] rlon (mo) output longitudes in degrees (if igdtnumo<0)
  !> @param[out] ibo (km) output bitmap flags
  !> @param[out] lo (mo,km) output bitmaps (always output)
  !> @param[out] go (mo,km) output fields interpolated
  !> @param[out] iret return code
  !> - 0 successful interpolation
  !> - 2 unrecognized input grid or no grid overlap
  !> - 3 unrecognized output grid
  !> - 41 invalid nonglobal input grid
  !> - 42 invalid spectral method parameters
  !>
  !>
  SUBROUTINE polates4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
       IGDTNUMO,IGDTMPLO,IGDTLENO, &
       MI,MO,KM,IBI,GI, &
       NO,RLAT,RLON,IBO,LO,GO,IRET)
    INTEGER,          INTENT(IN   ) :: IGDTNUMI, IGDTLENI
    INTEGER,          INTENT(IN   ) :: IGDTMPLI(IGDTLENI)
    INTEGER,          INTENT(IN   ) :: IGDTNUMO, IGDTLENO
    INTEGER,          INTENT(IN   ) :: IGDTMPLO(IGDTLENO)
    INTEGER,          INTENT(IN   ) :: IPOPT(20)
    INTEGER,          INTENT(IN   ) :: MI, MO
    INTEGER,          INTENT(IN   ) :: IBI(KM), KM
    INTEGER,          INTENT(  OUT) :: IBO(KM), IRET, NO
    !
    LOGICAL*1,        INTENT(  OUT) :: LO(MO,KM)
    !
    REAL,             INTENT(IN   ) :: GI(MI,KM)
    REAL,             INTENT(INOUT) :: RLAT(MO),RLON(MO)
    REAL,             INTENT(  OUT) :: GO(MO,KM)
    !
    REAL,             PARAMETER     :: FILL=-9999.
    REAL,             PARAMETER     :: PI=3.14159265358979
    REAL,             PARAMETER     :: DPR=180./pi
    !
    INTEGER                         :: IDRTI, IDRTO, IG, JG, IM, JM
    INTEGER                         :: IGO, JGO, IMO, JMO
    INTEGER                         :: ISCAN, JSCAN, NSCAN
    INTEGER                         :: ISCANO, JSCANO, NSCANO
    INTEGER                         :: ISKIPI, JSKIPI, ISCALE
    INTEGER                         :: IMAXI, JMAXI, ISPEC
    INTEGER                         :: IP, IPRIME, IPROJ, IROMB, K
    INTEGER                         :: MAXWV, N, NI, NJ, NPS
    !
    REAL                            :: DE, DR, DY
    REAL                            :: DLAT, DLON, DLATO, DLONO
    REAL                            :: GO2(MO,KM), H, HI, HJ
    REAL                            :: ORIENT, SLAT, RERTH, E2
    REAL                            :: RLAT1, RLON1, RLAT2, RLON2, RLATI
    REAL                            :: XMESH, XP, YP
    REAL                            :: XPTS(MO), YPTS(MO)
 
    type(grib2_descriptor) :: desc_in, desc_out
    class(ip_grid), allocatable :: grid_in, grid_out
 
    desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
    desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
 
    call init_grid(grid_in, desc_in)
    call init_grid(grid_out, desc_out)
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
    iret=0
    IF(igdtnumo.GE.0) THEN
       !CALL GDSWZD(IGDTNUMO,IGDTMPLO,IGDTLENO, 0,MO,FILL,XPTS,YPTS,RLON,RLAT,NO)
       CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
       IF(no.EQ.0) iret=3
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  AFFIRM APPROPRIATE INPUT GRID
    !    LAT/LON OR GAUSSIAN
    !    NO BITMAPS
    !    FULL ZONAL COVERAGE
    !    FULL MERIDIONAL COVERAGE
    idrti=igdtnumi
    IF(idrti==40) idrti=4
    IF(idrti==0.OR.idrti==4)THEN
       im=igdtmpli(8)
       jm=igdtmpli(9)
       iscale=igdtmpli(10)*igdtmpli(11)
       IF(iscale==0) iscale=10**6
       rlon1=float(igdtmpli(13))/float(iscale)
       rlon2=float(igdtmpli(16))/float(iscale)
       iscan=mod(igdtmpli(19)/128,2)
       jscan=mod(igdtmpli(19)/64,2)
       nscan=mod(igdtmpli(19)/32,2)
    ELSE
       iret=41
    ENDIF
    DO k=1,km
       IF(ibi(k).NE.0) iret=41
    ENDDO
    IF(iret.EQ.0) THEN
       IF(iscan.EQ.0) THEN
          dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
       ELSE
          dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
       ENDIF
       ig=nint(360/abs(dlon))
       iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
       imaxi=ig
       jmaxi=jm
       IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
    ENDIF
    IF(iret.EQ.0.AND.idrti.EQ.0) THEN
       iscale=igdtmpli(10)*igdtmpli(11)
       IF(iscale==0) iscale=10**6
       rlat1=float(igdtmpli(12))/float(iscale)
       rlat2=float(igdtmpli(15))/float(iscale)
       dlat=(rlat2-rlat1)/(jm-1)
       jg=nint(180/abs(dlat))
       IF(jm.EQ.jg) idrti=256
       IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
    ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
       jg=igdtmpli(18)*2
       IF(jm.NE.jg) iret=41
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  SET PARAMETERS
    IF(iret.EQ.0) THEN
       iromb=ipopt(1)
       maxwv=ipopt(2)
       IF(maxwv.EQ.-1) THEN
          IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
          IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
          IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
          IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
          IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
          IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
       ENDIF
       IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  INTERPOLATE
    IF(iret.EQ.0) THEN
       IF(nscan.EQ.0) THEN
          iskipi=1
          jskipi=im
       ELSE
          iskipi=jm
          jskipi=1
       ENDIF
       IF(iscan.EQ.1) iskipi=-iskipi
       IF(jscan.EQ.0) jskipi=-jskipi
       ispec=0
       !  SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
       IF((igdtnumo.EQ.0.OR.igdtnumo.EQ.40).AND. &
            mod(igdtmplo(8),2).EQ.0.AND.igdtmplo(13).EQ.0.AND.igdtmplo(19).EQ.0) THEN
          idrto=igdtnumo
          IF(idrto==40)idrto=4
          imo=igdtmplo(8)
          jmo=igdtmplo(9)
          iscale=igdtmplo(10)*igdtmplo(11)
          IF(iscale==0) iscale=10**6
          rlon2=float(igdtmplo(16))/float(iscale)
          dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
          igo=nint(360/abs(dlono))
          IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
             rlat1=float(igdtmplo(12))/float(iscale)
             rlat2=float(igdtmplo(15))/float(iscale)
             dlat=(rlat2-rlat1)/(jmo-1)
             jgo=nint(180/abs(dlat))
             IF(jmo.EQ.jgo) idrto=256
             IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
          ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
             jgo=igdtmplo(18)*2
             IF(jmo.EQ.jgo) ispec=1
          ENDIF
          IF(ispec.EQ.1) THEN
             CALL sptrun(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
                  km,iprime,iskipi,jskipi,mi,0,0,mo,0,gi,go)
          ENDIF
          !  SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
       ELSEIF(igdtnumo.EQ.20.AND. &
            igdtmplo(8).EQ.igdtmplo(9).AND.mod(igdtmplo(8),2).EQ.1.AND. &
            igdtmplo(15).EQ.igdtmplo(16).AND.igdtmplo(18).EQ.64) THEN
          nps=igdtmplo(8)
          rlat1=float(igdtmplo(10))*1.e-6
          rlon1=float(igdtmplo(11))*1.e-6
          orient=float(igdtmplo(14))*1.e-6
          xmesh=float(igdtmplo(15))*1.e-3
          iproj=mod(igdtmplo(17)/128,2)
          ip=(nps+1)/2
          h=(-1.)**iproj
          slat=float(abs(igdtmplo(13)))*1.e-6
          CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
          de=(1.+sin(slat/dpr))*rerth
          dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
          xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
          yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
          IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
             IF(iproj.EQ.0) THEN
                CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
                     iprime,iskipi,jskipi,mi,mo,0,0,0, &
                     slat,xmesh,orient,gi,go,go2)
             ELSE
                CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
                     iprime,iskipi,jskipi,mi,mo,0,0,0, &
                     slat,xmesh,orient,gi,go2,go)
             ENDIF
             ispec=1
          ENDIF
          !  SPECIAL CASE OF MERCATOR GRID
       ELSEIF(igdtnumo.EQ.10) THEN
          ni=igdtmplo(8)
          nj=igdtmplo(9)
          rlat1=float(igdtmplo(10))*1.0e-6
          rlon1=float(igdtmplo(11))*1.0e-6
          rlon2=float(igdtmplo(15))*1.0e-6
          rlati=float(igdtmplo(13))*1.0e-6
          iscano=mod(igdtmplo(16)/128,2)
          jscano=mod(igdtmplo(16)/64,2)
          nscano=mod(igdtmplo(16)/32,2)
          dy=float(igdtmplo(19))*1.0e-3
          hi=(-1.)**iscano
          hj=(-1.)**(1-jscano)
          CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
          dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
          dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
          IF(nscano.EQ.0) THEN
             CALL sptrunm(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
                  iprime,iskipi,jskipi,mi,mo,0,0,0, &
                  rlat1,rlon1,dlato,dlono,gi,go)
             ispec=1
          ENDIF
       ENDIF
       !  GENERAL SLOW CASE
       IF(ispec.EQ.0) THEN
          CALL sptrung(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
               iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon,gi,go)
       ENDIF
       DO k=1,km
          ibo(k)=0
          DO n=1,no
             lo(n,k)=.true.
          ENDDO
       ENDDO
    ELSE
       DO k=1,km
          ibo(k)=1
          DO n=1,no
             lo(n,k)=.false.
             go(n,k)=0.
          ENDDO
       ENDDO
    ENDIF
  SUBROUTINE polates4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, & …
  END SUBROUTINE polates4_grib2
 
  !> Interpolate scalar fields (spectral).
  !>
  !> This subprogram performs spectral interpolation from any grid to
  !> any grid for scalar fields.  It requires that the input fields be
  !> uniformly global.
  !>
  !> Options allow choices between triangular shape (ipopt(1)=0) and
  !> rhomboidal shape (ipopt(1)=1) which has no default; a second
  !> option is the truncation (ipopt(2)) which defaults to a sensible
  !> truncation for the input grid (if opt(2)=-1).
  !>
  !> @note If the output grid is not found in a special list, then the
  !> transform back to grid is not very fast. This special list
  !> contains global cylindrical grids, polar stereographic grids
  !> centered at the pole and mercator grids.
  !>
  !> Only horizontal interpolation is performed. The grids are defined
  !> by their grid description sections (passed in integer form as
  !> decoded by subprogram w3fi63()).
  !>
  !> The current code recognizes the following projections:
  !> - kgds(1) = 000 equidistant cylindrical
  !> - kgds(1) = 001 mercator cylindrical
  !> - kgds(1) = 003 lambert conformal conical
  !> - kgds(1) = 004 gaussian cylindrical (spectral native)
  !> - kgds(1) = 005 polar stereographic azimuthal
  !> - kgds(1) = 203 rotated equidistant cylindrical (e-stagger)
  !> - kgds(1) = 205 rotated equidistant cylindrical (b-stagger)
  !>
  !> Where kgds could be either input kgdsi or output kgdso. As an
  !> added bonus the number of output grid points and their latitudes
  !> and longitudes are also returned. On the other hand, the output
  !> can be a set of station points if kgdso(1)<0, in which case the
  !> number of points and their latitudes and longitudes must be
  !> input. Output bitmaps will not be created.
  !>        
  !> ### Program History Log
  !> Date | Programmer | Comments
  !> -----|------------|---------
  !> 96-04-10 | Iredell | Initial
  !> 2001-06-18 | Iredell | improve detection of special fast transform
  !> 2015-01-27 | Gayno | replace calls to gdswiz() with new merged version of gdswzd().
  !>
  !> @param[in] ipopt (20) interpolation options ipopt(1)=0 for
  !> triangular, ipopt(1)=1 for rhomboidal; ipopt(2) is truncation
  !> number (defaults to sensible if ipopt(2)=-1).
  !> @param[in] kgdsi (200) input gds parameters as decoded by w3fi63
  !> @param[in] kgdso (200) output gds parameters (kgdso(1)<0 implies random station points)
  !> @param[in] mi skip number between input grid fields if km>1 or
  !> dimension of input grid fields if km=1
  !> @param[in] mo skip number between output grid fields if km>1 or
  !> dimension of output grid fields if km=1
  !> @param[in] km number of fields to interpolate
  !> @param[in] ibi (km) input bitmap flags (must be all 0)
  !> @param[in] gi (mi,km) input fields to interpolate
  !> @param[out] no number of output points (only if kgdso(1)<0)
  !> @param[out] rlat (no) output latitudes in degrees (if kgdso(1)<0)
  !> @param[out] rlon (no) output longitudes in degrees (if kgdso(1)<0)
  !> @param[out] ibo (km) output bitmap flags
  !> @param[out] lo (mo,km) output bitmaps (always output)
  !> @param[out] go (mo,km) output fields interpolated
  !> @param[out] iret return code
  !> - 0 successful interpolation
  !> - 2 unrecognized input grid or no grid overlap
  !> - 3 unrecognized output grid
  !> - 41 invalid nonglobal input grid
  !> - 42 invalid spectral method parameters
  !>
  !> @author Iredell @date 96-04-10
  suBROUTINE polates4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,GI, &
       NO,RLAT,RLON,IBO,LO,GO,IRET)
    INTEGER,          INTENT(IN   ) :: IPOPT(20), KGDSI(200)
    INTEGER,          INTENT(IN   ) :: KGDSO(200), MI, MO
    INTEGER,          INTENT(IN   ) :: IBI(KM), KM
    INTEGER,          INTENT(  OUT) :: IBO(KM), IRET
    !
    LOGICAL*1,        INTENT(  OUT) :: LO(MO,KM)
    !
    REAL,             INTENT(IN   ) :: GI(MI,KM)
    REAL,             INTENT(INOUT) :: RLAT(MO),RLON(MO)
    REAL,             INTENT(  OUT) :: GO(MO,KM)
    !
    REAL,             PARAMETER     :: FILL=-9999.
    REAL,             PARAMETER     :: RERTH=6.3712e6
    REAL,             PARAMETER     :: PI=3.14159265358979
    REAL,             PARAMETER     :: DPR=180./pi
    !
    INTEGER                         :: IDRTI, IDRTO, IG, JG, IM, JM
    INTEGER                         :: IGO, JGO, IMO, JMO
    INTEGER                         :: ISCAN, JSCAN, NSCAN
    INTEGER                         :: ISCANO, JSCANO, NSCANO
    INTEGER                         :: ISKIPI, JSKIPI
    INTEGER                         :: IMAXI, JMAXI, ISPEC
    INTEGER                         :: IP, IPRIME, IPROJ, IROMB, K
    INTEGER                         :: MAXWV, N, NI, NJ, NPS, NO
    !
    REAL                            :: DE, DR, DY
    REAL                            :: DLAT, DLON, DLATO, DLONO
    REAL                            :: GO2(MO,KM), H, HI, HJ
    REAL                            :: ORIENT
    REAL                            :: RLAT1, RLON1, RLAT2, RLON2, RLATI
    REAL                            :: XMESH, XP, YP
    REAL                            :: XPTS(MO), YPTS(MO)
 
    type(grib1_descriptor) :: desc_in, desc_out
    class(ip_grid), allocatable :: grid_in, grid_out
 
    desc_in = init_descriptor(kgdsi)
    desc_out = init_descriptor(kgdso)
 
    call init_grid(grid_in, desc_in)
    call init_grid(grid_out, desc_out)
    
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
    iret=0
    IF(kgdso(1).GE.0) THEN
       CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
       IF(no.EQ.0) iret=3
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  AFFIRM APPROPRIATE INPUT GRID
    !    LAT/LON OR GAUSSIAN
    !    NO BITMAPS
    !    FULL ZONAL COVERAGE
    !    FULL MERIDIONAL COVERAGE
    idrti=kgdsi(1)
    im=kgdsi(2)
    jm=kgdsi(3)
    rlon1=kgdsi(5)*1.e-3
    rlon2=kgdsi(8)*1.e-3
    iscan=mod(kgdsi(11)/128,2)
    jscan=mod(kgdsi(11)/64,2)
    nscan=mod(kgdsi(11)/32,2)
    IF(idrti.NE.0.AND.idrti.NE.4) iret=41
    DO k=1,km
       IF(ibi(k).NE.0) iret=41
    ENDDO
    IF(iret.EQ.0) THEN
       IF(iscan.EQ.0) THEN
          dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
       ELSE
          dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
       ENDIF
       ig=nint(360/abs(dlon))
       iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
       imaxi=ig
       jmaxi=jm
       IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
    ENDIF
    IF(iret.EQ.0.AND.idrti.EQ.0) THEN
       rlat1=kgdsi(4)*1.e-3
       rlat2=kgdsi(7)*1.e-3
       dlat=(rlat2-rlat1)/(jm-1)
       jg=nint(180/abs(dlat))
       IF(jm.EQ.jg) idrti=256
       IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
    ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
       jg=kgdsi(10)*2
       IF(jm.NE.jg) iret=41
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  SET PARAMETERS
    IF(iret.EQ.0) THEN
       iromb=ipopt(1)
       maxwv=ipopt(2)
       IF(maxwv.EQ.-1) THEN
          IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
          IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
          IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
          IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
          IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
          IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
       ENDIF
       IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  INTERPOLATE
    IF(iret.EQ.0) THEN
       IF(nscan.EQ.0) THEN
          iskipi=1
          jskipi=im
       ELSE
          iskipi=jm
          jskipi=1
       ENDIF
       IF(iscan.EQ.1) iskipi=-iskipi
       IF(jscan.EQ.0) jskipi=-jskipi
       ispec=0
       !  SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
       IF((kgdso(1).EQ.0.OR.kgdso(1).EQ.4).AND. &
            mod(kgdso(2),2).EQ.0.AND.kgdso(5).EQ.0.AND.kgdso(11).EQ.0) THEN
          idrto=kgdso(1)
          imo=kgdso(2)
          jmo=kgdso(3)
          rlon2=kgdso(8)*1.e-3
          dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
          igo=nint(360/abs(dlono))
          IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
             rlat1=kgdso(4)*1.e-3
             rlat2=kgdso(7)*1.e-3
             dlat=(rlat2-rlat1)/(jmo-1)
             jgo=nint(180/abs(dlat))
             IF(jmo.EQ.jgo) idrto=256
             IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
          ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
             jgo=kgdso(10)*2
             IF(jmo.EQ.jgo) ispec=1
          ENDIF
          IF(ispec.EQ.1) THEN
             CALL sptrun(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
                  km,iprime,iskipi,jskipi,mi,0,0,mo,0,gi,go)
          ENDIF
          !  SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
       ELSEIF(kgdso(1).EQ.5.AND. &
            kgdso(2).EQ.kgdso(3).AND.mod(kgdso(2),2).EQ.1.AND. &
            kgdso(8).EQ.kgdso(9).AND.kgdso(11).EQ.64) THEN
          nps=kgdso(2)
          rlat1=kgdso(4)*1.e-3
          rlon1=kgdso(5)*1.e-3
          orient=kgdso(7)*1.e-3
          xmesh=kgdso(8)
          iproj=mod(kgdso(10)/128,2)
          ip=(nps+1)/2
          h=(-1.)**iproj
          de=(1.+sin(60./dpr))*rerth
          dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
          xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
          yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
          IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
             IF(iproj.EQ.0) THEN
                CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
                     iprime,iskipi,jskipi,mi,mo,0,0,0, &
                     60.,xmesh,orient,gi,go,go2)
             ELSE
                CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
                     iprime,iskipi,jskipi,mi,mo,0,0,0, &
                     60.,xmesh,orient,gi,go2,go)
             ENDIF
             ispec=1
          ENDIF
          !  SPECIAL CASE OF MERCATOR GRID
       ELSEIF(kgdso(1).EQ.1) THEN
          ni=kgdso(2)
          nj=kgdso(3)
          rlat1=kgdso(4)*1.e-3
          rlon1=kgdso(5)*1.e-3
          rlon2=kgdso(8)*1.e-3
          rlati=kgdso(9)*1.e-3
          iscano=mod(kgdso(11)/128,2)
          jscano=mod(kgdso(11)/64,2)
          nscano=mod(kgdso(11)/32,2)
          dy=kgdso(13)
          hi=(-1.)**iscano
          hj=(-1.)**(1-jscano)
          dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
          dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
          IF(nscano.EQ.0) THEN
             CALL sptrunm(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
                  iprime,iskipi,jskipi,mi,mo,0,0,0, &
                  rlat1,rlon1,dlato,dlono,gi,go)
             ispec=1
          ENDIF
       ENDIF
       !  GENERAL SLOW CASE
       IF(ispec.EQ.0) THEN
          CALL sptrung(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
               iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon,gi,go)
       ENDIF
       DO k=1,km
          ibo(k)=0
          DO n=1,no
             lo(n,k)=.true.
          ENDDO
       ENDDO
    ELSE
       DO k=1,km
          ibo(k)=1
          DO n=1,no
             lo(n,k)=.false.
             go(n,k)=0.
          ENDDO
       ENDDO
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  suBROUTINE polates4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,GI, & …
  END SUBROUTINE polates4_grib1
 
 
  !> Interpolate vector fields (spectral).
  !>
  !> This subprogram performs spectral interpolation from any grid to
  !> any grid for vector fields. It requires that the input fields be
  !> uniformly global. Options allow choices between triangular shape
  !> (ipopt(1)=0) and rhomboidal shape (ipopt(1)=1) which has no
  !> default; a second option is the truncation (ipopt(2)) which
  !> defaults to a sensible truncation for the input grid (if
  !> opt(2)=-1).
  !>
  !> @note If the output grid is not found in a special list, then the
  !> transform back to grid is not very fast.  This special list
  !> contains global cylindrical grids, polar stereographic grids
  !> centered at the pole and mercator grids. Only horizontal
  !> interpolation is performed.
  !>
  !> The input and output grids are defined by their grib 2 grid
  !> definition template as decoded by the ncep g2 library. The code
  !> recognizes the following projections, where "igdtnumi/o" is the
  !> grib 2 grid defintion template number for the input and output
  !> grids, respectively:
  !> - igdtnumi/o=00 equidistant cylindrical
  !> - igdtnumo  =01 rotated equidistant cylindrical. "e" and non-"e" staggered
  !> - igdtnumo  =10 mercator cylindrical
  !> - igdtnumo  =20 polar stereographic azimuthal
  !> - igdtnumo  =30 lambert conformal conical
  !> - igdtnumi/o=40 gaussian cylindrical
  !>
  !> The input and output vectors are rotated so that they are either
  !> resolved relative to the defined grid in the direction of
  !> increasing x and y coordinates or resolved relative to easterly
  !> and northerly directions, as designated by their respective grid
  !> description sections.
  !>
  !> As an added bonus the number of output grid points and their
  !> latitudes and longitudes are also returned along with their
  !> vector rotation parameters.  On the other hand, the output can be
  !> a set of station points if igdtnumo<0, in which case the number
  !> of points and their latitudes and longitudes must be input along
  !> with their vector rotation parameters.
  !>
  !> Output bitmaps will only be created when the output grid extends
  !> outside of the domain of the input grid.  the output field is set
  !> to 0 where the output bitmap is off.
  !>        
  !> ### Program History Log
  !> Date | Programmer | Comments
  !> -----|------------|---------
  !> 96-04-10 | iredell | initial
  !> 2001-06-18 | iredell | improve detection of special fast transform
  !> 2015-01-27 | gayno | replace calls to gdswiz() with new merged routine gdswzd().
  !> 2015-07-13 | gayno | convert to grib 2. replace grib 1 kgds arrays with grib 2 grid definition template arrays.
  !>
  !> @param[in] ipopt (20) interpolation options ipopt(1)=0 for
  !> triangular, ipopt(1)=1 for rhomboidal; ipopt(2) is truncation
  !> number (defaults to sensible if ipopt(2)=-1).
  !> @param[in] igdtnumi grid definition template number - input
  !> grid. Corresponds to the gfld%igdtnum component of the ncep g2
  !> library gridmod data structure:
  !> - 00 equidistant cylindrical
  !> - 01 rotated equidistant cylindrical.  "e" and non-"e" staggered
  !> - 10 mercator cyclindrical
  !> - 20 polar stereographic azimuthal
  !> - 30 lambert conformal conical
  !> - 40 gaussian equidistant cyclindrical
  !> @param[in] igdtmpli (igdtleni) grid definition template array - input
  !> grid. corresponds to the gfld%igdtmpl component of the ncep g2
  !> library gridmod data structure (section 3 info).  see comments in
  !> routine ipolatev for complete definition.
  !> @param[in] igdtleni number of elements of the grid definition
  !> template array - input grid.  corresponds to the gfld%igdtlen
  !> component of the ncep g2 library gridmod data structure.
  !> @param[in] igdtnumo grid definition template number - output
  !> grid. Corresponds to the gfld%igdtnum component of the ncep g2
  !> library gridmod data structure. igdtnumo<0 means interpolate to
  !> random station points. Otherwise, same definition as "igdtnumi".
  !> @param[in] igdtmplo (igdtleno) grid definition template array -
  !> output grid. corresponds to the gfld%igdtmpl component of the
  !> ncep g2 library gridmod data structure (section 3 info).  see
  !> comments in routine ipolatev() for complete definition.
  !> @param[in] igdtleno number of elements of the grid definition
  !> template array - output grid. Corresponds to the gfld%igdtlen
  !> component of the ncep g2 library gridmod data structure.
  !> @param[in] mi skip number between input grid fields if km>1 or
  !> dimension of input grid fields if km=1.
  !> @param[in] mo skip number between output grid fields if km>1 or
  !> dimension of output grid fields if km=1.
  !> @param[in] km number of fields to interpolate
  !> @param[in] ibi (km) input bitmap flags (must be all 0)
  !> @param[in] ui (mi,km) input u-component fields to interpolate
  !> @param[in] vi (mi,km) input v-component fields to interpolate
  !> @param[out] no number of output points (only if igdtnumo>=0)
  !> @param[inout] rlat (mo) output latitudes in degrees (if igdtnumo<0)
  !> @param[inout] rlon (mo) output longitudes in degrees (if igdtnumo<0)
  !> @param[inout] crot (mo) vector rotation cosines (if igdtnumo<0)
  !> @param[inout] srot (mo) vector rotation sines (if igdtnumo<0)
  !> (ugrid=crot*uearth-srot*vearth; vgrid=srot*uearth+crot*vearth)
  !> @param[out] ibo (km) output bitmap flags
  !> @param[out] lo (mo,km) output bitmaps (always output)
  !> @param[out] uo (mo,km) output u-component fields interpolated
  !> @param[out] vo (mo,km) output v-component fields interpolated
  !> @param[out] iret return code
  !> - 0 successful interpolation
  !> - 2 unrecognized input grid or no grid overlap
  !> - 3 unrecognized output grid
  !> - 41 invalid nonglobal input grid
  !> - 42 invalid spectral method parameters
  !>
  !> @author IREDELL @date 96-04-10
  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)
    INTEGER,          INTENT(IN   ) :: IPOPT(20), IBI(KM)
    INTEGER,          INTENT(IN   ) :: KM, MI, MO
    INTEGER,          INTENT(  OUT) :: IRET, IBO(KM), NO
    INTEGER,          INTENT(IN   ) :: IGDTNUMI, IGDTLENI
    INTEGER,          INTENT(IN   ) :: IGDTMPLI(IGDTLENI)
    INTEGER,          INTENT(IN   ) :: IGDTNUMO, IGDTLENO
    INTEGER,          INTENT(IN   ) :: IGDTMPLO(IGDTLENO)
    !
    LOGICAL*1,        INTENT(  OUT) :: LO(MO,KM)
    !
    REAL,             INTENT(IN   ) :: UI(MI,KM),VI(MI,KM)
    REAL,             INTENT(  OUT) :: UO(MO,KM),VO(MO,KM)
    REAL,             INTENT(INOUT) :: RLAT(MO),RLON(MO)
    REAL,             INTENT(  OUT) :: CROT(MO),SROT(MO)
    !
    REAL,                 PARAMETER :: FILL=-9999.
    REAL,                 PARAMETER :: PI=3.14159265358979
    REAL,                 PARAMETER :: DPR=180./pi
    !
    INTEGER                         :: IDRTO, IROMB, ISKIPI, ISPEC
    INTEGER                         :: IDRTI, IMAXI, JMAXI, IM, JM
    INTEGER                         :: IPRIME, IG, IMO, JMO, IGO, JGO
    INTEGER                         :: ISCAN, JSCAN, NSCAN
    INTEGER                         :: ISCANO, JSCANO, NSCANO
    INTEGER                         :: ISCALE, IP, IPROJ, JSKIPI, JG
    INTEGER                         :: K, MAXWV, N, NI, NJ, NPS
    !
    REAL                            :: DLAT, DLON, DLATO, DLONO, DE, DR, DY
    REAL                            :: E2, H, HI, HJ, DUMM(1)
    REAL                            :: ORIENT, RERTH, SLAT
    REAL                            :: RLAT1, RLON1, RLAT2, RLON2, RLATI
    REAL                            :: UROT, VROT, UO2(MO,KM),VO2(MO,KM)
    REAL                            :: XMESH, XP, YP, XPTS(MO),YPTS(MO)
 
    type(grib2_descriptor) :: desc_in, desc_out
    class(ip_grid), allocatable :: grid_in, grid_out
 
    desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
    desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
 
    call init_grid(grid_in, desc_in)
    call init_grid(grid_out, desc_out)
 
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
    iret=0
    IF(igdtnumo.GE.0) THEN
       CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts, &
            rlon,rlat,no,crot,srot)
       IF(no.EQ.0) iret=3
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  AFFIRM APPROPRIATE INPUT GRID
    !    LAT/LON OR GAUSSIAN
    !    NO BITMAPS
    !    FULL ZONAL COVERAGE
    !    FULL MERIDIONAL COVERAGE
    idrti=igdtnumi
    IF(idrti==40) idrti=4
    IF(idrti==0.OR.idrti==4)THEN
       im=igdtmpli(8)
       jm=igdtmpli(9)
       iscale=igdtmpli(10)*igdtmpli(11)
       IF(iscale==0) iscale=10**6
       rlon1=float(igdtmpli(13))/float(iscale)
       rlon2=float(igdtmpli(16))/float(iscale)
       iscan=mod(igdtmpli(19)/128,2)
       jscan=mod(igdtmpli(19)/64,2)
       nscan=mod(igdtmpli(19)/32,2)
    ELSE
       iret=41
    ENDIF
    DO k=1,km
       IF(ibi(k).NE.0) iret=41
    ENDDO
    IF(iret.EQ.0) THEN
       IF(iscan.EQ.0) THEN
          dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
       ELSE
          dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
       ENDIF
       ig=nint(360/abs(dlon))
       iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
       imaxi=ig
       jmaxi=jm
       IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
    ENDIF
    IF(iret.EQ.0.AND.idrti.EQ.0) THEN
       iscale=igdtmpli(10)*igdtmpli(11)
       IF(iscale==0) iscale=10**6
       rlat1=float(igdtmpli(12))/float(iscale)
       rlat2=float(igdtmpli(15))/float(iscale)
       dlat=(rlat2-rlat1)/(jm-1)
       jg=nint(180/abs(dlat))
       IF(jm.EQ.jg) idrti=256
       IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
    ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
       jg=igdtmpli(18)*2
       IF(jm.NE.jg) iret=41
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  SET PARAMETERS
    IF(iret.EQ.0) THEN
       iromb=ipopt(1)
       maxwv=ipopt(2)
       IF(maxwv.EQ.-1) THEN
          IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
          IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
          IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
          IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
          IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
          IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
       ENDIF
       IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  INTERPOLATE
    IF(iret.EQ.0) THEN
       IF(nscan.EQ.0) THEN
          iskipi=1
          jskipi=im
       ELSE
          iskipi=jm
          jskipi=1
       ENDIF
       IF(iscan.EQ.1) iskipi=-iskipi
       IF(jscan.EQ.0) jskipi=-jskipi
       ispec=0
       !  SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
       IF((igdtnumo.EQ.0.OR.igdtnumo.EQ.40).AND. &
            mod(igdtmplo(8),2).EQ.0.AND.igdtmplo(13).EQ.0.AND. &
            igdtmplo(19).EQ.0) THEN
          idrto=igdtnumo
          IF(idrto==40)idrto=4
          imo=igdtmplo(8)
          jmo=igdtmplo(9)
          iscale=igdtmplo(10)*igdtmplo(11)
          IF(iscale==0) iscale=10**6
          rlon2=float(igdtmplo(16))/float(iscale)
          dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
          igo=nint(360/abs(dlono))
          IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
             rlat1=float(igdtmplo(12))/float(iscale)
             rlat2=float(igdtmplo(15))/float(iscale)
             dlat=(rlat2-rlat1)/(jmo-1)
             jgo=nint(180/abs(dlat))
             IF(jmo.EQ.jgo) idrto=256
             IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
          ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
             jgo=igdtmplo(18)*2
             IF(jmo.EQ.jgo) ispec=1
          ENDIF
          IF(ispec.EQ.1) THEN
             CALL sptrunv(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
                  km,iprime,iskipi,jskipi,mi,0,0,mo,0,ui,vi, &
                  .true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
          ENDIF
          !  SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
       ELSEIF(igdtnumo.EQ.20.AND. &
            igdtmplo(8).EQ.igdtmplo(9).AND.mod(igdtmplo(8),2).EQ.1.AND. &
            igdtmplo(15).EQ.igdtmplo(16).AND.igdtmplo(18).EQ.64.AND. &
            mod(igdtmplo(12)/8,2).EQ.1) THEN
          nps=igdtmplo(8)
          rlat1=float(igdtmplo(10))*1.e-6
          rlon1=float(igdtmplo(11))*1.e-6
          orient=float(igdtmplo(14))*1.e-6
          xmesh=float(igdtmplo(15))*1.e-3
          iproj=mod(igdtmplo(17)/128,2)
          ip=(nps+1)/2
          h=(-1.)**iproj
          slat=float(abs(igdtmplo(13)))*1.e-6
          CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
          de=(1.+sin(slat/dpr))*rerth
          dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
          xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
          yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
          IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
             IF(iproj.EQ.0) THEN
                CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
                     iprime,iskipi,jskipi,mi,mo,0,0,0, &
                     slat,xmesh,orient,ui,vi,.true.,uo,vo,uo2,vo2, &
                     .false.,dumm,dumm,dumm,dumm, &
                     .false.,dumm,dumm,dumm,dumm)
             ELSE
                CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
                     iprime,iskipi,jskipi,mi,mo,0,0,0, &
                     slat,xmesh,orient,ui,vi,.true.,uo2,vo2,uo,vo, &
                     .false.,dumm,dumm,dumm,dumm, &
                     .false.,dumm,dumm,dumm,dumm)
             ENDIF
             ispec=1
          ENDIF
          !  SPECIAL CASE OF MERCATOR GRID
       ELSEIF(igdtnumo.EQ.10) THEN
          ni=igdtmplo(8)
          nj=igdtmplo(9)
          rlat1=float(igdtmplo(10))*1.0e-6
          rlon1=float(igdtmplo(11))*1.0e-6
          rlon2=float(igdtmplo(15))*1.0e-6
          rlati=float(igdtmplo(13))*1.0e-6
          iscano=mod(igdtmplo(16)/128,2)
          jscano=mod(igdtmplo(16)/64,2)
          nscano=mod(igdtmplo(16)/32,2)
          dy=float(igdtmplo(19))*1.0e-3
          hi=(-1.)**iscano
          hj=(-1.)**(1-jscano)
          CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
          dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
          dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
          IF(nscano.EQ.0) THEN
             CALL sptrunmv(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
                  iprime,iskipi,jskipi,mi,mo,0,0,0, &
                  rlat1,rlon1,dlato,dlono,ui,vi, &
                  .true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
             ispec=1
          ENDIF
       ENDIF
       !  GENERAL SLOW CASE
       IF(ispec.EQ.0) THEN
          CALL sptrungv(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
               iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon, &
               ui,vi,.true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
          DO k=1,km
             ibo(k)=0
             DO n=1,no
                lo(n,k)=.true.
                urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
                vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
                uo(n,k)=urot
                vo(n,k)=vrot
             ENDDO
          ENDDO
       ENDIF
    ELSE
       DO k=1,km
          ibo(k)=1
          DO n=1,no
             lo(n,k)=.false.
             uo(n,k)=0.
             vo(n,k)=0.
          ENDDO
       ENDDO
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  SUBROUTINE polatev4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, & …
  END SUBROUTINE polatev4_grib2
 
  !> Interpolate vector fields (spectral).
  !>
  !> This subprogram performs spectral interpolation from any grid to
  !> any grid for vector fields. It requires that the input fields be
  !> uniformly global. Options allow choices between triangular shape
  !> (ipopt(1)=0) and rhomboidal shape (ipopt(1)=1) which has no
  !> default; a second option is the truncation (ipopt(2)) which
  !> defaults to a sensible truncation for the input grid (if
  !> opt(2)=-1).
  !>
  !> @note If the output grid is not found in a special list, then the
  !> transform back to grid is not very fast.  This special list
  !> contains global cylindrical grids, polar stereographic grids
  !> centered at the pole and mercator grids.
  !>
  !> Only horizontal interpolation is performed. The grids are defined
  !> by their grid description sections (passed in integer form as
  !> decoded by subprogram w3fi63).
  !>
  !> The current code recognizes the following projections:
  !> - (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
  !> - (KGDS(1)=001) MERCATOR CYLINDRICAL
  !> - (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
  !> - (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
  !> - (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
  !> - (KGDS(1)=203) ROTATED EQUIDISTANT CYLINDRICAL (E-STAGGER)
  !> - (KGDS(1)=205) ROTATED EQUIDISTANT CYLINDRICAL (B-STAGGER)
  !>
  !> Where kgds could be either input kgdsi or output kgdso.
  !>
  !> The input and output vectors are rotated so that they are either
  !> resolved relative to the defined grid in the direction of
  !> increasing x and y coordinates or resolved relative to easterly
  !> and northerly directions, as designated by their respective grid
  !> description sections. As an added bonus the number of output grid
  !> points and their latitudes and longitudes are also returned along
  !> with their vector rotation parameters. On the other hand, the
  !> output can be a set of station points if kgdso(1)<0, in which
  !> case the number of points and their latitudes and longitudes must
  !> be input along with their vector rotation parameters.
  !>
  !> Output bitmaps will only be created when the output grid extends
  !> outside of the domain of the input grid. The output field is set
  !> to 0 where the output bitmap is off.
  !>        
  !> ### Program History Log
  !> Date | Programmer | Comments
  !> -----|------------|---------
  !> 96-04-10 | iredell | initial.
  !> 2001-06-18 | iredell |  improve detection of special fast transform
  !> 2015-01-27 | gayno | replace calls to gdswiz() with new merged routine gdswzd().
  !>
  !> @param[in] ipopt (20) interpolation options ipopt(1)=0 for
  !> triangular, ipopt(1)=1 for rhomboidal; ipopt(2) is truncation
  !> number (defaults to sensible if ipopt(2)=-1).
  !> @param[in] kgdsi (200) input gds parameters as decoded by w3fi63.
  !> @param[in] kgdso (200) output gds parameters (kgdso(1)<0 implies
  !> random station points).
  !> @param[in] mi skip number between input grid fields if km>1 or
  !> dimension of input grid fields if km=1.
  !> @param[in] mo skip number between output grid fields if km>1 or
  !> dimension of output grid fields if km=1.
  !> @param[in] km number of fields to interpolate
  !> @param[in] ibi (km) input bitmap flags (must be all 0)
  !> @param[in] ui (mi,km) input u-component fields to interpolate
  !> @param[in] vi (mi,km) input v-component fields to interpolate
  !> @param[out] no number of output points (only if kgdso(1)<0)
  !> @param[out] rlat (no) output latitudes in degrees (if kgdso(1)<0)
  !> @param[out] rlon (no) output longitudes in degrees (if kgdso(1)<0)
  !> @param[out] crot (no) vector rotation cosines (if kgdso(1)<0)
  !> @param[out] srot (no) vector rotation sines (if kgdso(1)<0)
  !> (ugrid=crot*uearth-srot*vearth; vgrid=srot*uearth+crot*vearth)
  !> @param[out] ibo (km) output bitmap flags
  !> @param[out] lo (mo,km) output bitmaps (always output)
  !> @param[out] uo (mo,km) output u-component fields interpolated
  !> @param[out] vo (mo,km) output v-component fields interpolated
  !> @param[out] iret return code
  !> - 0 successful interpolation
  !> - 2 unrecognized input grid or no grid overlap
  !> - 3 unrecognized output grid
  !> - 41 invalid nonglobal input grid
  !> - 42 invalid spectral method parameters
  !>
  !> @author IREDELL @date 96-04-10
  SUBROUTINE polatev4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,UI,VI, &
       NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
    INTEGER,          INTENT(IN   ) :: IPOPT(20), IBI(KM)
    INTEGER,          INTENT(IN   ) :: KM, MI, MO
    INTEGER,          INTENT(  OUT) :: IRET, IBO(KM)
    INTEGER,          INTENT(IN) :: KGDSI(200),KGDSO(200)
    !
    LOGICAL*1,        INTENT(  OUT) :: LO(MO,KM)
    !
    REAL,             INTENT(IN   ) :: UI(MI,KM),VI(MI,KM)
    REAL,             INTENT(  OUT) :: UO(MO,KM),VO(MO,KM)
    REAL,             INTENT(INOUT) :: RLAT(MO),RLON(MO)
    REAL,             INTENT(  OUT) :: CROT(MO),SROT(MO)
    !
    REAL,                 PARAMETER :: FILL=-9999.
    REAL,                 PARAMETER :: RERTH=6.3712e6
    REAL,                 PARAMETER :: PI=3.14159265358979
    REAL,                 PARAMETER :: DPR=180./pi
    !
    INTEGER                         :: IDRTO, IROMB, ISKIPI, ISPEC
    INTEGER                         :: IDRTI, IMAXI, JMAXI, IM, JM
    INTEGER                         :: IPRIME, IG, IMO, JMO, IGO, JGO
    INTEGER                         :: ISCAN, JSCAN, NSCAN
    INTEGER                         :: ISCANO, JSCANO, NSCANO
    INTEGER                         :: IP, IPROJ, JSKIPI, JG
    INTEGER                         :: K, MAXWV, N, NI, NJ, NO, NPS
    !
    REAL                            :: DLAT, DLON, DLATO, DLONO, DE, DR, DY
    REAL                            :: H, HI, HJ, DUMM(1)
    REAL                            :: ORIENT
    REAL                            :: RLAT1, RLON1, RLAT2, RLON2, RLATI
    REAL                            :: UROT, VROT, UO2(MO,KM),VO2(MO,KM)
    REAL                            :: XMESH, XP, YP, XPTS(MO),YPTS(MO)
 
    type(grib1_descriptor) :: desc_in, desc_out
    class(ip_grid), allocatable :: grid_in, grid_out
 
    desc_in = init_descriptor(kgdsi)
    desc_out = init_descriptor(kgdso)
 
    call init_grid(grid_in, desc_in)
    call init_grid(grid_out, desc_out)
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
    iret=0
    IF(kgdso(1).GE.0) THEN
       CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
       IF(no.EQ.0) iret=3
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  AFFIRM APPROPRIATE INPUT GRID
    !    LAT/LON OR GAUSSIAN
    !    NO BITMAPS
    !    FULL ZONAL COVERAGE
    !    FULL MERIDIONAL COVERAGE
    idrti=kgdsi(1)
    im=kgdsi(2)
    jm=kgdsi(3)
    rlon1=kgdsi(5)*1.e-3
    rlon2=kgdsi(8)*1.e-3
    iscan=mod(kgdsi(11)/128,2)
    jscan=mod(kgdsi(11)/64,2)
    nscan=mod(kgdsi(11)/32,2)
    IF(idrti.NE.0.AND.idrti.NE.4) iret=41
    DO k=1,km
       IF(ibi(k).NE.0) iret=41
    ENDDO
    IF(iret.EQ.0) THEN
       IF(iscan.EQ.0) THEN
          dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
       ELSE
          dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
       ENDIF
       ig=nint(360/abs(dlon))
       iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
       imaxi=ig
       jmaxi=jm
       IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
    ENDIF
    IF(iret.EQ.0.AND.idrti.EQ.0) THEN
       rlat1=kgdsi(4)*1.e-3
       rlat2=kgdsi(7)*1.e-3
       dlat=(rlat2-rlat1)/(jm-1)
       jg=nint(180/abs(dlat))
       IF(jm.EQ.jg) idrti=256
       IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
    ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
       jg=kgdsi(10)*2
       IF(jm.NE.jg) iret=41
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  SET PARAMETERS
    IF(iret.EQ.0) THEN
       iromb=ipopt(1)
       maxwv=ipopt(2)
       IF(maxwv.EQ.-1) THEN
          IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
          IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
          IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
          IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
          IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
          IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
       ENDIF
       IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  INTERPOLATE
    IF(iret.EQ.0) THEN
       IF(nscan.EQ.0) THEN
          iskipi=1
          jskipi=im
       ELSE
          iskipi=jm
          jskipi=1
       ENDIF
       IF(iscan.EQ.1) iskipi=-iskipi
       IF(jscan.EQ.0) jskipi=-jskipi
       ispec=0
       !  SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
       IF((kgdso(1).EQ.0.OR.kgdso(1).EQ.4).AND. &
            mod(kgdso(2),2).EQ.0.AND.kgdso(5).EQ.0.AND. &
            kgdso(11).EQ.0) THEN
          idrto=kgdso(1)
          imo=kgdso(2)
          jmo=kgdso(3)
          rlon2=kgdso(8)*1.e-3
          dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
          igo=nint(360/abs(dlono))
          IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
             rlat1=kgdso(4)*1.e-3
             rlat2=kgdso(7)*1.e-3
             dlat=(rlat2-rlat1)/(jmo-1)
             jgo=nint(180/abs(dlat))
             IF(jmo.EQ.jgo) idrto=256
             IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
          ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
             jgo=kgdso(10)*2
             IF(jmo.EQ.jgo) ispec=1
          ENDIF
          IF(ispec.EQ.1) THEN
             CALL sptrunv(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
                  km,iprime,iskipi,jskipi,mi,0,0,mo,0,ui,vi, &
                  .true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
          ENDIF
          !  SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
       ELSEIF(kgdso(1).EQ.5.AND. &
            kgdso(2).EQ.kgdso(3).AND.mod(kgdso(2),2).EQ.1.AND. &
            kgdso(8).EQ.kgdso(9).AND.kgdso(11).EQ.64.AND. &
            mod(kgdso(6)/8,2).EQ.1) THEN
          nps=kgdso(2)
          rlat1=kgdso(4)*1.e-3
          rlon1=kgdso(5)*1.e-3
          orient=kgdso(7)*1.e-3
          xmesh=kgdso(8)
          iproj=mod(kgdso(10)/128,2)
          ip=(nps+1)/2
          h=(-1.)**iproj
          de=(1.+sin(60./dpr))*rerth
          dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
          xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
          yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
          IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
             IF(iproj.EQ.0) THEN
                CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
                     iprime,iskipi,jskipi,mi,mo,0,0,0, &
                     60.,xmesh,orient,ui,vi,.true.,uo,vo,uo2,vo2, &
                     .false.,dumm,dumm,dumm,dumm, &
                     .false.,dumm,dumm,dumm,dumm)
             ELSE
                CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
                     iprime,iskipi,jskipi,mi,mo,0,0,0, &
                     60.,xmesh,orient,ui,vi,.true.,uo2,vo2,uo,vo, &
                     .false.,dumm,dumm,dumm,dumm, &
                     .false.,dumm,dumm,dumm,dumm)
             ENDIF
             ispec=1
          ENDIF
          !  SPECIAL CASE OF MERCATOR GRID
       ELSEIF(kgdso(1).EQ.1) THEN
          ni=kgdso(2)
          nj=kgdso(3)
          rlat1=kgdso(4)*1.e-3
          rlon1=kgdso(5)*1.e-3
          rlon2=kgdso(8)*1.e-3
          rlati=kgdso(9)*1.e-3
          iscano=mod(kgdso(11)/128,2)
          jscano=mod(kgdso(11)/64,2)
          nscano=mod(kgdso(11)/32,2)
          dy=kgdso(13)
          hi=(-1.)**iscano
          hj=(-1.)**(1-jscano)
          dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
          dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
          IF(nscano.EQ.0) THEN
             CALL sptrunmv(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
                  iprime,iskipi,jskipi,mi,mo,0,0,0, &
                  rlat1,rlon1,dlato,dlono,ui,vi, &
                  .true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
             ispec=1
          ENDIF
       ENDIF
       !  GENERAL SLOW CASE
       IF(ispec.EQ.0) THEN
          CALL sptrungv(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
               iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon, &
               ui,vi,.true.,uo,vo,.false.,dumm,dumm,.false.,dumm,dumm)
          DO k=1,km
             ibo(k)=0
             DO n=1,no
                lo(n,k)=.true.
                urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
                vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
                uo(n,k)=urot
                vo(n,k)=vrot
             ENDDO
          ENDDO
       ENDIF
    ELSE
       DO k=1,km
          ibo(k)=1
          DO n=1,no
             lo(n,k)=.false.
             uo(n,k)=0.
             vo(n,k)=0.
          ENDDO
       ENDDO
    ENDIF
  SUBROUTINE polatev4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,UI,VI, & …
  END SUBROUTINE polatev4_grib1
end module spectral_interp_mod