NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
ip_grid_descriptor_mod.F90
Go to the documentation of this file.
1
5
16 implicit none
17
18 private
19
20 public :: ip_grid_descriptor
23
24 public :: operator(==)
25
28 type, abstract :: ip_grid_descriptor
29 integer :: grid_num
30 contains
33 procedure :: is_same_grid
34 end type ip_grid_descriptor
35
40 integer :: gds(200)
41 contains
44 procedure :: is_same_grid_grib1
45 end type grib1_descriptor
46
50 integer :: gdt_num
51 integer :: gdt_len
52 integer, allocatable :: gdt_tmpl(:)
53 contains
56 procedure :: is_same_grid_grib2
57 end type grib2_descriptor
58
59 interface operator (==)
60 module procedure is_same_grid
61 end interface operator (==)
62
64 module procedure init_grib1_descriptor
65 module procedure init_grib2_descriptor
66 end interface init_descriptor
67
68contains
69
77 function init_grib1_descriptor(gds) result(desc)
78 type(grib1_descriptor) :: desc
79 integer, intent(in) :: gds(:)
80 desc%gds = gds
81 desc%grid_num = gds(1)
82
83 !call desc%decode_template()
84
85 end function init_grib1_descriptor
86
96 function init_grib2_descriptor(gdt_num, gdt_len, gdt_tmpl) result(desc)
97 type(grib2_descriptor) :: desc
98 integer, intent(in) :: gdt_num, gdt_len, gdt_tmpl(:)
99 desc%grid_num = gdt_num
100
101 desc%gdt_num = gdt_num
102 desc%gdt_len = gdt_len
103 allocate(desc%gdt_tmpl(gdt_len))
104 desc%gdt_tmpl = gdt_tmpl
105
106 !call desc%decode_template()
107
108 end function init_grib2_descriptor
109
118 logical function is_same_grid(grid1, grid2)
119 class(ip_grid_descriptor), intent(in) :: grid1, grid2
120
121 select type(grid1)
122 type is(grib1_descriptor)
123 select type(grid2)
124 type is(grib1_descriptor)
125 is_same_grid = grid1%is_same_grid_grib1(grid2)
126 class default
127 is_same_grid = .false.
128 end select
129 type is(grib2_descriptor)
130 select type(grid2)
131 type is(grib2_descriptor)
132 is_same_grid = grid1%is_same_grid_grib2(grid2)
133 class default
134 is_same_grid = .false.
135 end select
136 end select
137
138 end function is_same_grid
139
148 logical function is_same_grid_grib1(self, grid_desc) result(same_grid)
149 class(grib1_descriptor), intent(in) :: self, grid_desc
150
151 if (all(self%gds == grid_desc%gds)) then
152 same_grid = .true.
153 else
154 same_grid = .false.
155 end if
156
157 end function is_same_grid_grib1
158
167 logical function is_same_grid_grib2(self, grid_desc) result(same_grid)
168 class(grib2_descriptor), intent(in) :: self, grid_desc
169
170 same_grid = .false.
171 if (self%grid_num == grid_desc%grid_num) then
172 if (self%gdt_len == grid_desc%gdt_len) then
173 if (all(self%gdt_tmpl == grid_desc%gdt_tmpl)) then
174 same_grid = .true.
175 end if
176 end if
177 end if
178
179 end function is_same_grid_grib2
180
181
182
183 ! subroutine decode_template_grib1(self)
184 ! type(grib1_descriptor), intent(inout) :: self
185
186 ! integer :: im, jm, iwrap, jg
187 ! integer :: iscan, kscan, nscan, nscan_field_pos
188 ! integer :: jwrap1, jwrap2
189
190 ! real :: dlat, dlon
191 ! real :: rlat1, rlat2
192 ! real :: rlon1, rlon2
193
194 ! ! set a default value for this to check if it has changed
195 ! ! for rotated grids nscan is set = 3, but in other places it's not, so use this special value
196 ! ! just for field_position routine
197 ! nscan_field_pos = -1
198
199 ! im = self%gds(2)
200 ! jm = self%gds(3)
201 ! iwrap = 0
202 ! jwrap1 = 0
203 ! jwrap2 = 0
204 ! nscan = mod(self%gds(11) / 32, 2)
205 ! kscan = 0
206
207 ! select case(self%gds(1))
208 ! case(0)
209 ! rlon1=self%gds(5)*1.e-3
210 ! rlon2=self%gds(8)*1.e-3
211 ! iscan=mod(self%gds(11)/128,2)
212 ! if(iscan.eq.0) then
213 ! dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
214 ! else
215 ! dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
216 ! endif
217 ! iwrap=nint(360/abs(dlon))
218 ! if(im.lt.iwrap) iwrap=0
219 ! if(iwrap.gt.0.and.mod(iwrap,2).eq.0) then
220 ! rlat1=self%gds(4)*1.e-3
221 ! rlat2=self%gds(7)*1.e-3
222 ! dlat=abs(rlat2-rlat1)/(jm-1)
223 ! if(abs(rlat1).gt.90-0.25*dlat) then
224 ! jwrap1=2
225 ! elseif(abs(rlat1).gt.90-0.75*dlat) then
226 ! jwrap1=1
227 ! endif
228 ! if(abs(rlat2).gt.90-0.25*dlat) then
229 ! jwrap2=2*jm
230 ! elseif(abs(rlat2).gt.90-0.75*dlat) then
231 ! jwrap2=2*jm+1
232 ! endif
233 ! endif
234 ! case(1)
235 ! rlon1=self%gds(5)*1.e-3
236 ! rlon2=self%gds(8)*1.e-3
237 ! iscan=mod(self%gds(11)/128,2)
238 ! if(iscan.eq.0) then
239 ! dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
240 ! else
241 ! dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
242 ! endif
243 ! iwrap=nint(360/abs(dlon))
244 ! if(im.lt.iwrap) iwrap=0
245 ! case(4)
246 ! rlon1=self%gds(5)*1.e-3
247 ! rlon2=self%gds(8)*1.e-3
248 ! iscan=mod(self%gds(11)/128,2)
249 ! if(iscan.eq.0) then
250 ! dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
251 ! else
252 ! dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
253 ! endif
254 ! iwrap=nint(360/abs(dlon))
255 ! if(im.lt.iwrap) iwrap=0
256 ! if(iwrap.gt.0.and.mod(iwrap,2).eq.0) then
257 ! jg=self%gds(10)*2
258 ! if(jm.eq.jg) then
259 ! jwrap1=1
260 ! jwrap2=2*jm+1
261 ! endif
262 ! endif
263 ! case(203)
264 ! ! why is nscan set to 3 for this grid?
265 ! nscan_field_pos = 3
266 ! kscan=mod(self%gds(11) / 256, 2)
267 ! case default
268 ! print *, "grib1 grid type ", self%gds(1), " not recognized"
269 ! error stop
270 ! end select
271
272 ! self%im = im
273 ! self%jm = jm
274 ! self%nm = im * jm
275 ! self%iwrap = iwrap
276 ! self%jwrap1 = jwrap1
277 ! self%jwrap2 = jwrap2
278 ! self%nscan = nscan
279
280 ! if (nscan_field_pos == -1) then
281 ! ! just use regular value of nscan
282 ! self%nscan_field_pos = nscan
283 ! else
284 ! ! nscan = 3 for rotated grids and is set to 3 for use in field_position
285 ! self%nscan_field_pos = nscan_field_pos
286 ! end if
287
288 ! self%kscan = kscan
289
290 ! end subroutine decode_template_grib1
291
292 ! subroutine decode_template_grib2(self)
293 ! type(grib2_descriptor), intent(inout) :: self
294
295 ! integer :: im, jm, iwrap, jg
296 ! integer :: i_offset_odd, i_offset_even
297 ! integer :: iscan, kscan, nscan, nscan_field_pos
298 ! integer :: jwrap1, jwrap2, iscale
299
300 ! real :: dlat, dlon
301 ! real :: rlat1, rlat2
302 ! real :: rlon1, rlon2
303
304 ! nscan_field_pos = -1
305
306 ! select case(self%gdt_num)
307 ! ! EQUIDISTANT CYLINDRICAL
308 ! case(0)
309 ! im=self%gdt_tmpl(8)
310 ! jm=self%gdt_tmpl(9)
311 ! nscan=mod(self%gdt_tmpl(19)/32,2)
312 ! kscan=0
313 ! iscale=self%gdt_tmpl(10)*self%gdt_tmpl(11)
314 ! if(iscale==0) iscale=10**6
315 ! rlon1=float(self%gdt_tmpl(13))/float(iscale)
316 ! rlon2=float(self%gdt_tmpl(16))/float(iscale)
317 ! iscan=mod(self%gdt_tmpl(19)/128,2)
318 ! if(iscan.eq.0) then
319 ! dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
320 ! else
321 ! dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
322 ! endif
323 ! iwrap=nint(360/abs(dlon))
324 ! if(im.lt.iwrap) iwrap=0
325 ! jwrap1=0
326 ! jwrap2=0
327 ! if(iwrap.gt.0.and.mod(iwrap,2).eq.0) then
328 ! rlat1=float(self%gdt_tmpl(12))/float(iscale)
329 ! rlat2=float(self%gdt_tmpl(15))/float(iscale)
330 ! dlat=abs(rlat2-rlat1)/(jm-1)
331 ! if(abs(rlat1).gt.90-0.25*dlat) then
332 ! jwrap1=2
333 ! elseif(abs(rlat1).gt.90-0.75*dlat) then
334 ! jwrap1=1
335 ! endif
336 ! if(abs(rlat2).gt.90-0.25*dlat) then
337 ! jwrap2=2*jm
338 ! elseif(abs(rlat2).gt.90-0.75*dlat) then
339 ! jwrap2=2*jm+1
340 ! endif
341 ! endif
342 ! case(1)
343 ! i_offset_odd=mod(self%gdt_tmpl(19)/8,2)
344 ! i_offset_even=mod(self%gdt_tmpl(19)/4,2)
345 ! im=self%gdt_tmpl(8)
346 ! jm=self%gdt_tmpl(9)
347 ! iwrap=0
348 ! jwrap1=0
349 ! jwrap2=0
350 ! kscan=0
351 ! nscan=mod(self%gdt_tmpl(19)/32,2)
352 ! if(i_offset_odd/=i_offset_even)then
353 ! kscan=i_offset_odd
354 ! nscan_field_pos=3
355 ! endif
356 ! ! MERCATOR CYLINDRICAL
357 ! case(10)
358 ! im=self%gdt_tmpl(8)
359 ! jm=self%gdt_tmpl(9)
360 ! rlon1=float(self%gdt_tmpl(11))*1.0e-6
361 ! rlon2=float(self%gdt_tmpl(15))*1.0e-6
362 ! iscan=mod(self%gdt_tmpl(16)/128,2)
363 ! if(iscan.eq.0) then
364 ! dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
365 ! else
366 ! dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
367 ! endif
368 ! iwrap=nint(360/abs(dlon))
369 ! if(im.lt.iwrap) iwrap=0
370 ! jwrap1=0
371 ! jwrap2=0
372 ! kscan=0
373 ! nscan=mod(self%gdt_tmpl(16)/32,2)
374 ! ! POLAR STEREOGRAPHIC AZIMUTHAL
375 ! case(20)
376 ! im=self%gdt_tmpl(8)
377 ! jm=self%gdt_tmpl(9)
378 ! nscan=mod(self%gdt_tmpl(18)/32,2)
379 ! iwrap=0
380 ! jwrap1=0
381 ! jwrap2=0
382 ! kscan=0
383 ! ! LAMBERT CONFORMAL CONICAL
384 ! case(30)
385 ! im=self%gdt_tmpl(8)
386 ! jm=self%gdt_tmpl(9)
387 ! nscan=mod(self%gdt_tmpl(18)/32,2)
388 ! iwrap=0
389 ! jwrap1=0
390 ! jwrap2=0
391 ! kscan=0
392 ! ! GAUSSIAN CYLINDRICAL
393 ! case(40)
394 ! im=self%gdt_tmpl(8)
395 ! jm=self%gdt_tmpl(9)
396 ! iscale=self%gdt_tmpl(10)*self%gdt_tmpl(11)
397 ! if(iscale==0) iscale=10**6
398 ! rlon1=float(self%gdt_tmpl(13))/float(iscale)
399 ! rlon2=float(self%gdt_tmpl(16))/float(iscale)
400 ! iscan=mod(self%gdt_tmpl(19)/128,2)
401 ! if(iscan.eq.0) then
402 ! dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
403 ! else
404 ! dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
405 ! endif
406 ! iwrap=nint(360/abs(dlon))
407 ! if(im.lt.iwrap) iwrap=0
408 ! jwrap1=0
409 ! jwrap2=0
410 ! if(iwrap.gt.0.and.mod(iwrap,2).eq.0) then
411 ! jg=self%gdt_tmpl(18)*2
412 ! if(jm.eq.jg) then
413 ! jwrap1=1
414 ! jwrap2=2*jm+1
415 ! endif
416 ! endif
417 ! nscan=mod(self%gdt_tmpl(19)/32,2)
418 ! kscan=0
419 ! case default
420 ! print *, "gdt_num ", self%gdt_num, " not recognized"
421 ! error stop
422 ! end select
423
424 ! self%im = im
425 ! self%jm = jm
426 ! self%nm = im*jm
427 ! self%iwrap = iwrap
428 ! self%jwrap1 = jwrap1
429 ! self%jwrap2 = jwrap2
430 ! self%nscan = nscan
431 ! if (nscan_field_pos == -1) then
432 ! self%nscan_field_pos = nscan
433 ! else
434 ! self%nscan_field_pos = nscan_field_pos
435 ! end if
436 ! self%kscan = kscan
437
438 ! end subroutine decode_template_grib2
439
440 ! subroutine earth_radius(self, radius, eccen_squared)
441 ! class(ip_grid_descriptor), intent(in) :: self
442 ! real, intent(out) :: radius, eccen_squared
443
444 ! real :: flat, major_axis, minor_axis
445 ! logical :: elliptical
446
447 ! select type(self)
448 ! type is(grib1_descriptor)
449 ! associate(gds => self%gds)
450 ! select case(gds(1))
451 ! case(0)
452 ! radius = 6.3712d6
453 ! eccen_squared = 0d0
454 ! case(1)
455 ! radius = 6.3712d6
456 ! eccen_squared = 0d0
457 ! case(3)
458 ! radius = 6.3712d6
459 ! eccen_squared = 0d0
460 ! case(4)
461 ! radius = 6.3712d6
462 ! eccen_squared = 0d0
463 ! case(5)
464 ! elliptical = mod(gds(6) / 64, 2) == 1
465 ! if (.not. elliptical) then
466 ! radius = 6.3712d6
467 ! eccen_squared = 0d0
468 ! else
469 ! radius = 6.378137E6 ! WGS-84
470 ! eccen_squared = 0.00669437999013d0
471 ! end if
472 ! case(203)
473 ! radius = 6.3712d6
474 ! eccen_squared = 0d0
475 ! case(205)
476 ! radius = 6.3712d6
477 ! eccen_squared = 0d0
478 ! case default
479 ! print *, "grib1 grid not recognized. gds(1) = ", gds(1)
480 ! error stop
481 ! end select
482 ! end associate
483 ! type is(grib2_descriptor)
484 ! associate(gdt_tmpl => self%gdt_tmpl)
485 ! select case (gdt_tmpl(1))
486 ! case (0)
487 ! radius = 6367470.0
488 ! eccen_squared = 0.0
489 ! case (1) ! user specified spherical
490 ! radius = float(gdt_tmpl(3))/float(10**gdt_tmpl(2))
491 ! eccen_squared = 0.0
492 ! case (2) ! iau 1965
493 ! radius = 6378160.0 ! semi major axis
494 ! flat = 1.0/297.0 ! flattening
495 ! eccen_squared = (2.0*flat) - (flat**2)
496 ! case (3) ! user specified elliptical (km)
497 ! major_axis = float(gdt_tmpl(5))/float(10**gdt_tmpl(4))
498 ! major_axis = major_axis * 1000.0
499 ! minor_axis = float(gdt_tmpl(7))/float(10**gdt_tmpl(6))
500 ! minor_axis = minor_axis * 1000.0
501 ! eccen_squared = 1.0 - (minor_axis**2 / major_axis**2)
502 ! radius = major_axis
503 ! case (4) ! iag-grs80 model
504 ! radius = 6378137.0 ! semi major axis
505 ! flat = 1.0/298.2572 ! flattening
506 ! eccen_squared = (2.0*flat) - (flat**2)
507 ! case (5) ! wgs84 datum
508 ! radius = 6378137.0 ! semi major axis
509 ! eccen_squared = 0.00669437999013
510 ! case (6)
511 ! radius = 6371229.0
512 ! eccen_squared = 0.0
513 ! case (7) ! user specified elliptical (m)
514 ! major_axis = float(gdt_tmpl(5))/float(10**gdt_tmpl(4))
515 ! minor_axis = float(gdt_tmpl(7))/float(10**gdt_tmpl(6))
516 ! eccen_squared = 1.0 - (minor_axis**2 / major_axis**2)
517 ! radius = major_axis
518 ! case (8)
519 ! radius = 6371200.0
520 ! eccen_squared = 0.0
521 ! case default
522 ! radius = -9999.
523 ! eccen_squared = -9999.
524 ! error stop
525 ! end select
526 ! end associate
527 ! class default
528 ! print *, "unknown descriptor type"
529 ! error stop
530 ! end select
531 ! end subroutine earth_radius
532
533end module ip_grid_descriptor_mod
Users derived type grid descriptor objects to abstract away the raw GRIB1 and GRIB2 grid definitions.
type(grib2_descriptor) function, public init_grib2_descriptor(gdt_num, gdt_len, gdt_tmpl)
Initialize grib-2 descriptor from integer grid definition template (GDT).
logical function is_same_grid(grid1, grid2)
Test whether two grid descriptors are the same.
logical function is_same_grid_grib1(self, grid_desc)
Test whether two grib1_descriptors are the same.
logical function is_same_grid_grib2(self, grid_desc)
Test whether two grib2_descriptors are the same.
type(grib1_descriptor) function, public init_grib1_descriptor(gds)
Initialize grib-1 descriptor from integer grid definition section (GDS).
Descriptor representing a grib1 grib descriptor section (GDS) with an integer array.
Grib-2 descriptor containing a grib2 GDT represented by an integer array.
Abstract descriptor object which represents a grib1 or grib2 descriptor.