Commit 61e07bac authored by Espen Sollum's avatar Espen Sollum
Browse files

Adding files from CTBTO project wo_17

parent b5127f90
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
integer function detectformat()
!*****************************************************************************
! *
! This routine reads the 1st file with windfields to determine *
! the format. *
! *
! Authors: M. Harustak *
! *
! 6 May 2015 *
! *
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
! - Added routine to FP10 Flexpart distribution *
!*****************************************************************************
! *
! Variables: *
! fname file name of file to check *
! *
!*****************************************************************************
use par_mod
use com_mod
use class_gribfile
implicit none
character(len=255) :: filename
character(len=255) :: wfname1(maxwf)
integer :: metdata_format
! If no file is available
if ( maxwf.le.0 ) then
print*,'No wind file available'
detectformat = GRIBFILE_CENTRE_UNKNOWN
return
endif
! construct filename
filename = path(3)(1:length(3)) // trim(wfname(1))
! get format
detectformat = gribfile_centre(TRIM(filename))
end
MODULE class_gribfile
!*****************************************************************
! *
! This is a class-like structure for providing access to the *
! metadata and data within a GRIB met file *
! *
! *
!*****************************************************************
IMPLICIT NONE
PRIVATE ! The default is that everything is private, unless
! specified otherwise
! Note that all public interfaces and variables should have a
! GRIBFILE_ prefix
PUBLIC :: gribfile_object, &
gribfile_object_create, &
gribfile_printobj, &
gribfile_centre, &
gribfile_num_xlon, &
gribfile_num_ylat, &
gribfile_num_zlevel
PUBLIC :: GRIBFILE_TYPE_ECMWF_GRIB1, &
GRIBFILE_TYPE_ECMWF_GRIB2, &
GRIBFILE_TYPE_ECMWF_GRIB1_2, &
GRIBFILE_TYPE_NCEP_GRIB1, &
GRIBFILE_TYPE_NCEP_GRIB2, &
GRIBFILE_TYPE_UNKNOWN, &
GRIBFILE_CENTRE_NCEP, &
GRIBFILE_CENTRE_ECMWF, &
GRIBFILE_CENTRE_UNKNOWN
! These are codes for designating the type of GRIB file
! being looked at
INTEGER, PARAMETER :: GRIBFILE_TYPE_ECMWF_GRIB1 = 1, &
GRIBFILE_TYPE_ECMWF_GRIB2 = 2, &
GRIBFILE_TYPE_ECMWF_GRIB1_2 = 3, &
GRIBFILE_TYPE_NCEP_GRIB1 = 4, &
GRIBFILE_TYPE_NCEP_GRIB2 = 5, &
GRIBFILE_TYPE_UNKNOWN = -9999, &
GRIBFILE_CENTRE_NCEP = 1, &
GRIBFILE_CENTRE_ECMWF = 2, &
GRIBFILE_CENTRE_UNKNOWN = -9999
! These are the official centre codes for NCEP and ECMWF in grib files.
INTEGER, PARAMETER :: CENTRE_NCEP = 7, CENTRE_ECMWF = 98
TYPE gribfile_object
PRIVATE ! Make everything in here private so it's not directly manipulated outside
LOGICAL :: is_instantiated = .FALSE.
CHARACTER(LEN=256) :: file_path = ''
INTEGER :: grib_edition = 0 ! Not sure we want this, since it can vary on hybrid files
INTEGER :: grib_centre = GRIBFILE_CENTRE_UNKNOWN
INTEGER :: gribfile_type = GRIBFILE_TYPE_UNKNOWN
INTEGER :: num_xlon = -9999
INTEGER :: num_ylat = -9999
INTEGER :: num_zlevel = -9999
END TYPE gribfile_object
CONTAINS
SUBROUTINE gribfile_testhello()
PRINT *, 'Hello gribfile'
END SUBROUTINE gribfile_testhello
TYPE(gribfile_object) FUNCTION gribfile_object_create(filepath)
! This is the "constructor" for the pseudo gribfile object. Given the path to a gribfile,
! fills in attributes that can be accessed through methods in this
! module
USE grib_api
IMPLICIT NONE
CHARACTER(LEN=*), INTENT(IN) :: filepath ! full path to GRIB file
TYPE(gribfile_object) :: returned_object
INTEGER :: ifile, iret, igrib, grib_centre, gribfile_type
INTEGER :: xlon, ylat, zlev
CALL get_centre_and_type(filepath, grib_centre, gribfile_type)
returned_object%grib_centre = grib_centre
returned_object%gribfile_type = gribfile_type
! Get dimensions of 3d u field
CALL get_3d_u_dims(filepath, gribfile_type, xlon, ylat, zlev)
returned_object%num_xlon = xlon
returned_object%num_ylat = ylat
returned_object%num_zlevel = zlev
returned_object%is_instantiated = .TRUE.
returned_object%file_path = TRIM(filepath)
gribfile_object_create = returned_object
END FUNCTION gribfile_object_create
SUBROUTINE gribfile_printobj(gribfile_obj)
! Pretty prints the attributes of the gribfile pseudo-object
TYPE(gribfile_object), INTENT(IN) :: gribfile_obj
PRINT *, 'is_instantiated: ', gribfile_obj%is_instantiated
PRINT *, 'filepath: ', TRIM(gribfile_obj%file_path)
PRINT *, 'grib_centre: ', gribfile_obj%grib_centre
PRINT *, 'gribfile_type: ', gribfile_obj%gribfile_type
PRINT *, 'num_xlon: ', gribfile_obj%num_xlon
PRINT *, 'num_ylat: ', gribfile_obj%num_ylat
PRINT *, 'num_zlevel: ', gribfile_obj%num_zlevel
END SUBROUTINE gribfile_printobj
INTEGER FUNCTION gribfile_centre(filepath)
! Returns an integer constant denoting the grib centre (currently either ECMWF, NCEP or UNKNOWN)
! for the specified filepath. Returns one of the GRIBFILE_CENTRE_ constants defined at top of this
! module.
USE grib_api
IMPLICIT NONE
CHARACTER(LEN=*), INTENT(IN) :: filepath ! full path to GRIB file
INTEGER :: ifile, iret, igrib, grib_centre
CALL grib_open_file(ifile, filepath, 'r', iret)
IF (iret == 0) THEN
! Use first record to detect centre, which is assumed constant
! amongst all messages
CALL grib_new_from_file(ifile, igrib, iret)
CALL grib_get(igrib, 'centre', grib_centre)
CALL grib_close_file(ifile)
ELSE
PRINT *, "WARNING: problem opening GRIB file: ", filepath
grib_centre = -999
END IF
IF (grib_centre == CENTRE_NCEP) THEN
gribfile_centre = GRIBFILE_CENTRE_NCEP
ELSE IF (grib_centre == CENTRE_ECMWF) THEN
gribfile_centre = GRIBFILE_CENTRE_ECMWF
ELSE
gribfile_centre = GRIBFILE_CENTRE_UNKNOWN
END IF
END FUNCTION gribfile_centre
! This is currently a PRIVATE subroutine
SUBROUTINE get_centre_and_type(filepath, grib_centre, gribfile_type)
! Given specified grib file, passes back centre and gribfile
! type to the calling program. Numeric codes are defined as integer parameters
! in this module
! To get this information, we have to iterate through the entire file in order to
! determine if it is hybrid or not
!
USE grib_api
IMPLICIT NONE
CHARACTER(LEN=*), INTENT(IN) :: filepath ! full path to GRIB file
INTEGER, INTENT(OUT) :: grib_centre, gribfile_type
INTEGER :: fileptr, iret, igrib, centre, grib_edition
LOGICAL :: end_of_file
LOGICAL :: grib1_detected, grib2_detected
grib1_detected = .FALSE.
grib2_detected = .FALSE.
CALL grib_open_file(fileptr, filepath, 'r', iret)
IF (iret /= 0) THEN
PRINT *, 'class_gributils:get_centre_and_type()...'
PRINT *, ' unable to open filepath: ', filepath
STOP
END IF
! Use first record to detect centre and and grib version of first messages. We will
! then assume that all following messages have same centre, but not necessarily same
! GRIB version
CALL grib_new_from_file(fileptr, igrib, iret)
CALL grib_get(igrib, 'centre', grib_centre)
CALL grib_get(igrib, 'edition', grib_edition)
IF (grib_edition == 1) grib1_detected = .TRUE.
IF (grib_edition == 2) grib2_detected = .TRUE.
! Now, iterate through the rest of records to determine if this is a mixed edition file
end_of_file = .FALSE.
DO WHILE (.NOT. end_of_file)
CALL grib_new_from_file(fileptr, igrib, iret)
IF (iret .eq. GRIB_END_OF_FILE) THEN
end_of_file = .TRUE.
ELSE
! Get edition from file
CALL grib_get(igrib, 'edition', grib_edition)
IF (grib_edition .eq. 1) grib1_detected = .TRUE.
IF (grib_edition .eq. 2) grib2_detected = .TRUE.
END IF
END DO
CALL grib_close_file(fileptr)
! Determine the gribfile type depending on centre and edition(s)
IF (grib_centre == CENTRE_ECMWF) THEN
IF (grib1_detected .AND. grib2_detected) THEN
gribfile_type = GRIBFILE_TYPE_ECMWF_GRIB1_2
ELSE IF (grib1_detected .AND. .NOT. grib2_detected) THEN
gribfile_type = GRIBFILE_TYPE_ECMWF_GRIB1
ELSE IF (.NOT. grib1_detected .AND. grib2_detected) THEN
gribfile_type = GRIBFILE_TYPE_ECMWF_GRIB2
ELSE
gribfile_type = GRIBFILE_TYPE_UNKNOWN
ENDIF
ELSE IF (grib_centre == CENTRE_NCEP) THEN
IF (grib1_detected .AND. .NOT. grib2_detected) THEN
gribfile_type = GRIBFILE_TYPE_NCEP_GRIB1
ELSE IF (.NOT. grib1_detected .AND. grib2_detected) THEN
gribfile_type = GRIBFILE_TYPE_NCEP_GRIB2
ELSE
gribfile_type = GRIBFILE_TYPE_UNKNOWN
ENDIF
ELSE
gribfile_type = GRIBFILE_TYPE_UNKNOWN
ENDIF
END SUBROUTINE get_centre_and_type
SUBROUTINE get_3d_u_dims(filepath, gribfile_type, xlon, ylat, zlev)
! Looks at the 3d u fields in the GRIBFILE to get x and y dims, as well as number of levels
USE grib_api
IMPLICIT NONE
CHARACTER(LEN=*), INTENT(IN) :: filepath ! full path to GRIB file
INTEGER, INTENT(IN) :: gribfile_type
INTEGER, INTENT(OUT) :: xlon, ylat, zlev
INTEGER :: ifile, iret, igrib, grib_centre
LOGICAL :: end_of_file
! These will be assigned according to type of grib file, then used to filter
! for the 3d messages
! Name of the key being sought
CHARACTER(LEN=64) :: keyname_leveltype, keyname_shortname, keyname_level, &
keyname_xlon, keyname_ylat
! The key value being filtered for
CHARACTER(LEN=64) :: keyvalue_leveltype, keyvalue_shortname
! Actual values read in from the grib file
CHARACTER(LEN=64) :: value_leveltype, value_shortname
INTEGER :: value_level
INTEGER :: num_levels
! Get the field names to read, based on the type of grib file
!!! Note that ALL of these have the same key names, except that
!!! leveltype is 'hybrid' in ECMWF and 'isobaricInhPa' in NCEP.
!!! This could probably be consolidated, but because these are
!!! files that go through some preprocessing, and aren't
!!! necessarily standard (at least for ECMWF), I'm going to be
!!! safe and leave it as is for now, so it would be easier to
!!! modify for one type, if necessary.
IF (gribfile_type == GRIBFILE_TYPE_ECMWF_GRIB1) THEN
keyname_leveltype = 'typeOfLevel'
keyname_shortname = 'shortName'
keyname_level = 'level'
keyvalue_leveltype = 'hybrid'
keyvalue_shortname = 'u'
keyname_xlon = 'Ni'
keyname_ylat = 'Nj'
ELSE IF (gribfile_type == GRIBFILE_TYPE_ECMWF_GRIB1_2) THEN
keyname_leveltype = 'typeOfLevel'
keyname_shortname = 'shortName'
keyname_level = 'level'
keyvalue_leveltype = 'hybrid'
keyvalue_shortname = 'u'
keyname_xlon = 'Ni'
keyname_ylat = 'Nj'
ELSE IF (gribfile_type == GRIBFILE_TYPE_ECMWF_GRIB2) THEN
keyname_leveltype = 'typeOfLevel'
keyname_shortname = 'shortName'
keyname_level = 'level'
keyvalue_leveltype = 'hybrid'
keyvalue_shortname = 'u'
keyname_xlon = 'Ni'
keyname_ylat = 'Nj'
ELSE IF (gribfile_type == GRIBFILE_TYPE_NCEP_GRIB1) THEN
keyname_leveltype = 'typeOfLevel'
keyname_shortname = 'shortName'
keyname_level = 'level'
keyvalue_leveltype = 'isobaricInhPa'
keyvalue_shortname = 'u'
keyname_xlon = 'Ni'
keyname_ylat = 'Nj'
ELSE IF (gribfile_type == GRIBFILE_TYPE_NCEP_GRIB2) THEN
keyname_leveltype = 'typeOfLevel'
keyname_shortname = 'shortName'
keyname_level = 'level'
keyvalue_leveltype = 'isobaricInhPa'
keyvalue_shortname = 'u'
keyname_xlon = 'Ni'
keyname_ylat = 'Nj'
ELSE
PRINT *, 'class_gribfile:get_3d_u_dims(): Unsupported gribfile type: ', gribfile_type
STOP
ENDIF
CALL grib_open_file(ifile, filepath, 'r', iret)
IF (iret == 0) THEN
! Iterate through all messages to count 3d u messages (levels) and get x,y dimensions
end_of_file = .FALSE.
num_levels = 0
DO WHILE (.NOT. end_of_file)
CALL grib_new_from_file(ifile, igrib, iret)
IF (iret .eq. GRIB_END_OF_FILE) THEN
end_of_file = .TRUE.
ELSE
! Get relevant keys and filter for the 3d U wind
CALL grib_get(igrib, keyname_shortname, value_shortname)
CALL grib_get(igrib, keyname_leveltype, value_leveltype)
CALL grib_get(igrib, keyname_level, value_level)
IF ( TRIM(value_leveltype) == TRIM(keyvalue_leveltype) .AND. &
TRIM(value_shortname) == TRIM(keyvalue_shortname) ) THEN
! If this is first 3d U wind message, get dimensions
IF (num_levels == 0) THEN
CONTINUE
CALL grib_get(igrib, keyname_xlon, xlon)
CALL grib_get(igrib, keyname_ylat, ylat)
ENDIF
!PRINT *, TRIM(value_shortname), ' ', TRIM(value_leveltype), ' ', value_level
num_levels = num_levels + 1
END IF
END IF
END DO
ELSE
PRINT *, "ERROR: class_gribfile::get_3d_u_dims(): problem opening GRIB file: ", filepath
STOP
END IF
CALL grib_close_file(ifile)
zlev = num_levels
END SUBROUTINE get_3d_u_dims
!!! Getter methods
INTEGER FUNCTION gribfile_num_xlon(gribfile_obj)
! Returns x (lon) dimension of met data
TYPE(gribfile_object), INTENT(IN) :: gribfile_obj
IF (.NOT. gribfile_obj%is_instantiated) THEN
PRINT *, 'ERROR: class_gribfile: gribfile_obj not instantiated'
ENDIF
gribfile_num_xlon = gribfile_obj%num_xlon
END FUNCTION gribfile_num_xlon
INTEGER FUNCTION gribfile_num_ylat(gribfile_obj)
! Returns y (lat) dimension of met data
TYPE(gribfile_object), INTENT(IN) :: gribfile_obj
IF (.NOT. gribfile_obj%is_instantiated) THEN
PRINT *, 'ERROR: class_gribfile: gribfile_obj not instantiated'
ENDIF
gribfile_num_ylat = gribfile_obj%num_ylat
END FUNCTION gribfile_num_ylat
INTEGER FUNCTION gribfile_num_zlevel(gribfile_obj)
! Returns z (level) dimension of met data
TYPE(gribfile_object), INTENT(IN) :: gribfile_obj
IF (.NOT. gribfile_obj%is_instantiated) THEN
PRINT *, 'ERROR: class_gribfile: gribfile_obj not instantiated'
ENDIF
gribfile_num_zlevel = gribfile_obj%num_zlevel
END FUNCTION gribfile_num_zlevel
END MODULE class_gribfile
FC = gfortran
SRCDIR = ..
BINARY = testdrive
GRIBAPI = /opt/grib-api
MODOBJS = class_gribfile_mod.o
MODULES = class_gribfile_mod.mod
VPATH = ${SRCDIR}
FFLAGS = -mcmodel=medium
INCLUDES = -I${GRIBAPI}/include
LDFLAGS = -L${GRIBAPI}/lib -lgrib_api_f90 -lgrib_api
${BINARY} : ${BINARY}.o ${MODOBJS}
${FC} -o ${BINARY} ${BINARY}.o ${MODOBJS} ${LDFLAGS}
${BINARY}.o : ${BINARY}.f90 ${MODULES}
class_gribfile_mod.mod : class_gribfile_mod.o Makefile
${FC} -c ../class_gribfile_mod.f90 ${FFLAGS} ${INCLUDES}
test : ${BINARY}
( ./${BINARY} )
#------------------------------------------------------------------
clean :
rm -r *.o *.mod ${BINARY}
%.o: %.f90
${FC} -c ${FFLAGS} ${INCLUDES} $<
PROGRAM testdrive
USE class_gribfile
IMPLICIT NONE
CHARACTER(LEN=256) :: the_file_path
INTEGER, PARAMETER :: NCASES = 6
CHARACTER(LEN=256), DIMENSION(NCASES) :: file_paths
CHARACTER(LEN=256), DIMENSION(NCASES) :: descriptions
TYPE(gribfile_object) :: my_gribfile
INTEGER :: case_number
!!!!!!!!!!!!!!!!!!!
! These define the header and the grib file used for each test
descriptions(1) = "ECMWF GRIB1 on global 1.0 degree domain"
file_paths(1) = "../../../devtest/case_data/met_data/ecmwf/t1_03h_ec1p0d/EN13062503"
descriptions(2) = "ECMWF GRIB1/2 on tiny domain"
file_paths(2) = "../../../devtest/case_data/met_data/ecmwf/t1_33h_ec1p0d/EL14091909"
descriptions(3) = "ECMWF GRIB1/2 on global 1.0 degree domain"
file_paths(3) = "../../../devtest/case_data/met_data/ecmwf/t1_03h_ec1p0d_grib1-2/EE13110700"
descriptions(4) = "ECMWF GRIB2 on global 1.0 degree domain"
file_paths(4) = "../../../devtest/case_data/met_data/ecmwf/t1_03h_ec1p0d_grib2/EN13110700"
descriptions(5) = "NCEP GRIB1 on global 1.0 degree domain"
file_paths(5) = "../../../devtest/case_data/met_data/ncep/t1_06h_nc1p0d_grib1/GD05051406"
descriptions(6) = "NCEP GRIB2 on global 1.0 degree domain"
file_paths(6) = "../../../devtest/case_data/met_data/ncep/t1_03h_nc1p0d/GF15021603"
!!!!!!!!!!!!!!!!!!!
DO case_number = 1,NCASES
PRINT *,
PRINT *, TRIM( descriptions( case_number) )
my_gribfile = gribfile_object_create( file_paths(case_number) )
PRINT *,
PRINT *, 'Output from calling gribfile_printobj()...'
PRINT *, '++++++++++'
CALL gribfile_printobj(my_gribfile)
PRINT *, '++++++++++'
PRINT *,
PRINT *, 'Output from the getter methods...'
PRINT *, 'gribfile_center() : ', gribfile_centre( file_paths(case_number) )
PRINT *, 'gribfile_num_xlon : ', gribfile_num_xlon(my_gribfile)
PRINT *, 'gribfile_num_ylat : ', gribfile_num_ylat(my_gribfile)
PRINT *, 'gribfile_num_zlevel : ', gribfile_num_zlevel(my_gribfile)
PRINT *,
PRINT *, '----------------------------'
ENDDO
END PROGRAM testdrive
This diff is collapsed.
......@@ -11,31 +11,40 @@ SHELL = /bin/bash
# To use gfortran version 4.9, add "gcc=4.9" to the make command, e.g.
# 'make -j ecmwf gcc=4.9',
# also set environment variable LD_LIBRARY_PATH to point to compiler libraries