source: trunk/WRF.COMMON/WRFV2/external/io_grib2/g2lib/gdt2gds.F

Last change on this file was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 10.9 KB
Line 
1        subroutine gdt2gds(igds,igdstmpl,idefnum,ideflist,kgds,
2     &                     igrid,iret)
3C$$$  SUBPROGRAM DOCUMENTATION BLOCK
4C                .      .    .                                       .
5C SUBPROGRAM:    gdt2gds
6C   PRGMMR: Gilbert        ORG: W/NP11    DATE: 2003-06-17
7C
8C ABSTRACT: This routine converts grid information from a GRIB2
9C   Grid Description Section as well as its
10C   Grid Definition Template to GRIB1 GDS info.  In addition,
11C   a check is made to determine if the grid is an NCEP
12C   predefined grid.
13C
14C PROGRAM HISTORY LOG:
15C 2003-06-17  Gilbert
16C 2004-04-27  Gilbert - Added support for gaussian grids.
17C
18C USAGE:    CALL gdt2gds(igds,igdstmpl,idefnum,ideflist,kgds,igrid,iret)
19C   INPUT ARGUMENT LIST:
20C     igds()   - Contains information read from the appropriate GRIB Grid
21C                Definition Section 3 for the field being returned.
22C                Must be dimensioned >= 5.
23C                igds(1)=Source of grid definition (see Code Table 3.0)
24C                igds(2)=Number of grid points in the defined grid.
25C                igds(3)=Number of octets needed for each
26C                            additional grid points definition.
27C                            Used to define number of
28C                            points in each row ( or column ) for
29C                            non-regular grids.
30C                            = 0, if using regular grid.
31C                igds(4)=Interpretation of list for optional points
32C                            definition.  (Code Table 3.11)
33C                igds(5)=Grid Definition Template Number (Code Table 3.1)
34C     igdstmpl() - Grid Definition Template values for GDT 3.igds(5)
35C     idefnum    - The number of entries in array ideflist. 
36C                  i.e. number of rows ( or columns )
37C                  for which optional grid points are defined.
38C     ideflist() - Optional integer array containing
39C                  the number of grid points contained in each row (or column).
40C
41C   OUTPUT ARGUMENT LIST:      (INCLUDING WORK ARRAYS)
42C     kgds()   - GRIB1 GDS as described in w3fi63 format.
43C     igrid    - NCEP predifined GRIB1 grid number
44C                set to 255, if not NCEP grid
45C     iret     - Error return value:
46C                  0  = Successful
47C                  1  = Unrecognized GRIB2 GDT number 3.igds(5)
48C
49C REMARKS: LIST CAVEATS, OTHER HELPFUL HINTS OR INFORMATION
50C
51C ATTRIBUTES:
52C   LANGUAGE: INDICATE EXTENSIONS, COMPILER OPTIONS
53C   MACHINE:  IBM SP
54C
55C$$$
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
Note: See TracBrowser for help on using the repository browser.