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