| blob | 1e233ea1c10fa20e9d8b84962623c15b882a9b24 |
1 !*****************************************************************************!
2 program g2print !
3 ! !
4 use table
5 use gridinfo
6 use filelist
7 implicit none
8 interface
9 subroutine parse_args(err, a1, h1, i1, l1, a2, h2, i2, l2,&
10 a3, h3, i3, l3, hlast)
11 integer :: err
12 character(len=*) , optional :: a1, a2, a3
13 character(len=*), optional :: h1, h2, h3
14 integer , optional :: i1, i2, i3
15 logical, optional :: l1, l2, l3
16 character(len=*), optional :: hlast
17 end subroutine parse_args
18 end interface
20 integer :: nunit1 = 12
21 character(LEN=120) :: gribflnm
23 integer :: iprint
25 integer , parameter :: maxlvl = 150
27 real :: startlat, startlon, deltalat, deltalon
28 real :: level
29 character (LEN=9) :: field
30 character (LEN=3) :: out_format
32 logical :: readit
34 integer, dimension(255) :: iuarr = 0
36 character (LEN=19) :: HSTART, HEND, HDATE
37 character(LEN=19) :: hsave = '0000-00-00_00:00:00'
38 integer :: itime
39 integer :: ntimes
40 integer :: interval
41 integer :: ierr
42 logical :: ordered_by_date
43 integer :: debug_level
44 integer :: grib_version
45 integer :: vtable_columns
47 character(len=30) :: hopt
48 logical :: ivb = .FALSE.
49 logical :: idb = .FALSE.
51 ! -----------------
53 gribflnm = ' '
54 call parse_args(ierr, a1='v', l1=ivb, a2='V', l2=idb, hlast=gribflnm)
55 if (ierr.ne.0) then
56 call getarg(0, hopt)
57 write(*,'(//,"Usage: ", A, " [-v] [-V] file",/)') trim(hopt)
58 write(*,'(" -v : Print more information about the GRIB records")')
59 write(*,'(" -V : Print way too much information about the GRIB&
60 & records")')
61 write(*,'(" file : GRIB file to read"//)')
62 stop
63 endif
65 ! -----------------
66 ! Determine GRIB Edition number
67 grib_version=0
68 call edition_num(nunit1, trim(gribflnm), grib_version, ierr)
69 if (ierr.eq.3) STOP 'GRIB file problem'
72 debug_level = 0
73 if (ivb) debug_level = 51
74 if (idb) debug_level = 101
75 write(6,*) 'reading from grib file = ',gribflnm
77 LOOP1 : DO
78 ! At the beginning of LOOP1, we are at a new time period.
79 ! Clear the storage arrays and associated level information.
81 ! If we need to read a new grib record, then read one.
83 if (grib_version.ne.2) then
84 ! write(6,*) 'calling r_grib1 with iunit ', nunit1
85 ! write(6,*) 'flnm = ',gribflnm
86 write(6,*) 'This is a Grib1 file. Please use g1print.\n'
87 stop
88 ! Read one record at a time from GRIB1 (and older Editions)
89 ! call r_grib1(nunit1, gribflnm, level, field, &
90 ! hdate, debug_level, ierr, iuarr, iprint)
91 else
93 ! Read one file of records from GRIB2.
94 if (debug_level .gt. 100) write(6,*) 'calling r_grib2'
95 call r_grib2(nunit1, gribflnm, hdate, &
96 grib_version, debug_level, ierr)
98 endif
100 if (ierr.eq.1) then
101 ! We have hit the end of a file. Exit LOOP1.
102 exit LOOP1
103 endif
105 enddo LOOP1
107 if (grib_version.ne.2) then
108 call c_close(iuarr(nunit1), iprint, ierr)
109 iuarr(nunit1) = 0
110 endif
112 ! And Now we are done:
114 print*,' '
115 print*,' '
116 print*,' Successful completion of g2print '
118 contains
119 subroutine sort_filedates
120 implicit none
122 integer :: n
123 logical :: done
124 if (nfiles > 1) then
125 done = .FALSE.
126 do while ( .not. done)
127 done = .TRUE.
128 do n = 1, nfiles-1
129 if (filedates(n) > filedates(n+1)) then
130 filedates(size(filedates)) = filedates(n)
131 filedates(n) = filedates(n+1)
132 filedates(n+1) = filedates(size(filedates))
133 filedates(size(filedates)) = '0000-00-00 00:00:00.0000'
134 done = .FALSE.
135 endif
136 enddo
137 enddo
138 endif
139 end subroutine sort_filedates
141 end program g2print
143 !*****************************************************************************!
145 SUBROUTINE r_grib2(junit, gribflnm, hdate, &
146 grib_edition, debug_level, ireaderr)
148 use grib_mod
149 use params
150 use table ! Included to define g2code
151 use gridinfo ! Included to define map%
153 real, allocatable, dimension(:) :: hold_array
154 parameter(msk1=32000,msk2=4000)
155 character(len=1),allocatable,dimension(:) :: cgrib
156 integer :: listsec0(3)
157 integer :: listsec1(13)
158 integer year, month, day, hour, minute, second, fcst
159 character(len=*) :: gribflnm
160 character(len=*) :: hdate
161 character(len=8) :: pabbrev
162 character(len=20) :: labbrev
163 character(len=80) :: tabbrev
164 integer :: lskip, lgrib
165 integer :: junit, itot, icount, iseek
166 integer :: grib_edition
167 integer :: i, j, ireaderr, ith
168 integer :: currlen
169 logical :: unpack, expand
170 type(gribfield) :: gfld
171 real :: level
172 real :: scale_factor
173 integer :: iplvl, lvl2
174 ! For subroutine output
175 integer , parameter :: maxlvl = 150
176 real , dimension(maxlvl) :: plvl
177 integer :: nlvl
178 integer , dimension(maxlvl) :: level_array
179 logical :: verbose=.false.
180 logical :: first = .true.
181 integer :: debug_level
182 character(len=4) :: tmp4
183 character(len=40) :: string
184 character(len=13) :: pstring = ',t50,":",i14)'
185 character(len=15) :: rstring = ',t50,":",f14.5)'
186 character(len=13) :: astring = ',t50,":",a14)'
187 character(len=15) :: estring = ',t50,":",e14.5)'
189 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190 ! SET ARGUMENTS
192 if (debug_level .gt. 50 ) then
193 unpack=.true.
194 else
195 unpack=.false.
196 endif
197 expand=.true.
198 hdate = '0000-00-00_00:00:00'
199 ierr=0
200 itot=0
201 icount=0
202 iseek=0
203 lskip=0
204 lgrib=0
205 currlen=0
206 ith=1
207 scale_factor = 1e6
208 ! do j = 1,10
209 ! write(6,'("j = ",i4," level1 = ",i8," level2 = ",i8)') j, &
210 ! level1(j), level2(j)
211 ! enddo
213 !/* IOS Return Codes from BACIO: */
214 !/* 0 All was well */
215 !/* -1 Tried to open read only _and_ write only */
216 !/* -2 Tried to read and write in the same call */
217 !/* -3 Internal failure in name processing */
218 !/* -4 Failure in opening file */
219 !/* -5 Tried to read on a write-only file */
220 !/* -6 Failed in read to find the 'start' location */
221 !/* -7 Tried to write to a read only file */
222 !/* -8 Failed in write to find the 'start' location */
223 !/* -9 Error in close */
224 !/* -10 Read or wrote fewer data than requested */
226 !if ireaderr =1 we have hit the end of a file.
227 !if ireaderr =2 we have hit the end of all the files.
230 if ( debug_level .gt. 100 ) verbose = .true.
231 if (verbose) write(6,*) 'begin r_grib2, flnm = ',gribflnm
232 ! Open a byte-addressable file.
233 CALL BAOPENR(junit,gribflnm,IOS)
234 first = .true.
235 if (verbose) write(6,*) 'back from baopenr, ios = ',ios
236 if (ios.eq.0) then
237 VERSION: do
239 ! Search opend file for the next GRIB2 messege (record).
240 if (verbose) write(6,*) 'calling skgb'
241 call skgb(junit,iseek,msk1,lskip,lgrib)
243 ! Check for EOF, or problem
244 if (lgrib.eq.0) then
245 exit
246 endif
248 ! Check size, if needed allocate more memory.
249 if (lgrib.gt.currlen) then
250 if (allocated(cgrib)) deallocate(cgrib)
251 allocate(cgrib(lgrib),stat=is)
252 !print *,'G2 allocate(cgrib(lgrib)) status: ',IS
253 currlen=lgrib
254 endif
256 ! Read a given number of bytes from unblocked file.
257 call baread(junit,lskip,lgrib,lengrib,cgrib)
259 if (lgrib.ne.lengrib) then
260 print *,'G2 r_grib2: IO Error.',lgrib,".ne.",lengrib
261 stop 9
262 endif
263 iseek=lskip+lgrib
264 icount=icount+1
266 if (verbose) PRINT *,'G2 GRIB MESSAGE ',icount,' starts at',lskip+1
268 ! Unpack GRIB2 field
269 call gb_info(cgrib,lengrib,listsec0,listsec1, &
270 numfields,numlocal,maxlocal,ierr)
271 if (ierr.ne.0) then
272 write(*,*) ' ERROR querying GRIB2 message = ',ierr
273 stop 10
274 endif
275 itot=itot+numfields
277 grib_edition=listsec0(2)
278 if (grib_edition.ne.2) then
279 exit VERSION
280 endif
282 ! Additional print statments for developer.
283 if (verbose) then
284 print *,'G2 SECTION 0: ',(listsec0(j),j=1,3)
285 print *,'G2 SECTION 1: ',(listsec1(j),j=1,13)
286 print *,'G2 Contains ',numlocal,' Local Sections ', &
287 ' and ',numfields,' data fields.'
288 endif
290 ! ----
291 ! Once per file fill in date, model and projection values.
293 if (first) then
294 first = .false.
296 ! Build the 19-character date string, based on GRIB2 header date
297 ! and time information, including forecast time information:
299 n=1
300 call gf_getfld(cgrib,lengrib,n,unpack,expand,gfld,ierr)
302 if (debug_level .gt. 100 ) then
303 write(6,*) 'gfld%version = ',gfld%version
304 if (gfld%discipline .eq. 0) then
305 string = 'Meteorological products'
306 else if (gfld%discipline .eq. 1) then
307 string = 'Hydrological products'
308 else if (gfld%discipline .eq. 2) then
309 string = 'Land Surface products'
310 else
311 string = 'See code table 0.0'
312 endif
313 write(6,*) 'Discipline = ',gfld%discipline,' ',string
314 write(6,*) 'gfld%idsect(1) = ',gfld%idsect(1)
315 write(6,*) 'gfld%idsect(2) = ',gfld%idsect(2)
316 write(6,*) 'gfld%idsect(3) = ',gfld%idsect(3)
317 write(6,*) 'gfld%idsect(4) = ',gfld%idsect(4)
318 write(6,*) 'gfld%idsect(5) = ',gfld%idsect(5)
319 write(6,*) 'gfld%idsect(6) = ',gfld%idsect(6)
320 write(6,*) 'gfld%idsect(7) = ',gfld%idsect(7)
321 write(6,*) 'gfld%idsect(8) = ',gfld%idsect(8)
322 write(6,*) 'gfld%idsect(9) = ',gfld%idsect(9)
323 write(6,*) 'gfld%idsect(10) = ',gfld%idsect(10)
324 write(6,*) 'gfld%idsect(11) = ',gfld%idsect(11)
325 write(6,*) 'gfld%idsect(12) = ',gfld%idsect(12)
326 write(6,*) 'gfld%idsect(13) = ',gfld%idsect(13)
328 write(6,*) 'gfld%idsectlen = ',gfld%idsectlen
329 write(6,*) 'gfld%locallen = ',gfld%locallen
330 write(6,*) 'gfld%ifldnum = ',gfld%ifldnum
331 write(6,*) 'gfld%ngrdpts = ',gfld%ngrdpts
332 write(6,*) 'gfld%numoct_opt = ',gfld%numoct_opt
333 write(6,*) 'gfld%interp_opt = ',gfld%interp_opt
335 write(6,*) 'gfld%griddef = ',gfld%griddef
336 if (gfld%igdtnum .eq. 0) then
337 string = 'Lat/Lon cylindrical equidistant'
338 else if (gfld%igdtnum .eq. 1) then
339 string = 'Rotated Lat/Lon'
340 else if (gfld%igdtnum .eq. 2) then
341 string = 'Stretched Lat/Lon'
342 else if (gfld%igdtnum .eq. 20) then
343 string = 'Polar Stereographic'
344 else if (gfld%igdtnum .eq. 30) then
345 string = 'Lambert Conformal'
346 else if (gfld%igdtnum .eq. 40) then
347 string = 'Gaussian Lat/Lon'
348 else if (gfld%igdtnum .eq. 50) then
349 string = 'Spherical harmonic coefficients'
350 else
351 string = 'see code table 3.1'
352 endif
353 write(6,*) 'Grid Template number = ',gfld%igdtnum,' ',string
354 write(6,*) 'gfld%igdtlen = ',gfld%igdtlen
355 do i = 1, gfld%igdtlen
356 write(6,*) 'gfld%igdtmpl(',i,') = ',gfld%igdtmpl(i)
357 enddo
359 write(6,*) 'gfld%ipdtnum = ',gfld%ipdtnum
360 write(6,*) 'gfld%ipdtlen = ',gfld%ipdtlen
361 if ( gfld%ipdtnum .eq. 0 ) then
362 do i = 1, gfld%ipdtlen
363 write(6,*) 'gfld%ipdtmpl(',i,') = ',gfld%ipdtmpl(i)
364 enddo
365 endif
366 write(6,*) 'gfld%num_coord = ',gfld%num_coord
367 write(6,*) 'gfld%ndpts = ',gfld%ndpts
368 write(6,*) 'gfld%idrtnum = ',gfld%idrtnum
369 write(6,*) 'gfld%idrtlen = ',gfld%idrtlen
370 write(6,*) 'gfld%expanded = ',gfld%expanded
371 write(6,*) 'gfld%ibmap = ',gfld%ibmap
372 endif
374 year =gfld%idsect(6) !(FOUR-DIGIT) YEAR OF THE DATA
375 month =gfld%idsect(7) ! MONTH OF THE DATA
376 day =gfld%idsect(8) ! DAY OF THE DATA
377 hour =gfld%idsect(9) ! HOUR OF THE DATA
378 minute=gfld%idsect(10) ! MINUTE OF THE DATA
379 second=gfld%idsect(11) ! SECOND OF THE DATA
381 fcst = 0
383 ! Extract forecast time.
384 if ( gfld%ipdtmpl(8) .eq. 1 ) then ! time units are hours
385 fcst = gfld%ipdtmpl(9)
386 else if ( gfld%ipdtmpl(8) .eq. 0 ) then ! minutes
387 fcst = gfld%ipdtmpl(9) / 60.
388 else if ( gfld%ipdtmpl(8) .eq. 2 ) then ! days
389 fcst = gfld%ipdtmpl(9) * 24.
390 else
391 fcst = 999
392 endif
394 ! Compute valid time.
396 if (verbose) then
397 print *, 'ymd',gfld%idsect(6),gfld%idsect(7),gfld%idsect(8)
398 print *, 'hhmm ',gfld%idsect(9),gfld%idsect(10)
399 endif
401 call build_hdate(hdate,year,month,day,hour,minute,second)
402 if (verbose) print *, 'G2 hdate = ',hdate
403 ! call geth_newdate(hdate,hdate,3600*fcst) ! no need for this in print
404 ! print *, 'G2 hdate (fcst?) = ',hdate
406 !--
408 ! Indicator of the source (center) of the data.
409 icenter = gfld%idsect(1)
411 ! Indicator of model (or whatever) which generated the data.
412 iprocess = gfld%ipdtmpl(5)
415 if (icenter.eq.7) then
416 ! Values obtained from http://www.nco.ncep.noaa.gov/pmb/docs/on388/tablea.html
417 ! Note that NCEP recycles process numbers. This may cause labelling issues for
418 ! ancient datasets.
419 if (iprocess.eq.81) then
420 map%source = 'NCEP GFS Analysis'
421 elseif (iprocess.eq.82) then
422 map%source = 'NCEP GFS GDAS/FNL'
423 elseif (iprocess.eq.83) then
424 map%source = 'NCEP HRRR Model'
425 elseif (iprocess.eq.84) then
426 map%source = 'NCEP MESO NAM Model'
427 elseif (iprocess.eq.89) then
428 map%source = 'NCEP NMM '
429 elseif (iprocess.eq.96) then
430 map%source = 'NCEP GFS Model'
431 elseif (iprocess.eq.86 .or. iprocess.eq.100) then
432 map%source = 'NCEP RUC Model' ! 60 km
433 elseif (iprocess.eq.101) then
434 map%source = 'NCEP RUC Model' ! 40 km
435 elseif (iprocess.eq.105) then
436 if (year .gt. 2011) then
437 map%source = 'NCEP RAP Model'
438 else
439 map%source = 'NCEP RUC Model' ! 20 km
440 endif
441 elseif (iprocess.eq.107) then
442 map%source = 'NCEP GEFS'
443 elseif (iprocess.eq.109) then
444 map%source = 'NCEP RTMA'
445 elseif (iprocess.eq.140) then
446 map%source = 'NCEP NARR'
447 elseif (iprocess.eq.44) then
448 map%source = 'NCEP SST Analysis'
449 elseif (iprocess.eq.70) then
450 map%source = 'GFDL Hurricane Model'
451 elseif (iprocess.eq.80) then
452 map%source = 'NCEP GFS Ensemble'
453 elseif (iprocess.eq.107) then ! renumbered as of 23 Feb 2010
454 map%source = 'NCEP GFS Ensemble'
455 elseif (iprocess.eq.111) then
456 map%source = 'NCEP NMMB Model'
457 elseif (iprocess.eq.112) then
458 map%source = 'NCEP WRF-NMM Model'
459 elseif (iprocess.eq.116) then
460 map%source = 'NCEP WRF-ARW Model'
461 elseif (iprocess.eq.129) then
462 map%source = 'NCEP GODAS'
463 elseif (iprocess.eq.197) then
464 map%source = 'NCEP CDAS CFSV2'
465 elseif (iprocess.eq.25) then
466 map%source = 'NCEP SNOW COVER ANALYSIS'
467 else
468 map%source = 'unknown model from NCEP'
469 write (6,*) 'unknown NCEP model, iprocess = ',iprocess
470 end if
471 else if (icenter .eq. 57) then
472 if (iprocess .eq. 87) then
473 map%source = 'AFWA AGRMET'
474 else
475 map%source = 'AFWA'
476 endif
477 else if (icenter .eq. 58) then
478 map%source = 'US Navy FNOC'
479 else if (icenter .eq. 98) then
480 map%source = 'ECMWF'
481 else if (icenter .eq. 34) then
482 map%source = 'JMA'
483 else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
484 map%source = 'UKMO'
485 else
486 map%source = 'unknown model and orig center'
487 end if
488 write (6,*) ' ',map%source
490 if (debug_level .le. 50) then
491 write(6,'(a)') '---------------------------------------------------------------------------------------'
492 write(6,'(a)') ' rec Prod Cat Param Lvl Lvl Lvl Prod Name Time Fcst'
493 write(6,'(a)') ' num Disc num code one two Templ hour'
494 write(6,'(a)') '---------------------------------------------------------------------------------------'
495 endif
498 !--
500 ! Store information about the grid on which the data is.
501 ! This stuff gets stored in the MAP variable, as defined in
502 ! module GRIDINFO.
504 map%startloc = 'SWCORNER'
506 if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
507 map%igrid = 0
508 map%nx = gfld%igdtmpl(8)
509 map%ny = gfld%igdtmpl(9)
510 map%dx = gfld%igdtmpl(17)
511 map%dy = gfld%igdtmpl(18)
512 map%lat1 = gfld%igdtmpl(12)
513 map%lon1 = gfld%igdtmpl(13)
515 if ((gfld%igdtmpl(10) .eq. 0).OR. &
516 (gfld%igdtmpl(10) .eq. 255)) THEN
517 ! Scale lat/lon values to 0-180, default range is 1e6.
518 map%lat1 = map%lat1/scale_factor
519 map%lon1 = map%lon1/scale_factor
520 ! Scale dx/dy values to degrees, default range is 1e6.
521 map%dx = map%dx/scale_factor
522 map%dy = map%dy/scale_factor
523 else
524 ! Basic angle and subdivisions are non-zero (not tested)
525 map%lat1 = map%lat1 * &
526 (gfld%igdtmpl(11)/gfld%igdtmpl(10))
527 map%lon1 = map%lon1 * &
528 (gfld%igdtmpl(11)/gfld%igdtmpl(10))
529 map%dx = map%dx * &
530 (gfld%igdtmpl(11)/gfld%igdtmpl(10))
531 map%dy = map%dy * &
532 (gfld%igdtmpl(11)/gfld%igdtmpl(10))
533 endif
535 elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
536 map%igrid = 3
537 map%nx = gfld%igdtmpl(8)
538 map%ny = gfld%igdtmpl(9)
539 map%lov = gfld%igdtmpl(14)
540 map%truelat1 = gfld%igdtmpl(19)
541 map%truelat2 = gfld%igdtmpl(20)
542 map%dx = gfld%igdtmpl(15)
543 map%dy = gfld%igdtmpl(16)
544 map%lat1 = gfld%igdtmpl(10)
545 map%lon1 = gfld%igdtmpl(11)
547 elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
548 map%igrid = 0
549 map%nx = gfld%igdtmpl(8)
550 map%ny = gfld%igdtmpl(9)
551 map%dx = gfld%igdtmpl(17)
552 map%dy = gfld%igdtmpl(18) ! ?not in Grid Definition Template 3.40 doc
553 map%lat1 = gfld%igdtmpl(12)
554 map%lon1 = gfld%igdtmpl(13)
556 ! Scale dx/dy values to degrees, default range is 1e6.
557 if (map%dx.gt.10000) then
558 map%dx = map%dx/scale_factor
559 endif
560 if (map%dy.gt.10000) then
561 map%dy = (map%dy/scale_factor)*(-1)
562 endif
564 ! Scale lat/lon values to 0-180, default range is 1e6.
565 if (abs(map%lat1).ge.scale_factor) then
566 map%lat1 = map%lat1/scale_factor
567 endif
568 if (map%lon1.ge.scale_factor) then
569 map%lon1 = map%lon1/scale_factor
570 endif
571 print *,'Gaussian Grid: Dx,Dy,lat,lon',map%dx,map%dy, &
572 map%lat1,map%lon1
574 elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
575 map%igrid = 5
576 map%nx = gfld%igdtmpl(8)
577 map%ny = gfld%igdtmpl(9)
578 !map%lov = gfld%igdtmpl(x) ! ?not in Grid Definition Template 3.20 doc
579 map%truelat1 = 60.
580 map%truelat2 = 91.
581 !map%dx = gfld%igdtmpl(x)
582 !map%dy = gfld%igdtmpl(x)
583 map%lat1 = gfld%igdtmpl(10)
584 map%lon1 = gfld%igdtmpl(11)
586 else
587 print*, 'GRIB2 Unknown Projection: ',gfld%igdtnum
588 print*, 'see Code Table 3.1: Grid Definition Template No'
589 endif
591 call gf_free(gfld)
592 endif
594 ! ----
596 ! Continue to unpack GRIB2 field.
597 if (debug_level .gt. 100) write(6,*) 'numfields = ',numfields
598 NUM_FIELDS: do n = 1, numfields
599 ! e.g. U and V would =2, otherwise its usually =1
600 call gf_getfld(cgrib,lengrib,n,unpack,expand,gfld,ierr)
601 if (ierr.ne.0) then
602 write(*,*) ' ERROR extracting field gf_getfld = ',ierr
603 cycle
604 endif
606 ! ------------------------------------
607 ! Additional print information for developer.
608 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1), &
609 gfld%ipdtmpl(2))
610 if (debug_level .gt. 50 ) then
611 print *
612 ! print *,'G2 FIELD ',n
613 if (n==1) then
614 write(*,'(/,"GRIB2 SECTION 0 - INDICATOR SECTION:")')
615 write(*,'(5x,"Discipline"'//pstring) gfld%discipline
616 write(*,'(5x,"GRIB Edition Number"'//pstring) gfld%version
617 write(*,'(5x,"GRIB length"'//pstring) lengrib
618 write(*,'(/,"GRIB2 SECTION 1 - IDENTIFICATION SECTION:")')
619 write(*,'(5x,"Length of Section"'//pstring) gfld%idsectlen
620 write(*,'(5x,"Originating Center ID"'//pstring) &
621 gfld%idsect(1)
622 write(*,'(5x,"Subcenter ID"'//pstring) gfld%idsect(2)
623 write(*,'(5x,"GRIB Master Table Version"'//pstring) &
624 gfld%idsect(3)
625 write(*,'(5x,"GRIB Local Table Version"'//pstring) &
626 gfld%idsect(4)
627 write(*,'(5x,"Significance of Reference Time"'//pstring) &
628 gfld%idsect(5)
629 write(*,'(5x,"Year"'//pstring) gfld%idsect(6)
630 write(*,'(5x,"Month"'//pstring) gfld%idsect(7)
631 write(*,'(5x,"Day"'//pstring) gfld%idsect(8)
632 write(*,'(5x,"Hour"'//pstring) gfld%idsect(9)
633 write(*,'(5x,"Minute"'//pstring) gfld%idsect(10)
634 write(*,'(5x,"Second"'//pstring) gfld%idsect(11)
635 write(*,'(5x,"Production Status of data"'//pstring) &
636 gfld%idsect(12)
637 write(*,'(5x,"Type of processed data"'//pstring) &
638 gfld%idsect(13)
639 ! print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
640 endif
641 write(*,'(/,"GRIB2 SECTION 2 - LOCAL SECTION:")')
642 write(*,'(5x,"Length of Section 2"'//pstring) gfld%locallen
643 if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
644 do j = 1, gfld%locallen
645 write(*,'(5x,"Local value "'//astring) gfld%local(j)
646 enddo
647 ! print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
648 endif
649 write(*,'(/,"GRIB2 SECTION 3 - GRID DEFINITION SECTION:")')
650 ! write(*,'(5x,"Length of Section 3"'//pstring) gfld%unknown
651 write(*,'(5x,"Source of grid definition"'&
652 //pstring) gfld%griddef
653 write(*,'(5x,"Number of grid points"'//pstring) gfld%ngrdpts
654 write(*,'(5x,"Number of octets for addnl points"'//pstring) &
655 gfld%numoct_opt
656 write(*,'(5x,"Interpretation list"'//pstring) &
657 gfld%interp_opt
658 write(*,'(5x,"Grid Definition Template Number"'//pstring) &
659 gfld%igdtnum
660 if (gfld%igdtnum .eq. 0 .or. gfld%igdtnum .eq. 1 .or. &
661 gfld%igdtnum .eq. 2 .or. gfld%igdtnum .eq. 3 ) then
662 if (gfld%igdtnum .eq. 0 ) then
663 write(*,'(5x,"Lat/Lon or Cylindrical Equidistant Grid")')
664 else if (gfld%igdtnum .eq. 1 ) then
665 write(*,'(5x,"Rotated Lat/Lon or Cylind. Equi. Grid")')
666 else if (gfld%igdtnum .eq. 2 ) then
667 write(*,'(5x,"Stretched Lat/Lon or Cylind. Equi. Grid")')
668 else if (gfld%igdtnum .eq. 3 ) then
669 write(*,'(5x,"Stretched and Rotated Lat/Lon Grid")')
670 endif
671 write(*,'(10x,"Shape of the Earth"'//pstring) &
672 gfld%igdtmpl(1)
673 write(*,'(10x,"Scale factor of spher. Earth"'//pstring) &
674 gfld%igdtmpl(2)
675 write(*,'(10x,"Scaled value of spher. Earth"'//pstring) &
676 gfld%igdtmpl(3)
677 write(*,'(10x,"Scale factor of major axis"'//pstring) &
678 gfld%igdtmpl(4)
679 write(*,'(10x,"Scaled value of major axis"'//pstring) &
680 gfld%igdtmpl(5)
681 write(*,'(10x,"Scale factor of minor axis"'//pstring) &
682 gfld%igdtmpl(6)
683 write(*,'(10x,"Scaled value of minor axis"'//pstring) &
684 gfld%igdtmpl(7)
685 write(*,'(10x,"Ni - points along a parallel"'//pstring) &
686 gfld%igdtmpl(8)
687 write(*,'(10x,"Nj - points along a meridian"'//pstring) &
688 gfld%igdtmpl(9)
689 write(*,'(10x,"Basic angle of initial domain"'//pstring)&
690 gfld%igdtmpl(10)
691 write(*,'(10x,"Subdivisions of basic angle"'//pstring) &
692 gfld%igdtmpl(11)
693 write(*,'(10x,"La1"'//pstring) gfld%igdtmpl(12)
694 write(*,'(10x,"Lo1"'//pstring) gfld%igdtmpl(13)
695 write(*,'(10x,"Resolution and Component",t50,":",B14.8)')&
696 gfld%igdtmpl(14)
697 write(*,'(10x,"La2"'//pstring) gfld%igdtmpl(15)
698 write(*,'(10x,"Lo2"'//pstring) gfld%igdtmpl(16)
699 write(*,'(10x,"Di - i-dir increment"'//pstring) &
700 gfld%igdtmpl(17)
701 write(*,'(10x,"Dj - j-dir increment"'//pstring) &
702 gfld%igdtmpl(18)
703 write(*,'(10x,"Scanning mode"'//pstring) &
704 gfld%igdtmpl(19)
705 if (gfld%igdtnum .eq. 1) then
706 write(*,'(10x,"Lat of southern pole of project"'//pstring)&
707 gfld%igdtmpl(20)
708 write(*,'(10x,"Lon of southern pole of project"'//pstring)&
709 gfld%igdtmpl(21)
710 write(*,'(10x,"Angle of rotation of projection"'//pstring)&
711 gfld%igdtmpl(22)
712 else if (gfld%igdtnum .eq. 2) then
713 write(*,'(10x,"Lat of the pole of stretching "'//pstring)&
714 gfld%igdtmpl(20)
715 write(*,'(10x,"Lon of the pole of stretching "'//pstring)&
716 gfld%igdtmpl(21)
717 write(*,'(10x,"Stretching factor"'//pstring) &
718 gfld%igdtmpl(22)
719 else if (gfld%igdtnum .eq. 3) then
720 write(*,'(10x,"Lat of southern pole of project"'//pstring)&
721 gfld%igdtmpl(20)
722 write(*,'(10x,"Lon of southern pole of project"'//pstring)&
723 gfld%igdtmpl(21)
724 write(*,'(10x,"Angle of rotation of projection"'//pstring)&
725 gfld%igdtmpl(22)
726 write(*,'(10x,"Lat of the pole of stretching "'//pstring)&
727 gfld%igdtmpl(23)
728 write(*,'(10x,"Lon of the pole of stretching "'//pstring)&
729 gfld%igdtmpl(24)
730 write(*,'(10x,"Stretching factor"'//pstring) &
731 gfld%igdtmpl(25)
732 endif
733 else if (gfld%igdtnum .eq. 10) then
734 write(*,'(5x,"Mercator Grid")')
735 else if (gfld%igdtnum .eq. 20 .or. gfld%igdtnum .eq. 30) then
736 if (gfld%igdtnum .eq. 20) then
737 write(*,'(5x,"Polar Stereographic Grid")')
738 else if (gfld%igdtnum .eq. 30) then
739 write(*,'(5x,"Lambert Conformal Grid")')
740 endif
741 write(*,'(10x,"Shape of the Earth"'//pstring) &
742 gfld%igdtmpl(1)
743 write(*,'(10x,"Scale factor of spher. Earth"'//pstring) &
744 gfld%igdtmpl(2)
745 write(*,'(10x,"Scaled value of spher. Earth"'//pstring) &
746 gfld%igdtmpl(3)
747 write(*,'(10x,"Scale factor of major axis"'//pstring) &
748 gfld%igdtmpl(4)
749 write(*,'(10x,"Scaled value of major axis"'//pstring) &
750 gfld%igdtmpl(5)
751 write(*,'(10x,"Scale factor of minor axis"'//pstring) &
752 gfld%igdtmpl(6)
753 write(*,'(10x,"Scaled value of minor axis"'//pstring) &
754 gfld%igdtmpl(7)
755 write(*,'(10x,"Nx"'//pstring) gfld%igdtmpl(8)
756 write(*,'(10x,"Ny"'//pstring) gfld%igdtmpl(9)
757 write(*,'(10x,"La1"'//pstring) gfld%igdtmpl(10)
758 write(*,'(10x,"Lo1"'//pstring) gfld%igdtmpl(11)
759 write(*,'(10x,"Resolution and Component",t50,":",B14.8)')&
760 gfld%igdtmpl(12)
761 write(*,'(10x,"LaD"'//pstring) gfld%igdtmpl(13)
762 write(*,'(10x,"LoV"'//pstring) gfld%igdtmpl(14)
763 write(*,'(10x,"Dx"'//pstring) gfld%igdtmpl(15)
764 write(*,'(10x,"Dy"'//pstring) gfld%igdtmpl(16)
765 write(*,'(10x,"Projection Center Flag"'//pstring) &
766 gfld%igdtmpl(17)
767 write(*,'(10x,"Scanning mode"'//pstring) &
768 gfld%igdtmpl(18)
769 if (gfld%igdtnum .eq. 30) then
770 write(*,'(10x,"Latin 1 "'//pstring) &
771 gfld%igdtmpl(19)
772 write(*,'(10x,"Latin 2 "'//pstring) &
773 gfld%igdtmpl(20)
774 write(*,'(10x,"Lat of southern pole of project"'//pstring)&
775 gfld%igdtmpl(21)
776 write(*,'(10x,"Lon of southern pole of project"'//pstring)&
777 gfld%igdtmpl(22)
778 endif
779 else if (gfld%igdtnum .eq. 40 .or. gfld%igdtnum .eq. 41) then
780 if (gfld%igdtnum .eq. 40) then
781 write(*,'(5x,"Gaussian Lat/Lon Grid")')
782 else if (gfld%igdtnum .eq. 41) then
783 write(*,'(5x,"Rotated Gaussian Lat/Lon Grid")')
784 else if (gfld%igdtnum .eq. 42) then
785 write(*,'(5x,"Stretched Gaussian Lat/Lon Grid")')
786 else if (gfld%igdtnum .eq. 43) then
787 write(*,'(5x,"Stretched and Rotated Gaussian Lat/Lon ")')
788 endif
789 else
790 do j = 1, gfld%igdtlen
791 write(*,'(5x,"Grid Definition Template entry "'//pstring) &
792 gfld%igdtmpl(j)
793 enddo
794 endif
795 ! print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts, &
796 ! gfld%numoct_opt,gfld%interp_opt, &
797 ! gfld%igdtnum
798 ! print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ', &
799 ! (gfld%igdtmpl(j),j=1,gfld%igdtlen)
800 if ( gfld%num_opt .eq. 0 ) then
801 ! print *,'G2 NO Section 3 List Defining No. of Data Points.'
802 else
803 print *,'G2 Section 3 Optional List: ', &
804 (gfld%list_opt(j),j=1,gfld%num_opt)
805 endif
806 write(*,'(/,"GRIB2 SECTION 4 - PRODUCT DEFINITION SECTION:")')
807 ! write(*,'(5x,"Length of Section 4"'//pstring) gfld%unknown
808 write(*,'(5x,"Product Definition Template Number"'//pstring)&
809 gfld%ipdtnum
810 do j = 1, gfld%ipdtlen
811 write(tmp4,'(i4)') j
812 string = '(5x,"Template Entry '//tmp4 // '"'
813 write(*,string//pstring) gfld%ipdtmpl(j)
814 enddo
815 ! print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ', &
816 ! (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
818 !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
819 !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
820 write(*,'(5x,"Product Abbreviated Name",t50,":",a14)')&
821 pabbrev
822 ! print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
824 if ( gfld%num_coord .eq. 0 ) then
825 ! print *,'G2 NO Optional Vertical Coordinate List.'
826 else
827 print *,'G2 Section 4 Optional Coordinates: ', &
828 (gfld%coord_list(j),j=1,gfld%num_coord)
829 endif
830 ! if ( gfld%ibmap .ne. 255 ) then
831 ! print *,'G2 Num. of Data Points = ',gfld%ndpts, &
832 ! ' with BIT-MAP ',gfld%ibmap
833 ! else
834 ! print *,'G2 Num. of Data Points = ',gfld%ndpts, &
835 ! ' NO BIT-MAP '
836 ! endif
837 write(*,'(/,"GRIB2 SECTION 5 - DATA REPRESENTATION SECTION:")')
838 write(*,'(5x,"Data Representation Template Number"'//pstring)&
839 gfld%idrtnum
840 do j = 1, gfld%idrtlen
841 write(tmp4,'(i4)') j
842 string = '(5x,"Template Entry '//tmp4 // '"'
843 write(*,string//pstring) gfld%idrtmpl(j)
844 enddo
845 ! print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ', &
846 ! (gfld%idrtmpl(j),j=1,gfld%idrtlen)
848 ! if ( gfld%ipdtnum .eq. 0 ) then
849 ! if (gfld%ipdtmpl(1) .eq. 0 ) then
850 ! write(6,*) 'Temperature'
851 ! else if (gfld%ipdtmpl(1) .eq. 1 ) then
852 ! write(6,*) 'Moisture'
853 ! else if (gfld%ipdtmpl(1) .eq. 2 ) then
854 ! write(6,*) 'Momentum'
855 ! else if (gfld%ipdtmpl(1) .eq. 3 ) then
856 ! write(6,*) 'Mass'
857 ! endif
858 ! endif
860 write(*,'(/,"GRIB2 SECTION 6 - BIT-MAP SECTION:")')
861 write(*,'(5x,"Bit-map indicator"'//pstring) &
862 gfld%ibmap
864 fldmax=gfld%fld(1)
865 fldmin=gfld%fld(1)
866 sum=gfld%fld(1)
867 do j=2,gfld%ndpts
868 if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
869 if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
870 sum=sum+gfld%fld(j)
871 enddo ! gfld%ndpts
873 write(*,'(/,"GRIB2 SECTION 7 - DATA SECTION:")')
875 if ( fldmax .lt. -1.e10 ) then
876 write(*,'(5x,"Minimum Data Value "'//estring)&
877 fldmin
878 else
879 write(*,'(5x,"Minimum Data Value "'//rstring)&
880 fldmin
881 endif
882 if ( fldmax .gt. 1.e10 ) then
883 write(*,'(5x,"Maximum Data Value "'//estring)&
884 fldmax
885 else
886 write(*,'(5x,"Maximum Data Value "'//rstring)&
887 fldmax
888 endif
889 ! print *,'G2 Data Values:'
890 ! write(*,fmt='("G2 MIN=",f21.8," AVE=",f21.8, &
891 ! " MAX=",f21.8)') fldmin,sum/gfld%ndpts,fldmax
892 if (debug_level .gt. 100 ) then
893 print*, 'gfld%fld = ',gfld%fld
894 ! do j=1,gfld%ndpts
895 ! write(*,*) j, gfld%fld(j)
896 ! enddo
897 endif
898 endif ! Additional Print information
899 ! ------------------------------------
900 if ( debug_level .le. 50) then
901 if(gfld%ipdtmpl(10).eq.100) then ! pressure level
902 level=gfld%ipdtmpl(12) * &
903 (10. ** (-1. * gfld%ipdtmpl(11)))
904 else if(gfld%ipdtmpl(10).eq.101 .or.& ! sea level, sfc, or trop
905 gfld%ipdtmpl(10).eq.1 .or. gfld%ipdtmpl(10).eq.7) then
906 level = 0
907 else if(gfld%ipdtmpl(10).eq.106) then ! below ground sfc is in cm in Vtable
908 level= 100. * gfld%ipdtmpl(12)*(10.**(-1.*gfld%ipdtmpl(11)))
909 else
910 level=gfld%ipdtmpl(12) * 10.** (-1.*gfld%ipdtmpl(11))
911 endif
912 if (gfld%ipdtmpl(13) .eq. 255) then
913 lvl2 = 0
914 else if(gfld%ipdtmpl(10).eq.106) then ! below ground sfc is in cm in Vtable
915 lvl2 = 100. * gfld%ipdtmpl(15) * 10.** (-1.*gfld%ipdtmpl(14))
916 else
917 lvl2 = gfld%ipdtmpl(15) * 10.** (-1.*gfld%ipdtmpl(14))
918 endif
919 ! Account for multiple forecast hours in one file
920 if (gfld%ipdtnum.eq.0 .or. gfld%ipdtnum.eq.1 .or. gfld%ipdtnum.eq. 8) then
921 ! Product Definition Template 4.0, 4.1, 4.8
922 ! Extract forecast time.
923 if ( gfld%ipdtmpl(8) .eq. 1 ) then ! time units are hours
924 fcst = gfld%ipdtmpl(9)
925 else if ( gfld%ipdtmpl(8) .eq. 0 ) then ! minutes
926 fcst = gfld%ipdtmpl(9) / 60.
927 else if ( gfld%ipdtmpl(8) .eq. 2 ) then ! days
928 fcst = gfld%ipdtmpl(9) * 24.
929 else
930 fcst = 999
931 endif
932 endif
934 ! Non-standard Product Definition Templates need to be reported
935 string = ' '
936 if ( gfld%ipdtnum .eq. 8 ) then
937 string = ' PDT4.8'
938 else if ( gfld%ipdtnum .eq. 1 ) then
939 string = ' PDT4.1'
940 endif
941 write(6,987) itot,gfld%discipline,gfld%ipdtmpl(1), &
942 gfld%ipdtmpl(2),gfld%ipdtmpl(10),int(level),&
943 lvl2,gfld%ipdtnum,pabbrev,hdate,fcst,string
944 987 format(2i4,i5,i4,i8,i8,i8,i8,3x,a10,a20,i5.2,a10)
946 endif
948 ! Deallocate arrays decoding GRIB2 record.
949 call gf_free(gfld)
951 enddo NUM_FIELDS
953 enddo VERSION ! skgb
955 if (debug_level .gt. 50) &
956 print *, 'G2 total number of fields found = ',itot
958 CALL BACLOSE(junit,IOS)
960 ireaderr=1
961 else
962 print *,'open status failed because',ios
963 hdate = '9999-99-99_99:99:99'
964 ireaderr=2
965 endif ! ireaderr check
967 END subroutine r_grib2
969 !*****************************************************************************!
970 ! Subroutine edition_num !
971 ! !
972 ! Purpose: !
973 ! Read one record from the input GRIB2 file. Based on the information in !
974 ! the GRIB2 header and the user-defined Vtable, decide whether the field in!
975 ! the GRIB2 record is one to process or to skip. If the field is one we !
976 ! want to keep, extract the data from the GRIB2 record, and pass the data !
977 ! back to the calling routine. !
978 ! !
979 ! Argument list: !
980 ! Input: !
981 ! JUNIT : "Unit Number" to open and read from. Not really a Fortran !
982 ! unit number, since we do not do Fortran I/O for the GRIB2 !
983 ! files. Nor is it a UNIX File Descriptor returned from a C !
984 ! OPEN statement. It is really just an array index to the !
985 ! array (IUARR) where the UNIX File Descriptor values are !
986 ! stored. !
987 ! GRIB2FILE: File name to open, if it is not already open. !
988 ! !
989 ! Output: !
990 ! GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2 !
991 ! IERR : Error flag: 0 - no error on read from GRIB2 file. !
992 ! 1 - Hit the end of the GRIB2 file. !
993 ! 2 - The file GRIBFLNM we tried to open does !
994 ! not exist. !
995 ! Author: Paula McCaslin !
996 ! NOAA/FSL !
997 ! Sept 2004 !
998 !*****************************************************************************!
1000 SUBROUTINE edition_num(junit, gribflnm, grib_edition, ireaderr)
1002 use grib_mod
1003 use params
1005 parameter(msk1=32000,msk2=4000)
1006 character(len=1),allocatable,dimension(:) :: cgrib
1007 integer :: listsec0(3)
1008 integer :: listsec1(13)
1009 character(len=*) :: gribflnm
1010 integer :: lskip, lgrib
1011 integer :: junit
1012 integer :: grib_edition
1013 integer :: i, j, ireaderr
1014 integer :: currlen
1016 character(len=4) :: ctemp
1017 character(len=4),parameter :: grib='GRIB',c7777='7777'
1019 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1020 ! SET ARGUMENTS
1022 itot=0
1023 icount=0
1024 iseek=0
1025 lskip=0
1026 lgrib=0
1027 currlen=0
1029 !/* IOS Return Codes from BACIO: */
1030 !/* 0 All was well */
1031 !/* -1 Tried to open read only _and_ write only */
1032 !/* -2 Tried to read and write in the same call */
1033 !/* -3 Internal failure in name processing */
1034 !/* -4 Failure in opening file */
1035 !/* -5 Tried to read on a write-only file */
1036 !/* -6 Failed in read to find the 'start' location */
1037 !/* -7 Tried to write to a read only file */
1038 !/* -8 Failed in write to find the 'start' location */
1039 !/* -9 Error in close */
1040 !/* -10 Read or wrote fewer data than requested */
1042 !if ireaderr =1 we have hit the end of a file.
1043 !if ireaderr =2 we have hit the end of all the files.
1044 !if ireaderr =3 beginning characters 'GRIB' not found
1046 ! write(6,*) 'junit = ',junit,' gribflnm = ',gribflnm
1048 ! Open a byte-addressable file.
1049 CALL BAOPENR(junit,gribflnm,IOS) ! from w3lib
1050 ! write(6,*) 'ios = ',ios
1051 if (ios.eq.0) then
1053 ! Search opend file for the next GRIB2 messege (record).
1054 call skgb(junit,iseek,msk1,lskip,lgrib)
1056 ! Check for EOF, or problem
1057 if (lgrib.eq.0) then
1058 write(*,'("\n\tThere is a problem with the input file.")')
1059 write(*,'("\tPerhaps it is not a Grib2 file?\n")')
1060 STOP "Grib2 file or date problem, stopping in edition_num."
1061 endif
1063 ! Check size, if needed allocate more memory.
1064 if (lgrib.gt.currlen) then
1065 if (allocated(cgrib)) deallocate(cgrib)
1066 allocate(cgrib(lgrib),stat=is)
1067 currlen=lgrib
1068 endif
1070 ! Read a given number of bytes from unblocked file.
1071 call baread(junit,lskip,lgrib,lengrib,cgrib)
1073 ! Check for beginning of GRIB message in the first 100 bytes
1074 istart=0
1075 do j=1,100
1076 ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
1077 if (ctemp.eq.grib ) then
1078 istart=j
1079 exit
1080 endif
1081 enddo
1082 if (istart.eq.0) then
1083 ireaderr=3
1084 print*, "The beginning 4 characters >GRIB< were not found."
1085 endif
1087 ! Unpack Section 0 - Indicator Section to extract GRIB edition field
1088 iofst=8*(istart+5)
1089 call gbyte(cgrib,discipline,iofst,8) ! Discipline
1090 iofst=iofst+8
1091 call gbyte(cgrib,grib_edition,iofst,8) ! GRIB edition number
1093 print *, 'ungrib - grib edition num', grib_edition
1094 CALL BACLOSE(junit,IOS)
1095 ireaderr=1
1096 else if (ios .eq. -4) then
1097 print *,'edition_num: unable to open ',gribflnm
1098 stop 'edition_num'
1099 else
1100 print *,'edition_num: open status failed because',ios,gribflnm
1101 ireaderr=2
1102 endif ! ireaderr check
1104 END subroutine edition_num
1105 subroutine parse_args(err, a1, h1, i1, l1, a2, h2, i2, l2, a3, h3, i3, l3, &
1106 hlast)
1107 integer :: err
1108 character(len=*) , optional :: a1, a2, a3
1109 character(len=*), optional :: h1, h2, h3
1110 integer , optional :: i1, i2, i3
1111 logical, optional :: l1, l2, l3
1112 character(len=*), optional :: hlast
1114 character(len=100) :: hold
1115 integer :: ioff = 0
1117 if (present(hlast)) then
1118 ioff = -1
1119 endif
1121 err = 0
1123 narg = iargc()
1124 numarg = narg + ioff
1126 i = 1
1127 LOOP : do while ( i <= numarg)
1129 ierr = 1
1130 if (present(i1)) then
1131 call checkiarg(i, a1, i1, ierr)
1132 elseif (present(h1)) then
1133 call checkharg(i, a1, h1, ierr)
1134 elseif (present(l1)) then
1135 call checklarg(i, a1, l1, ierr)
1136 endif
1137 if (ierr.eq.0) cycle LOOP
1139 if (present(i2)) then
1140 call checkiarg(i, a2, i2, ierr)
1141 elseif (present(h2)) then
1142 call checkharg(i, a2, h2, ierr)
1143 elseif (present(l2)) then
1144 call checklarg(i, a2, l2, ierr)
1145 endif
1146 if (ierr.eq.0) cycle LOOP
1148 if (present(i3)) then
1149 call checkiarg(i, a3, i3, ierr)
1150 elseif (present(h3)) then
1151 call checkharg(i, a3, h3, ierr)
1152 elseif (present(l3)) then
1153 call checklarg(i, a3, l3, ierr)
1154 endif
1155 if (ierr.eq.0) cycle LOOP
1157 err = 1
1158 call getarg(1, hold)
1159 write(*, '("arg = ", A)') trim(hold)
1161 exit LOOP
1163 enddo LOOP
1165 if (present(hlast)) then
1166 if (narg.eq.0) then
1167 err = 1
1168 else
1169 call getarg(narg, hlast)
1170 endif
1171 endif
1173 contains
1174 subroutine checkiarg(c, a, i, ierr)
1175 integer :: c
1176 character(len=*) :: a
1177 integer :: i
1179 character(len=100) :: hold
1180 ierr = 1
1182 call getarg(c, hold)
1184 if ('-'//a.eq.trim(hold)) then
1185 c = c + 1
1186 call getarg(c, hold)
1187 read(hold, *) i
1188 c = c + 1
1189 ierr = 0
1190 elseif ('-'//a .eq. hold(1:len_trim(a)+1)) then
1191 hold = hold(len_trim(a)+2: len(hold))
1192 read(hold, *) i
1193 c = c + 1
1194 ierr = 0
1195 endif
1197 end subroutine checkiarg
1198 subroutine checkharg(c, a, h, ierr)
1199 integer :: c
1200 character(len=*) :: a
1201 character(len=*) :: h
1203 character(len=100) :: hold
1204 ierr = 1
1206 call getarg(c, hold)
1208 if ('-'//a.eq.trim(hold)) then
1209 c = c + 1
1210 call getarg(c, hold)
1211 h = trim(hold)
1212 c = c + 1
1213 ierr = 0
1214 elseif ('-'//a .eq. hold(1:len_trim(a)+1)) then
1215 hold = hold(len_trim(a)+2: len(hold))
1216 h = trim(hold)
1217 c = c + 1
1218 ierr = 0
1219 endif
1221 end subroutine checkharg
1223 subroutine checklarg(c, a, l, ierr)
1224 integer :: c
1225 character(len=*) :: a
1226 logical :: l
1228 character(len=100) :: hold
1229 ierr = 1
1231 call getarg(c, hold)
1232 if ('-'//a.eq.trim(hold)) then
1233 l = .TRUE.
1234 c = c + 1
1235 ierr = 0
1236 endif
1238 end subroutine checklarg
1240 end subroutine parse_args