WPS/ungrib/src/ngl/g2/gdt2gds.f

   1         subroutine gdt2gds(igds,igdstmpl,idefnum,ideflist,kgds,
   2      &                     igrid,iret)
   3 C$$$  SUBPROGRAM DOCUMENTATION BLOCK
   4 C                .      .    .                                       .
   5 C SUBPROGRAM:    gdt2gds
   6 C   PRGMMR: Gilbert        ORG: W/NP11    DATE: 2003-06-17
   7 C
   8 C ABSTRACT: This routine converts grid information from a GRIB2
   9 C   Grid Description Section as well as its
  10 C   Grid Definition Template to GRIB1 GDS info.  In addition,
  11 C   a check is made to determine if the grid is an NCEP
  12 C   predefined grid.
  13 C
  14 C PROGRAM HISTORY LOG:
  15 C 2003-06-17  Gilbert
  16 C 2004-04-27  Gilbert - Added support for gaussian grids.
  17 C 2007-04-16  Vuong   - Added Curvilinear Orthogonal grids.
  18 C
  19 C USAGE:    CALL gdt2gds(igds,igdstmpl,idefnum,ideflist,kgds,igrid,iret)
  20 C   INPUT ARGUMENT LIST:
  21 C     igds()   - Contains information read from the appropriate GRIB Grid
  22 C                Definition Section 3 for the field being returned.
  23 C                Must be dimensioned >= 5.
  24 C                igds(1)=Source of grid definition (see Code Table 3.0)
  25 C                igds(2)=Number of grid points in the defined grid.
  26 C                igds(3)=Number of octets needed for each
  27 C                            additional grid points definition.
  28 C                            Used to define number of
  29 C                            points in each row ( or column ) for
  30 C                            non-regular grids.
  31 C                            = 0, if using regular grid.
  32 C                igds(4)=Interpretation of list for optional points
  33 C                            definition.  (Code Table 3.11)
  34 C                igds(5)=Grid Definition Template Number (Code Table 3.1)
  35 C     igdstmpl() - Grid Definition Template values for GDT 3.igds(5)
  36 C     idefnum    - The number of entries in array ideflist.
  37 C                  i.e. number of rows ( or columns )
  38 C                  for which optional grid points are defined.
  39 C     ideflist() - Optional integer array containing
  40 C                  the number of grid points contained in each row (or column).
  41 C
  42 C   OUTPUT ARGUMENT LIST:      (INCLUDING WORK ARRAYS)
  43 C     kgds()   - GRIB1 GDS as described in w3fi63 format.
  44 C     igrid    - NCEP predefined GRIB1 grid number
  45 C                set to 255, if not NCEP grid
  46 C     iret     - Error return value:
  47 C                  0  = Successful
  48 C                  1  = Unrecognized GRIB2 GDT number 3.igds(5)
  49 C
  50 C REMARKS: LIST CAVEATS, OTHER HELPFUL HINTS OR INFORMATION
  51 C
  52 C ATTRIBUTES:
  53 C   LANGUAGE: INDICATE EXTENSIONS, COMPILER OPTIONS
  54 C   MACHINE:  IBM SP
  55 C
  56 C$$$
  57 !
  58         integer,intent(in) :: idefnum
  59         integer,intent(in) :: igds(*),igdstmpl(*),ideflist(*)
  60         integer,intent(out) :: kgds(*),igrid,iret
  61
  62         integer :: kgds72(200),kgds71(200),idum(200),jdum(200)
  63
  64         iret=0
  65         if (igds(5).eq.0) then       !  Lat/Lon grid
  66            kgds(1)=0
  67            kgds(2)=igdstmpl(8)            ! Ni
  68            kgds(3)=igdstmpl(9)            ! Nj
  69            kgds(4)=igdstmpl(12)/1000      ! Lat of 1st grid point
  70            kgds(5)=igdstmpl(13)/1000      ! Long of 1st grid point
  71            kgds(6)=0                      ! resolution and component flags
  72            if (igdstmpl(1)==2 ) kgds(6)=64
  73            if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) )
  74      &         kgds(6)=kgds(6)+128
  75            if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
  76            kgds(7)=igdstmpl(15)/1000      ! Lat of last grid point
  77            kgds(8)=igdstmpl(16)/1000      ! Long of last grid point
  78            kgds(9)=igdstmpl(17)/1000      ! Di
  79            kgds(10)=igdstmpl(18)/1000     ! Dj
  80            kgds(11)=igdstmpl(19)          ! Scanning mode
  81            kgds(12)=0
  82            kgds(13)=0
  83            kgds(14)=0
  84            kgds(15)=0
  85            kgds(16)=0
  86            kgds(17)=0
  87            kgds(18)=0
  88            kgds(19)=0
  89            kgds(20)=255
  90            kgds(21)=0
  91            kgds(22)=0
  92            !
  93            !  Process irreg grid stuff, if necessary
  94            !
  95            if ( idefnum.ne.0 ) then
  96               if ( igdstmpl(8).eq.-1 ) then
  97                  kgds(2)=65535
  98                  kgds(9)=65535
  99               endif
 100               if ( igdstmpl(9).eq.-1 ) then
 101                  kgds(3)=65535
 102                  kgds(10)=65535
 103               endif
 104               kgds(19)=0
 105               kgds(20)=33
 106               if ( kgds(1).eq.1.OR.kgds(1).eq.3 ) kgds(20)=43
 107               kgds(21)=igds(2)                   ! num of grid points
 108               do j=1,idefnum
 109                  kgds(21+j)=ideflist(j)
 110               enddo
 111            endif
 112         elseif (igds(5).eq.10) then       !  Mercator grid
 113            kgds(1)=1                 ! Grid Definition Template number
 114            kgds(2)=igdstmpl(8)            ! Ni
 115            kgds(3)=igdstmpl(9)            ! Nj
 116            kgds(4)=igdstmpl(10)/1000      ! Lat of 1st grid point
 117            kgds(5)=igdstmpl(11)/1000      ! Long of 1st grid point
 118            kgds(6)=0                      ! resolution and component flags
 119            if (igdstmpl(1)==2 ) kgds(6)=64
 120            if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) )
 121      &         kgds(6)=kgds(6)+128
 122            if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8
 123            kgds(7)=igdstmpl(14)/1000      ! Lat of last grid point
 124            kgds(8)=igdstmpl(15)/1000      ! Long of last grid point
 125            kgds(9)=igdstmpl(13)/1000      ! Lat intersects earth
 126            kgds(10)=0
 127            kgds(11)=igdstmpl(16)          ! Scanning mode
 128            kgds(12)=igdstmpl(18)/1000     ! Di
 129            kgds(13)=igdstmpl(19)/1000     ! Dj
 130            kgds(14)=0
 131            kgds(15)=0
 132            kgds(16)=0
 133            kgds(17)=0
 134            kgds(18)=0
 135            kgds(19)=0
 136            kgds(20)=255
 137            kgds(21)=0
 138            kgds(22)=0
 139         elseif (igds(5).eq.30) then       ! Lambert Conformal Grid
 140            kgds(1)=3
 141            kgds(2)=igdstmpl(8)            ! Nx
 142            kgds(3)=igdstmpl(9)            ! Ny
 143            kgds(4)=igdstmpl(10)/1000      ! Lat of 1st grid point
 144            kgds(5)=igdstmpl(11)/1000      ! Long of 1st grid point
 145            kgds(6)=0                      ! resolution and component flags
 146            if (igdstmpl(1)==2 ) kgds(6)=64
 147            if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) )
 148      &         kgds(6)=kgds(6)+128
 149            if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8
 150            kgds(7)=igdstmpl(14)/1000      ! Lon of orientation
 151            kgds(8)=igdstmpl(15)/1000      ! Dx
 152            kgds(9)=igdstmpl(16)/1000      ! Dy
 153            kgds(10)=igdstmpl(17)          ! Projection Center Flag
 154            kgds(11)=igdstmpl(18)          ! Scanning mode
 155            kgds(12)=igdstmpl(19)/1000     ! Lat in 1
 156            kgds(13)=igdstmpl(20)/1000     ! Lat in 2
 157            kgds(14)=igdstmpl(21)/1000     ! Lat of S. Pole of projection
 158            kgds(15)=igdstmpl(22)/1000     ! Lon of S. Pole of projection
 159            kgds(16)=0
 160            kgds(17)=0
 161            kgds(18)=0
 162            kgds(19)=0
 163            kgds(20)=255
 164            kgds(21)=0
 165            kgds(22)=0
 166         elseif (igds(5).eq.40) then       !  Gaussian Lat/Lon grid
 167            kgds(1)=4
 168            kgds(2)=igdstmpl(8)            ! Ni
 169            kgds(3)=igdstmpl(9)            ! Nj
 170            kgds(4)=igdstmpl(12)/1000      ! Lat of 1st grid point
 171            kgds(5)=igdstmpl(13)/1000      ! Long of 1st grid point
 172            kgds(6)=0                      ! resolution and component flags
 173            if (igdstmpl(1)==2 ) kgds(6)=64
 174            if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) )
 175      &         kgds(6)=kgds(6)+128
 176            if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
 177            kgds(7)=igdstmpl(15)/1000      ! Lat of last grid point
 178            kgds(8)=igdstmpl(16)/1000      ! Long of last grid point
 179            kgds(9)=igdstmpl(17)/1000      ! Di
 180            kgds(10)=igdstmpl(18)          ! N - Number of parallels
 181            kgds(11)=igdstmpl(19)          ! Scanning mode
 182            kgds(12)=0
 183            kgds(13)=0
 184            kgds(14)=0
 185            kgds(15)=0
 186            kgds(16)=0
 187            kgds(17)=0
 188            kgds(18)=0
 189            kgds(19)=0
 190            kgds(20)=255
 191            kgds(21)=0
 192            kgds(22)=0
 193         elseif (igds(5).eq.20) then       ! Polar Stereographic Grid
 194            kgds(1)=5
 195            kgds(2)=igdstmpl(8)            ! Nx
 196            kgds(3)=igdstmpl(9)            ! Ny
 197            kgds(4)=igdstmpl(10)/1000      ! Lat of 1st grid point
 198            kgds(5)=igdstmpl(11)/1000      ! Long of 1st grid point
 199            kgds(6)=0                      ! resolution and component flags
 200            if (igdstmpl(1)==2 ) kgds(6)=64
 201            if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) )
 202      &         kgds(6)=kgds(6)+128
 203            if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8
 204            kgds(7)=igdstmpl(14)/1000      ! Lon of orientation
 205            kgds(8)=igdstmpl(15)/1000      ! Dx
 206            kgds(9)=igdstmpl(16)/1000      ! Dy
 207            kgds(10)=igdstmpl(17)          ! Projection Center Flag
 208            kgds(11)=igdstmpl(18)          ! Scanning mode
 209            kgds(12)=0
 210            kgds(13)=0
 211            kgds(14)=0
 212            kgds(15)=0
 213            kgds(16)=0
 214            kgds(17)=0
 215            kgds(18)=0
 216            kgds(19)=0
 217            kgds(20)=255
 218            kgds(21)=0
 219            kgds(22)=0
 220         elseif (igds(5).eq.204) then      ! Curvilinear Orthogonal
 221            kgds(1)=204
 222            kgds(2)=igdstmpl(8)            ! Ni
 223            kgds(3)=igdstmpl(9)            ! Nj
 224            kgds(4)=0
 225            kgds(5)=0
 226            kgds(6)=0                      ! resolution and component flags
 227            if (igdstmpl(1)==2 ) kgds(6)=64
 228            if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) )
 229      &         kgds(6)=kgds(6)+128
 230            if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
 231            kgds(7)=0
 232            kgds(8)=0
 233            kgds(9)=0
 234            kgds(10)=0
 235            kgds(11)=igdstmpl(19)          ! Scanning mode
 236            kgds(12)=0
 237            kgds(13)=0
 238            kgds(14)=0
 239            kgds(15)=0
 240            kgds(16)=0
 241            kgds(17)=0
 242            kgds(18)=0
 243            kgds(19)=0
 244            kgds(20)=255
 245            kgds(21)=0
 246            kgds(22)=0
 247            !
 248            !  Process irreg grid stuff, if necessary
 249            !
 250            if ( idefnum.ne.0 ) then
 251               if ( igdstmpl(8).eq.-1 ) then
 252                  kgds(2)=65535
 253                  kgds(9)=65535
 254               endif
 255               if ( igdstmpl(9).eq.-1 ) then
 256                  kgds(3)=65535
 257                  kgds(10)=65535
 258               endif
 259               kgds(19)=0
 260               kgds(20)=33
 261               if ( kgds(1).eq.1.OR.kgds(1).eq.3 ) kgds(20)=43
 262               kgds(21)=igds(2)                   ! num of grid points
 263               do j=1,idefnum
 264                  kgds(21+j)=ideflist(j)
 265               enddo
 266            endif
 267         else
 268            Print *,'gdt2gds: Unrecognized GRIB2 GDT = 3.',igds(5)
 269            iret=1
 270            kgds(1:22)=0
 271            return
 272         endif
 273 !
 274 !   Can we determine NCEP grid number ?
 275 !
 276         igrid=255
 277         do j=254,1,-1
 278         !do j=225,225
 279            kgds71=0
 280            kgds72=0
 281            call w3fi71(j,kgds71,ierr)
 282            if ( ierr.ne.0 ) cycle
 283            ! convert W to E for longitudes
 284            if ( kgds71(3).eq.0 ) then    ! lat/lon
 285               if ( kgds71(7).lt.0 ) kgds71(7)=360000+kgds71(7)
 286               if ( kgds71(10).lt.0 ) kgds71(10)=360000+kgds71(10)
 287            elseif ( kgds71(3).eq.1 ) then    ! mercator
 288               if ( kgds71(7).lt.0 ) kgds71(7)=360000+kgds71(7)
 289               if ( kgds71(10).lt.0 ) kgds71(10)=360000+kgds71(10)
 290            elseif ( kgds71(3).eq.3 ) then     ! lambert conformal
 291               if ( kgds71(7).lt.0 ) kgds71(7)=360000+kgds71(7)
 292               if ( kgds71(9).lt.0 ) kgds71(9)=360000+kgds71(9)
 293               if ( kgds71(18).lt.0 ) kgds71(18)=360000+kgds71(18)
 294            elseif ( kgds71(3).eq.4 ) then     ! Guassian lat/lon
 295               if ( kgds71(7).lt.0 ) kgds71(7)=360000+kgds71(7)
 296               if ( kgds71(10).lt.0 ) kgds71(10)=360000+kgds71(10)
 297            elseif ( kgds71(3).eq.5 ) then     ! polar stereographic
 298               if ( kgds71(7).lt.0 ) kgds71(7)=360000+kgds71(7)
 299               if ( kgds71(9).lt.0 ) kgds71(9)=360000+kgds71(9)
 300            endif
 301            call r63w72(idum,kgds,jdum,kgds72)
 302            if ( kgds72(3).eq.3 ) kgds72(14)=0    ! lambert conformal fix
 303            if ( kgds72(3).eq.1 ) kgds72(15:18)=0    ! mercator fix
 304            if ( kgds72(3).eq.5 ) kgds72(14:18)=0    ! polar str fix
 305 c           print *,' kgds71(',j,')= ', kgds71(1:30)
 306 c           print *,' kgds72       = ', kgds72(1:30)
 307            if ( all(kgds71.eq.kgds72) ) then
 308               igrid=j
 309               return
 310            endif
 311         enddo
 312
 313         return
 314         end