source: flexpart.git/flexpart_code/class_vtable_mod.F90 @ 7e6dc50

FPv9.3.1FPv9.3.2fp9.3.1-20161214-nc4grib2nc4_repair
Last change on this file since 7e6dc50 was 7e6dc50, checked in by Don Morton <Don.Morton@…>, 8 years ago

Removed a fair amount of annoying print statements.

  • Property mode set to 100644
File size: 17.3 KB
Line 
1
2module class_vtable
3
4    implicit none
5    private
6
7    ! Note that all public interfaces and variables should have a
8    ! "vtable" prefix
9    public :: Vtable,                 &
10              Vtable_record,          &
11              vtable_load_by_name,    &
12              vtable_get_fpname,      &
13              vtable_detect_gribfile_type, &
14              Vtable_GRIBFILE_TYPE_ECMWF_GRIB1,   &
15              Vtable_GRIBFILE_TYPE_ECMWF_GRIB2,   &
16              Vtable_GRIBFILE_TYPE_ECMWF_GRIB1_2, &
17              Vtable_GRIBFILE_TYPE_NCEP_GRIB1,    &
18              Vtable_GRIBFILE_TYPE_NCEP_GRIB2,    &
19              Vtable_GRIBFILE_TYPE_UNKNOWN   
20
21
22    integer, parameter :: VTABLE_MISSING_ENTRY = -9999
23
24    ! These are codes for designating the type of GRIB file being
25    ! looked at
26    integer, parameter :: Vtable_GRIBFILE_TYPE_ECMWF_GRIB1 = 1,   &
27                          Vtable_GRIBFILE_TYPE_ECMWF_GRIB2 = 2,   &
28                          Vtable_GRIBFILE_TYPE_ECMWF_GRIB1_2 = 3, &
29                          Vtable_GRIBFILE_TYPE_NCEP_GRIB1 = 4,    &
30                          Vtable_GRIBFILE_TYPE_NCEP_GRIB2 = 5,    &
31                          Vtable_GRIBFILE_TYPE_UNKNOWN = -99 
32                         
33    ! These codes depict the origin/model of the GRIB File
34    integer, parameter :: GRIB_CENTRE_NCEP = 7, &
35                          GRIB_CENTRE_ECMWF = 98
36
37    type Vtable_record
38        character(len=15) :: fpname
39        integer :: paramid
40        integer :: indicator_of_parameter
41        integer :: discipline
42        integer :: category
43        integer :: number
44        integer :: typesurface
45        character(len=25) :: typelevel
46        character(len=15) :: units
47        character(len=10) :: shortname
48        character(len=25) :: description
49        integer :: grib_version
50        character(len=10) :: center
51    end type Vtable_record
52
53    type Vtable
54        logical :: initialized=.FALSE.
55        character(len=255) :: path_to_vtable_file
56        integer :: num_entries = 0
57        type(Vtable_record), allocatable :: the_entries(:)
58    end type Vtable
59
60contains
61
62   
63
64
65    integer function vtable_detect_gribfile_type(gribfilename)
66        ! Given specified grib file, returns an integer code indicating its
67        ! type to the calling program.  Numeric codes are defined as integer parameters
68        ! in this module
69
70        use grib_api       
71        implicit none
72        character(len=255), intent(in) :: gribfilename  ! Full path to grib file
73       
74        integer :: ifile, iret, igrib, grib_centre, grib_edition
75        logical :: end_of_file
76        logical :: grib1_detected = .FALSE., grib2_detected = .FALSE.
77       
78        call grib_open_file(ifile, gribfilename, 'r', iret)
79
80        ! Use first record to detect centre and and grib version of first messages.  We will
81        ! then assume that all following messages have same centre, but not necessarily same
82        ! GRIB version
83        call grib_new_from_file(ifile, igrib, iret)
84        call grib_get(igrib, 'centre', grib_centre)
85        call grib_get(igrib, 'edition', grib_edition)
86       
87        if (grib_edition .eq. 1) grib1_detected = .TRUE.
88        if (grib_edition .eq. 2) grib2_detected = .TRUE.
89       
90        ! Now, iterate through the rest of records to determine if this is a mixed edition file
91        end_of_file = .FALSE.
92        do while (.NOT. end_of_file)   
93            call grib_new_from_file(ifile, igrib, iret)
94            if (iret .eq. GRIB_END_OF_FILE) then
95                end_of_file = .TRUE.
96            else
97
98                ! Get edition from file
99                call grib_get(igrib, 'edition', grib_edition)
100                if (grib_edition .eq. 1) grib1_detected = .TRUE.
101                if (grib_edition .eq. 2) grib2_detected = .TRUE.
102            end if
103        end do
104       
105    call grib_close_file(ifile)
106       
107    ! Determine the gribfile type depending on centre and edition(s)
108   
109    if (grib_centre .eq. GRIB_CENTRE_ECMWF) then
110        if (grib1_detected .and. grib2_detected) then
111            vtable_detect_gribfile_type = Vtable_GRIBFILE_TYPE_ECMWF_GRIB1_2
112        else if (grib1_detected .and. .not. grib2_detected) then
113            vtable_detect_gribfile_type = Vtable_GRIBFILE_TYPE_ECMWF_GRIB1
114        else if (.not. grib1_detected .and. grib2_detected) then
115            vtable_detect_gribfile_type = Vtable_GRIBFILE_TYPE_ECMWF_GRIB2
116        else
117            vtable_detect_gribfile_type = Vtable_GRIBFILE_TYPE_UNKNOWN                       
118        endif
119    else if (grib_centre .eq. GRIB_CENTRE_NCEP) then
120        if (grib1_detected .and. .not. grib2_detected) then
121             vtable_detect_gribfile_type = Vtable_GRIBFILE_TYPE_NCEP_GRIB1       
122        else if (.not. grib1_detected .and. grib2_detected) then
123            vtable_detect_gribfile_type = Vtable_GRIBFILE_TYPE_NCEP_GRIB2
124        else
125            vtable_detect_gribfile_type = Vtable_GRIBFILE_TYPE_UNKNOWN                       
126        endif
127    else
128        vtable_detect_gribfile_type = Vtable_GRIBFILE_TYPE_UNKNOWN
129    endif
130
131    end function vtable_detect_gribfile_type
132
133
134
135    subroutine vtable_load_by_name(vtable_name, the_vtable_data)
136        implicit none
137
138        character(len=255), intent(in) :: vtable_name  ! Full path to vtable file
139
140        logical :: lexist
141        integer :: ierr
142        integer :: num_vrecs = 0
143        integer :: vrec_idx
144        character(len=255) :: file_line = ' '
145
146        type(Vtable), intent(out) :: the_vtable_data    ! Data structure holding the vtable
147
148        type(Vtable_record) :: vrecord
149
150        ! Make sure the file exists
151        inquire(file=trim(vtable_name), exist=lexist)
152        if (.not. lexist) then
153            print *, 'file: ', trim(vtable_name), ' does not exist...'
154            stop
155        endif
156
157        ! Open file
158        open(10, file=trim(vtable_name), status='old', form='formatted', iostat=ierr)
159        if (ierr .ne. 0) then
160            print *, 'file: ', trim(vtable_name), ' open failed...'
161            stop
162        endif
163
164        ! Go through the file once and count the vtable_records
165        ! Read past headers
166        file_line = ' '
167        do while(file_line(1:1) .ne. '-')
168            read(10, '(A255)', iostat=ierr) file_line
169        enddo
170
171        ! Now we are at the '----------' line - process everything between
172        ! here and the next '----------' line.  In this case, we just want to
173        ! count
174        file_line = ' '
175        num_vrecs = 0
176        do while(file_line(1:1) .ne. '-')
177            read(10, '(A255)', iostat=ierr) file_line
178            !print *, file_line
179            num_vrecs = num_vrecs + 1
180        enddo
181        num_vrecs = num_vrecs - 1
182
183        !print *, 'num_vrecs: ', num_vrecs
184
185        ! Rewind
186        rewind(unit=10)
187
188        ! Allocate array for storing the vtable records, and store
189        ! num_entries
190        !print *, 'Ready to allocate the_vtable_data'
191        allocate(the_vtable_data%the_entries(num_vrecs))
192        !print *, 'Allocated the_vtable_data'
193        the_vtable_data%num_entries = num_vrecs
194
195        ! Read, parse and store the vtable records
196        ! Read past headers
197        file_line = ' '
198        do while(file_line(1:1) .ne. '-')
199            read(10, '(A255)', iostat=ierr) file_line
200            !print *, file_line
201        enddo
202
203        ! Now we are at the '----------' line - process everything between
204        ! here and the next '----------' line.  In this case, we just want to
205        ! count
206        file_line = ' '
207        vrec_idx = 0
208        do while(file_line(1:1) .ne. '-')
209            read(10, '(A255)', iostat=ierr) file_line
210            if (file_line(1:1) .ne. '-') then
211                ! PROCESS THE LINE
212                vrec_idx = vrec_idx + 1
213
214                ! Parse the line and put it in the vtable structure
215                the_vtable_data%the_entries(vrec_idx) = vtable_parse_record(file_line)
216             
217                !print *, the_vtable_data%the_entries(vrec_idx)
218                !print *, file_line
219                !print *, 'hello'
220            endif
221        enddo
222        num_vrecs = num_vrecs - 1
223
224
225        ! Close the file
226        close(unit=10)
227
228        the_vtable_data%initialized = .TRUE.
229       
230        !print *, the_vtable_data%the_entries(1)
231    end subroutine vtable_load_by_name
232
233
234
235
236    type(Vtable_record) function vtable_parse_record(vtable_line)
237
238
239    !!! Using a vtable line as input argument, parses into a Vtable_record, and returns that
240    !!! record
241        implicit none
242        character(LEN=255), intent(in) :: vtable_line
243
244        !!! This is a sample of what a Vtable line will look like
245        !!! THIS NEEDS TO BE MODIFIED FOR NEW STRUCTURE (23 Nov 2015)
246        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
247        !vtable_line = 'UU       |   131   |     0      |    2     |    2   | &
248        !           &   105      |  hybrid            | ms**-1    |    u    &
249        !           &    | u wind             |     2      |'
250        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
251
252
253        ! Storage for Vtable record tokens
254        character(25) :: token_fpname='',&
255                         token_paramid='', &
256                         token_indofparam='', &
257                         token_discipline='', &
258                         token_category='', &
259                         token_number='', &
260                         token_typesurface='', &
261                         token_typelevel='', &
262                         token_units='', &
263                         token_shortname='', &
264                         token_description='', &
265                         token_gribversion='', &
266                         token_center=''
267
268        ! These indices mark the locations of the '|' delimiter in a Vtable record
269        integer :: delim_fpname_idx, &
270                   delim_paramid_idx, &
271                   delim_indofparam_idx, &
272                   delim_disc_idx, &
273                   delim_cat_idx, &
274                   delim_numb_idx, &
275                   delim_typesurf_idx, &
276                   delim_typelevel_idx, &
277                   delim_units_idx, &
278                   delim_shortname_idx, &
279                   delim_descr_idx, &
280                   delim_version_idx, &
281                   delim_center_idx
282
283        type(Vtable_record) :: vrecord
284
285        integer :: istat    ! Error indicator for some I/O routines
286
287        ! Calculate the indices of each field so we can extract later
288        delim_fpname_idx = index(vtable_line, '|')
289        delim_paramid_idx = index(vtable_line(delim_fpname_idx+1:), '|') &
290                               + delim_fpname_idx
291        delim_indofparam_idx = index(vtable_line(delim_paramid_idx+1:), '|') &
292                               + delim_paramid_idx
293        delim_disc_idx = index(vtable_line(delim_indofparam_idx+1:), '|') &
294                            + delim_indofparam_idx
295        delim_cat_idx = index(vtable_line(delim_disc_idx+1:), '|') &
296                           + delim_disc_idx
297        delim_numb_idx = index(vtable_line(delim_cat_idx+1:), '|') &
298                            + delim_cat_idx
299        delim_typesurf_idx = index(vtable_line(delim_numb_idx+1:), '|') &
300                            + delim_numb_idx
301        delim_typelevel_idx = index(vtable_line(delim_typesurf_idx+1:), '|') &
302                            + delim_typesurf_idx
303        delim_units_idx = index(vtable_line(delim_typelevel_idx+1:), '|') &
304                            + delim_typelevel_idx
305        delim_shortname_idx = index(vtable_line(delim_units_idx+1:), '|') &
306                            + delim_units_idx
307        delim_descr_idx = index(vtable_line(delim_shortname_idx+1:), '|') &
308                            + delim_shortname_idx
309        delim_version_idx = index(vtable_line(delim_descr_idx+1:), '|') &
310                            + delim_descr_idx
311        delim_center_idx = index(vtable_line(delim_version_idx+1:), '|') &
312                            + delim_version_idx
313
314        ! Extract the tokens
315        token_fpname = vtable_line(:delim_fpname_idx-1)
316        token_paramid = vtable_line(delim_fpname_idx+1:delim_paramid_idx-1)
317        token_indofparam = vtable_line(delim_paramid_idx+1:delim_indofparam_idx-1)
318        token_discipline = vtable_line(delim_indofparam_idx+1:delim_disc_idx-1)
319        token_category = vtable_line(delim_disc_idx+1:delim_cat_idx-1)
320        token_number = vtable_line(delim_cat_idx+1:delim_numb_idx-1)
321        token_typesurface = vtable_line(delim_numb_idx+1:delim_typesurf_idx-1)
322        token_typelevel = vtable_line(delim_typesurf_idx+1:delim_typelevel_idx-1)
323        token_units = vtable_line(delim_typelevel_idx+1:delim_units_idx-1)
324        token_shortname = vtable_line(delim_units_idx+1:delim_shortname_idx-1)
325        token_description = vtable_line(delim_shortname_idx+1:delim_descr_idx-1)
326        token_gribversion = vtable_line(delim_descr_idx+1:delim_version_idx-1)
327        token_center = vtable_line(delim_version_idx+1:delim_center_idx-1)
328
329        ! Jam the data in the record for return
330        vrecord%fpname = token_fpname
331
332        read(token_paramid, *, iostat=istat) vrecord%paramid
333        if (istat .ne. 0) vrecord%paramid = VTABLE_MISSING_ENTRY
334
335        read(token_indofparam, *, iostat=istat) vrecord%indicator_of_parameter
336        if (istat .ne. 0) vrecord%indicator_of_parameter = VTABLE_MISSING_ENTRY
337
338
339        read(token_discipline, *, iostat=istat) vrecord%discipline
340        if (istat .ne. 0) vrecord%discipline = VTABLE_MISSING_ENTRY
341
342        read(token_category, *, iostat=istat) vrecord%category
343        if (istat .ne. 0) vrecord%category = VTABLE_MISSING_ENTRY
344
345        read(token_number, *, iostat=istat) vrecord%number
346        if (istat .ne. 0) vrecord%number = VTABLE_MISSING_ENTRY
347
348        read(token_typesurface, *, iostat=istat) vrecord%typesurface
349        if (istat .ne. 0) vrecord%typesurface = VTABLE_MISSING_ENTRY
350
351        vrecord%typelevel = token_typelevel
352        vrecord%units = token_units
353        vrecord%shortname = token_shortname
354        vrecord%description = token_description
355
356        read(token_gribversion, *, iostat=istat) vrecord%grib_version
357        if (istat .ne. 0) vrecord%grib_version = VTABLE_MISSING_ENTRY
358
359        vrecord%center = token_center
360
361
362        !vrecord%fpname = 'UU'
363        vtable_parse_record = vrecord
364
365        !print *, "Hello vtable_parse_record()"
366        !print *, vrecord
367    end function vtable_parse_record
368
369    character(len=15) function vtable_get_fpname(igrib, vtable_object)
370   
371        !!! Assumes that a calling routine has opened up a GRIB file and obtained the
372        !!! grib id for a specific message.
373        !!! Given a grib message and a Vtable, looks up the message parameters in the Vtable
374        !!! and, if found, returns the fpname
375       
376        use grib_api
377        implicit none
378       
379        integer, intent(in) :: igrib
380        type(Vtable), intent(in) :: vtable_object
381       
382        integer :: parameter_id, category, number, discipline, edition, surface_type, &
383                   level, indicator_of_parameter
384        character(len=10) :: center
385       
386        integer :: idx
387        logical :: record_match
388       
389        call grib_get(igrib, 'editionNumber', edition)
390        call grib_get(igrib, 'level', level)
391               
392        if (edition .eq. 1) then
393            call grib_get(igrib, 'indicatorOfParameter', indicator_of_parameter)
394            call grib_get(igrib, 'indicatorOfTypeOfLevel', surface_type)
395            !print *, '(edition, indicator_of_parameter, surftype, level): ', edition, indicator_of_parameter, surface_type,&
396            !          level
397        else if (edition .eq. 2) then
398            call grib_get(igrib, 'parameterNumber', number)
399            call grib_get(igrib, 'parameterCategory', category)
400            call grib_get(igrib, 'discipline', discipline)
401            call grib_get(igrib, 'typeOfFirstFixedSurface', surface_type)
402            !print *, '(edition, number, cat, disc, surftype, level): ', edition, number, &
403            !          category, discipline, surface_type, level
404        else
405            print *, 'Illegal edition: ', edition
406            stop
407        endif       
408
409        ! Iterate through Vtable and look for a match
410        vtable_get_fpname = 'None'
411        record_match = .FALSE.
412        idx = 1
413        do while (.NOT. record_match .AND. idx .LE. vtable_object%num_entries)
414
415
416
417            if (edition .eq. 1) then
418                if (vtable_object%the_entries(idx)%indicator_of_parameter .eq. indicator_of_parameter .and. &
419                    vtable_object%the_entries(idx)%typesurface .eq. surface_type) then
420                    vtable_get_fpname = vtable_object%the_entries(idx)%fpname
421                    record_match = .TRUE.
422                end if                               
423            else if (edition .eq. 2) then
424                if (vtable_object%the_entries(idx)%discipline .eq. discipline .and.    &
425                    vtable_object%the_entries(idx)%number .eq. number .and.   &
426                    vtable_object%the_entries(idx)%category .eq. category .and.   &
427                    vtable_object%the_entries(idx)%typesurface .eq. surface_type) then
428               
429                    vtable_get_fpname = vtable_object%the_entries(idx)%fpname
430                    record_match = .TRUE.
431                end if
432            else
433                print *, 'Illegal edition: ', edition
434                stop
435            endif                 
436            idx = idx + 1   
437        end do
438       
439       
440    end function vtable_get_fpname
441
442   
443
444end module class_vtable
445
446
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG