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