Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

g2print.F90 @ 87

Last change on this file since 87 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 40.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format