Context Navigation

class_gribfile_mod.f90 @ 027e844

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie

Last change on this file since 027e844 was 61e07ba, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago
Adding files from CTBTO project wo_17
Property mode set to `100644`
File size: 15.2 KB

Line
1	MODULE class_gribfile
2
3	!*****************************************************************
4	! *
5	! This is a class-like structure for providing access to the *
6	! metadata and data within a GRIB met file *
7	! *
8	! *
9	!*****************************************************************
10
11	IMPLICIT NONE
12	PRIVATE ! The default is that everything is private, unless
13	! specified otherwise
14
15	! Note that all public interfaces and variables should have a
16	! GRIBFILE_ prefix
17
18	PUBLIC :: gribfile_object, &
19	gribfile_object_create, &
20	gribfile_printobj, &
21	gribfile_centre, &
22	gribfile_num_xlon, &
23	gribfile_num_ylat, &
24	gribfile_num_zlevel
25
26
27	PUBLIC :: GRIBFILE_TYPE_ECMWF_GRIB1, &
28	GRIBFILE_TYPE_ECMWF_GRIB2, &
29	GRIBFILE_TYPE_ECMWF_GRIB1_2, &
30	GRIBFILE_TYPE_NCEP_GRIB1, &
31	GRIBFILE_TYPE_NCEP_GRIB2, &
32	GRIBFILE_TYPE_UNKNOWN, &
33	GRIBFILE_CENTRE_NCEP, &
34	GRIBFILE_CENTRE_ECMWF, &
35	GRIBFILE_CENTRE_UNKNOWN
36
37	! These are codes for designating the type of GRIB file
38	! being looked at
39	INTEGER, PARAMETER :: GRIBFILE_TYPE_ECMWF_GRIB1 = 1, &
40	GRIBFILE_TYPE_ECMWF_GRIB2 = 2, &
41	GRIBFILE_TYPE_ECMWF_GRIB1_2 = 3, &
42	GRIBFILE_TYPE_NCEP_GRIB1 = 4, &
43	GRIBFILE_TYPE_NCEP_GRIB2 = 5, &
44	GRIBFILE_TYPE_UNKNOWN = -9999, &
45	GRIBFILE_CENTRE_NCEP = 1, &
46	GRIBFILE_CENTRE_ECMWF = 2, &
47	GRIBFILE_CENTRE_UNKNOWN = -9999
48
49	! These are the official centre codes for NCEP and ECMWF in grib files.
50	INTEGER, PARAMETER :: CENTRE_NCEP = 7, CENTRE_ECMWF = 98
51
52	TYPE gribfile_object
53	PRIVATE ! Make everything in here private so it's not directly manipulated outside
54	LOGICAL :: is_instantiated = .FALSE.
55	CHARACTER(LEN=256) :: file_path = ''
56	INTEGER :: grib_edition = 0 ! Not sure we want this, since it can vary on hybrid files
57	INTEGER :: grib_centre = GRIBFILE_CENTRE_UNKNOWN
58	INTEGER :: gribfile_type = GRIBFILE_TYPE_UNKNOWN
59	INTEGER :: num_xlon = -9999
60	INTEGER :: num_ylat = -9999
61	INTEGER :: num_zlevel = -9999
62	END TYPE gribfile_object
63
64
65
66	CONTAINS
67
68	SUBROUTINE gribfile_testhello()
69	PRINT *, 'Hello gribfile'
70	END SUBROUTINE gribfile_testhello
71
72
73
74	TYPE(gribfile_object) FUNCTION gribfile_object_create(filepath)
75
76	! This is the "constructor" for the pseudo gribfile object. Given the path to a gribfile,
77	! fills in attributes that can be accessed through methods in this
78	! module
79
80	USE grib_api
81
82	IMPLICIT NONE
83
84	CHARACTER(LEN=*), INTENT(IN) :: filepath ! full path to GRIB file
85
86	TYPE(gribfile_object) :: returned_object
87
88	INTEGER :: ifile, iret, igrib, grib_centre, gribfile_type
89	INTEGER :: xlon, ylat, zlev
90
91
92	CALL get_centre_and_type(filepath, grib_centre, gribfile_type)
93	returned_object%grib_centre = grib_centre
94	returned_object%gribfile_type = gribfile_type
95
96	! Get dimensions of 3d u field
97	CALL get_3d_u_dims(filepath, gribfile_type, xlon, ylat, zlev)
98	returned_object%num_xlon = xlon
99	returned_object%num_ylat = ylat
100	returned_object%num_zlevel = zlev
101
102
103	returned_object%is_instantiated = .TRUE.
104	returned_object%file_path = TRIM(filepath)
105
106	gribfile_object_create = returned_object
107
108	END FUNCTION gribfile_object_create
109
110
111	SUBROUTINE gribfile_printobj(gribfile_obj)
112
113	! Pretty prints the attributes of the gribfile pseudo-object
114	TYPE(gribfile_object), INTENT(IN) :: gribfile_obj
115
116	PRINT *, 'is_instantiated: ', gribfile_obj%is_instantiated
117	PRINT *, 'filepath: ', TRIM(gribfile_obj%file_path)
118	PRINT *, 'grib_centre: ', gribfile_obj%grib_centre
119	PRINT *, 'gribfile_type: ', gribfile_obj%gribfile_type
120	PRINT *, 'num_xlon: ', gribfile_obj%num_xlon
121	PRINT *, 'num_ylat: ', gribfile_obj%num_ylat
122	PRINT *, 'num_zlevel: ', gribfile_obj%num_zlevel
123
124	END SUBROUTINE gribfile_printobj
125
126	INTEGER FUNCTION gribfile_centre(filepath)
127
128	! Returns an integer constant denoting the grib centre (currently either ECMWF, NCEP or UNKNOWN)
129	! for the specified filepath. Returns one of the GRIBFILE_CENTRE_ constants defined at top of this
130	! module.
131
132
133	USE grib_api
134
135	IMPLICIT NONE
136
137	CHARACTER(LEN=*), INTENT(IN) :: filepath ! full path to GRIB file
138
139
140
141	INTEGER :: ifile, iret, igrib, grib_centre
142
143	CALL grib_open_file(ifile, filepath, 'r', iret)
144	IF (iret == 0) THEN
145	! Use first record to detect centre, which is assumed constant
146	! amongst all messages
147	CALL grib_new_from_file(ifile, igrib, iret)
148	CALL grib_get(igrib, 'centre', grib_centre)
149	CALL grib_close_file(ifile)
150	ELSE
151	PRINT *, "WARNING: problem opening GRIB file: ", filepath
152	grib_centre = -999
153	END IF
154
155
156
157
158
159	IF (grib_centre == CENTRE_NCEP) THEN
160	gribfile_centre = GRIBFILE_CENTRE_NCEP
161	ELSE IF (grib_centre == CENTRE_ECMWF) THEN
162	gribfile_centre = GRIBFILE_CENTRE_ECMWF
163	ELSE
164	gribfile_centre = GRIBFILE_CENTRE_UNKNOWN
165	END IF
166
167	END FUNCTION gribfile_centre
168
169
170	! This is currently a PRIVATE subroutine
171	SUBROUTINE get_centre_and_type(filepath, grib_centre, gribfile_type)
172	! Given specified grib file, passes back centre and gribfile
173	! type to the calling program. Numeric codes are defined as integer parameters
174	! in this module
175
176	! To get this information, we have to iterate through the entire file in order to
177	! determine if it is hybrid or not
178	!
179
180	USE grib_api
181	IMPLICIT NONE
182	CHARACTER(LEN=*), INTENT(IN) :: filepath ! full path to GRIB file
183	INTEGER, INTENT(OUT) :: grib_centre, gribfile_type
184
185	INTEGER :: fileptr, iret, igrib, centre, grib_edition
186	LOGICAL :: end_of_file
187	LOGICAL :: grib1_detected, grib2_detected
188
189
190	grib1_detected = .FALSE.
191	grib2_detected = .FALSE.
192
193	CALL grib_open_file(fileptr, filepath, 'r', iret)
194	IF (iret /= 0) THEN
195	PRINT *, 'class_gributils:get_centre_and_type()...'
196	PRINT *, ' unable to open filepath: ', filepath
197	STOP
198	END IF
199
200
201	! Use first record to detect centre and and grib version of first messages. We will
202	! then assume that all following messages have same centre, but not necessarily same
203	! GRIB version
204	CALL grib_new_from_file(fileptr, igrib, iret)
205	CALL grib_get(igrib, 'centre', grib_centre)
206	CALL grib_get(igrib, 'edition', grib_edition)
207
208	IF (grib_edition == 1) grib1_detected = .TRUE.
209	IF (grib_edition == 2) grib2_detected = .TRUE.
210
211	! Now, iterate through the rest of records to determine if this is a mixed edition file
212	end_of_file = .FALSE.
213	DO WHILE (.NOT. end_of_file)
214	CALL grib_new_from_file(fileptr, igrib, iret)
215	IF (iret .eq. GRIB_END_OF_FILE) THEN
216	end_of_file = .TRUE.
217	ELSE
218
219	! Get edition from file
220	CALL grib_get(igrib, 'edition', grib_edition)
221	IF (grib_edition .eq. 1) grib1_detected = .TRUE.
222	IF (grib_edition .eq. 2) grib2_detected = .TRUE.
223	END IF
224	END DO
225
226	CALL grib_close_file(fileptr)
227
228	! Determine the gribfile type depending on centre and edition(s)
229	IF (grib_centre == CENTRE_ECMWF) THEN
230	IF (grib1_detected .AND. grib2_detected) THEN
231	gribfile_type = GRIBFILE_TYPE_ECMWF_GRIB1_2
232	ELSE IF (grib1_detected .AND. .NOT. grib2_detected) THEN
233	gribfile_type = GRIBFILE_TYPE_ECMWF_GRIB1
234	ELSE IF (.NOT. grib1_detected .AND. grib2_detected) THEN
235	gribfile_type = GRIBFILE_TYPE_ECMWF_GRIB2
236	ELSE
237	gribfile_type = GRIBFILE_TYPE_UNKNOWN
238	ENDIF
239	ELSE IF (grib_centre == CENTRE_NCEP) THEN
240	IF (grib1_detected .AND. .NOT. grib2_detected) THEN
241	gribfile_type = GRIBFILE_TYPE_NCEP_GRIB1
242	ELSE IF (.NOT. grib1_detected .AND. grib2_detected) THEN
243	gribfile_type = GRIBFILE_TYPE_NCEP_GRIB2
244	ELSE
245	gribfile_type = GRIBFILE_TYPE_UNKNOWN
246	ENDIF
247	ELSE
248	gribfile_type = GRIBFILE_TYPE_UNKNOWN
249	ENDIF
250
251	END SUBROUTINE get_centre_and_type
252
253
254	SUBROUTINE get_3d_u_dims(filepath, gribfile_type, xlon, ylat, zlev)
255
256	! Looks at the 3d u fields in the GRIBFILE to get x and y dims, as well as number of levels
257	USE grib_api
258
259	IMPLICIT NONE
260
261	CHARACTER(LEN=*), INTENT(IN) :: filepath ! full path to GRIB file
262	INTEGER, INTENT(IN) :: gribfile_type
263	INTEGER, INTENT(OUT) :: xlon, ylat, zlev
264
265	INTEGER :: ifile, iret, igrib, grib_centre
266	LOGICAL :: end_of_file
267
268	! These will be assigned according to type of grib file, then used to filter
269	! for the 3d messages
270	! Name of the key being sought
271	CHARACTER(LEN=64) :: keyname_leveltype, keyname_shortname, keyname_level, &
272	keyname_xlon, keyname_ylat
273
274	! The key value being filtered for
275	CHARACTER(LEN=64) :: keyvalue_leveltype, keyvalue_shortname
276
277	! Actual values read in from the grib file
278	CHARACTER(LEN=64) :: value_leveltype, value_shortname
279	INTEGER :: value_level
280
281	INTEGER :: num_levels
282
283	! Get the field names to read, based on the type of grib file
284
285	!!! Note that ALL of these have the same key names, except that
286	!!! leveltype is 'hybrid' in ECMWF and 'isobaricInhPa' in NCEP.
287	!!! This could probably be consolidated, but because these are
288	!!! files that go through some preprocessing, and aren't
289	!!! necessarily standard (at least for ECMWF), I'm going to be
290	!!! safe and leave it as is for now, so it would be easier to
291	!!! modify for one type, if necessary.
292	IF (gribfile_type == GRIBFILE_TYPE_ECMWF_GRIB1) THEN
293	keyname_leveltype = 'typeOfLevel'
294	keyname_shortname = 'shortName'
295	keyname_level = 'level'
296	keyvalue_leveltype = 'hybrid'
297	keyvalue_shortname = 'u'
298	keyname_xlon = 'Ni'
299	keyname_ylat = 'Nj'
300	ELSE IF (gribfile_type == GRIBFILE_TYPE_ECMWF_GRIB1_2) THEN
301	keyname_leveltype = 'typeOfLevel'
302	keyname_shortname = 'shortName'
303	keyname_level = 'level'
304	keyvalue_leveltype = 'hybrid'
305	keyvalue_shortname = 'u'
306	keyname_xlon = 'Ni'
307	keyname_ylat = 'Nj'
308	ELSE IF (gribfile_type == GRIBFILE_TYPE_ECMWF_GRIB2) THEN
309	keyname_leveltype = 'typeOfLevel'
310	keyname_shortname = 'shortName'
311	keyname_level = 'level'
312	keyvalue_leveltype = 'hybrid'
313	keyvalue_shortname = 'u'
314	keyname_xlon = 'Ni'
315	keyname_ylat = 'Nj'
316	ELSE IF (gribfile_type == GRIBFILE_TYPE_NCEP_GRIB1) THEN
317	keyname_leveltype = 'typeOfLevel'
318	keyname_shortname = 'shortName'
319	keyname_level = 'level'
320	keyvalue_leveltype = 'isobaricInhPa'
321	keyvalue_shortname = 'u'
322	keyname_xlon = 'Ni'
323	keyname_ylat = 'Nj'
324	ELSE IF (gribfile_type == GRIBFILE_TYPE_NCEP_GRIB2) THEN
325	keyname_leveltype = 'typeOfLevel'
326	keyname_shortname = 'shortName'
327	keyname_level = 'level'
328	keyvalue_leveltype = 'isobaricInhPa'
329	keyvalue_shortname = 'u'
330	keyname_xlon = 'Ni'
331	keyname_ylat = 'Nj'
332	ELSE
333	PRINT *, 'class_gribfile:get_3d_u_dims(): Unsupported gribfile type: ', gribfile_type
334	STOP
335	ENDIF
336
337	CALL grib_open_file(ifile, filepath, 'r', iret)
338	IF (iret == 0) THEN
339
340	! Iterate through all messages to count 3d u messages (levels) and get x,y dimensions
341	end_of_file = .FALSE.
342	num_levels = 0
343	DO WHILE (.NOT. end_of_file)
344	CALL grib_new_from_file(ifile, igrib, iret)
345	IF (iret .eq. GRIB_END_OF_FILE) THEN
346	end_of_file = .TRUE.
347	ELSE
348
349	! Get relevant keys and filter for the 3d U wind
350	CALL grib_get(igrib, keyname_shortname, value_shortname)
351	CALL grib_get(igrib, keyname_leveltype, value_leveltype)
352	CALL grib_get(igrib, keyname_level, value_level)
353	IF ( TRIM(value_leveltype) == TRIM(keyvalue_leveltype) .AND. &
354	TRIM(value_shortname) == TRIM(keyvalue_shortname) ) THEN
355
356	! If this is first 3d U wind message, get dimensions
357	IF (num_levels == 0) THEN
358	CONTINUE
359	CALL grib_get(igrib, keyname_xlon, xlon)
360	CALL grib_get(igrib, keyname_ylat, ylat)
361	ENDIF
362	!PRINT *, TRIM(value_shortname), ' ', TRIM(value_leveltype), ' ', value_level
363	num_levels = num_levels + 1
364	END IF
365	END IF
366	END DO
367
368
369	ELSE
370	PRINT *, "ERROR: class_gribfile::get_3d_u_dims(): problem opening GRIB file: ", filepath
371	STOP
372	END IF
373
374	CALL grib_close_file(ifile)
375
376	zlev = num_levels
377
378	END SUBROUTINE get_3d_u_dims
379
380
381
382	!!! Getter methods
383	INTEGER FUNCTION gribfile_num_xlon(gribfile_obj)
384
385	! Returns x (lon) dimension of met data
386	TYPE(gribfile_object), INTENT(IN) :: gribfile_obj
387
388	IF (.NOT. gribfile_obj%is_instantiated) THEN
389	PRINT *, 'ERROR: class_gribfile: gribfile_obj not instantiated'
390	ENDIF
391	gribfile_num_xlon = gribfile_obj%num_xlon
392
393	END FUNCTION gribfile_num_xlon
394
395	INTEGER FUNCTION gribfile_num_ylat(gribfile_obj)
396
397	! Returns y (lat) dimension of met data
398	TYPE(gribfile_object), INTENT(IN) :: gribfile_obj
399
400	IF (.NOT. gribfile_obj%is_instantiated) THEN
401	PRINT *, 'ERROR: class_gribfile: gribfile_obj not instantiated'
402	ENDIF
403	gribfile_num_ylat = gribfile_obj%num_ylat
404
405	END FUNCTION gribfile_num_ylat
406
407	INTEGER FUNCTION gribfile_num_zlevel(gribfile_obj)
408
409	! Returns z (level) dimension of met data
410	TYPE(gribfile_object), INTENT(IN) :: gribfile_obj
411
412	IF (.NOT. gribfile_obj%is_instantiated) THEN
413	PRINT *, 'ERROR: class_gribfile: gribfile_obj not instantiated'
414	ENDIF
415	gribfile_num_zlevel = gribfile_obj%num_zlevel
416
417	END FUNCTION gribfile_num_zlevel
418
419	END MODULE class_gribfile

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

Download in other formats:

Original Format