Index: src/detectformat.f90
===================================================================
--- src/detectformat.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/detectformat.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,67 @@
+!**********************************************************************
+! 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 
Index: src/gributils/class_gribfile_mod.f90
===================================================================
--- src/gributils/class_gribfile_mod.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/gributils/class_gribfile_mod.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,419 @@
+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
Index: src/gributils/test/.gitignore
===================================================================
--- src/gributils/test/.gitignore	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/gributils/test/.gitignore	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,3 @@
+*.mod 
+*.o
+testdrive
Index: src/gributils/test/Makefile
===================================================================
--- src/gributils/test/Makefile	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/gributils/test/Makefile	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,39 @@
+
+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} $<
+
+
Index: src/gributils/test/testdrive.f90
===================================================================
--- src/gributils/test/testdrive.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/gributils/test/testdrive.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,61 @@
+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
Index: src/gridcheck_ecmwf.f90
===================================================================
--- src/gridcheck_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/gridcheck_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,587 @@
+!**********************************************************************
+! 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/>.   *
+!**********************************************************************
+
+subroutine gridcheck_ecmwf
+
+  !**********************************************************************
+  !                                                                     *
+  !             FLEXPART MODEL SUBROUTINE GRIDCHECK                     *
+  !                                                                     *
+  !**********************************************************************
+  !                                                                     *
+  !             AUTHOR:      G. WOTAWA                                  *
+  !             DATE:        1997-08-06                                 *
+  !             LAST UPDATE: 1997-10-10                                 *
+  !                                                                     *
+  !             Update:      1999-02-08, global fields allowed, A. Stohl*
+  !             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
+  !                                 ECMWF grib_api                      *
+  !             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
+  !                                 ECMWF grib_api                      *
+  !                                                                     *
+  !   Unified ECMWF and GFS builds                                      *
+  !   Marian Harustak, 12.5.2017                                        *
+  !     - Renamed from gridcheck to gridcheck_ecmwf                     *
+  !                                                                     *
+  !**********************************************************************
+  !                                                                     *
+  ! DESCRIPTION:                                                        *
+  !                                                                     *
+  ! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT      *
+  ! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST-   *
+  ! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE       *
+  ! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES    *
+  ! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT    *
+  ! ANY CALL.                                                           *
+  !                                                                     *
+  ! XLON0                geographical longitude of lower left gridpoint *
+  ! YLAT0                geographical latitude of lower left gridpoint  *
+  ! NX                   number of grid points x-direction              *
+  ! NY                   number of grid points y-direction              *
+  ! DX                   grid distance x-direction                      *
+  ! DY                   grid distance y-direction                      *
+  ! NUVZ                 number of grid points for horizontal wind      *
+  !                      components in z direction                      *
+  ! NWZ                  number of grid points for vertical wind        *
+  !                      component in z direction                       *
+  ! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
+  !                      points of the polar stereographic grid):       *
+  !                      used to check the CFL criterion                *
+  ! UVHEIGHT(1)-         heights of gridpoints where u and v are        *
+  ! UVHEIGHT(NUVZ)       given                                          *
+  ! WHEIGHT(1)-          heights of gridpoints where w is given         *
+  ! WHEIGHT(NWZ)                                                        *
+  !                                                                     *
+  !**********************************************************************
+
+  use grib_api
+  use par_mod
+  use com_mod
+  use conv_mod
+  use cmapf_mod, only: stlmbr,stcm2p
+
+  implicit none
+
+  !HSO  parameters for grib_api
+  integer :: ifile
+  integer :: iret
+  integer :: igrib
+  integer :: gotGrid
+  real(kind=4) :: xaux1,xaux2,yaux1,yaux2
+  real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
+  !HSO  end
+  integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
+  real :: sizesouth,sizenorth,xauxa,pint,conversion_factor
+
+  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
+
+  ! dimension of isec2 at least (22+n), where n is the number of parallels or
+  ! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
+
+  ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
+  ! coordinate parameters
+
+  integer :: isec1(56),isec2(22+nxmax+nymax)
+  real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+  character(len=1) :: opt
+
+  !HSO  grib api error messages
+  character(len=24) :: gribErrorMsg = 'Error reading grib file'
+  character(len=20) :: gribFunction = 'gridcheck'
+
+
+  iumax=0
+  iwmax=0
+
+  if(ideltas.gt.0) then
+    ifn=1
+  else
+    ifn=numbwf
+  endif
+  !
+  ! OPENING OF DATA FILE (GRIB CODE)
+  !
+5 call grib_open_file(ifile,path(3)(1:length(3)) &
+       //trim(wfname(ifn)),'r',iret)
+  if (iret.ne.GRIB_SUCCESS) then
+    goto 999   ! ERROR DETECTED
+  endif
+  !turn on support for multi fields messages
+  !call grib_multi_support_on
+
+  gotGrid=0
+  ifield=0
+10 ifield=ifield+1
+
+  !
+  ! GET NEXT FIELDS
+  !
+  call grib_new_from_file(ifile,igrib,iret)
+  if (iret.eq.GRIB_END_OF_FILE )  then
+    goto 30    ! EOF DETECTED
+  elseif (iret.ne.GRIB_SUCCESS) then
+    goto 999   ! ERROR DETECTED
+  endif
+
+  !first see if we read GRIB1 or GRIB2
+  call grib_get_int(igrib,'editionNumber',gribVer,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+
+  if (gribVer.eq.1) then ! GRIB Edition 1
+
+    !print*,'GRiB Edition 1'
+    !read the grib2 identifiers
+    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',isec1(8),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !change code for etadot to code for omega
+    if (isec1(6).eq.77) then
+      isec1(6)=135
+    endif
+
+    !print*,isec1(6),isec1(8)
+
+  else
+
+    !print*,'GRiB Edition 2'
+    !read the grib2 identifiers
+    call grib_get_int(igrib,'discipline',discipl,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterCategory',parCat,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterNumber',parNum,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',valSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'paramId',parId,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !print*,discipl,parCat,parNum,typSurf,valSurf
+
+    !convert to grib1 identifiers
+    isec1(6)=-1
+    isec1(7)=-1
+    isec1(8)=-1
+    isec1(8)=valSurf     ! level
+    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
+      isec1(6)=130         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
+      isec1(6)=131         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
+      isec1(6)=132         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
+      isec1(6)=133         ! indicatorOfParameter
+      !ZHG FOR CLOUDS FROM GRIB
+    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
+      isec1(6)=246         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
+      isec1(6)=247         ! indicatorOfParameter
+      !ZHG end
+      ! ESO qc(=clwc+ciwc)
+    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
+      isec1(6)=201031      ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
+      isec1(6)=134         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
+      isec1(6)=151         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
+      isec1(6)=165         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
+      isec1(6)=166         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
+      isec1(6)=167         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
+      isec1(6)=168         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
+      isec1(6)=141         ! indicatorOfParameter
+    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
+      isec1(6)=164         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
+      isec1(6)=142         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
+      isec1(6)=143         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
+      isec1(6)=146         ! indicatorOfParameter
+    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
+      isec1(6)=176         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
+      isec1(6)=180         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
+      isec1(6)=181         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
+      isec1(6)=129         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
+      isec1(6)=160         ! indicatorOfParameter
+    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
+         (typSurf.eq.1)) then ! LSM
+      isec1(6)=172         ! indicatorOfParameter
+    else
+      print*,'***ERROR: undefined GRiB2 message found!',discipl, &
+           parCat,parNum,typSurf
+    endif
+    if(parId .ne. isec1(6) .and. parId .ne. 77) then
+      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
+      !    stop
+    endif
+
+  endif
+
+  CALL grib_get_int(igrib,'numberOfPointsAlongAParallel', &
+       isec2(2),iret)
+  ! nx=isec2(2)
+  ! WRITE(*,*) nx,nxmax
+  IF (isec2(2).GT.nxmax) THEN
+    WRITE(*,*) 'FLEXPART error: Too many grid points in x direction.'
+    WRITE(*,*) 'Reduce resolution of wind fields.'
+    WRITE(*,*) 'Or change parameter settings in file ecmwf_mod.'
+    WRITE(*,*) nx,nxmax
+!    STOP
+  ENDIF
+
+  !get the size and data of the values array
+  if (isec1(6).ne.-1) then
+    call grib_get_real4_array(igrib,'values',zsec4,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+  endif
+
+  if (ifield.eq.1) then
+
+    !HSO  get the required fields from section 2 in a gribex compatible manner
+    call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
+         isec2(2),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
+         isec2(3),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
+         xaux1in,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
+         isec2(12),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    !  get the size and data of the vertical coordinate array
+    call grib_get_real4_array(igrib,'pv',zsec2,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+    nxfield=isec2(2)
+    ny=isec2(3)
+    nlev_ec=isec2(12)/2-1
+  endif
+
+  !HSO  get the second part of the grid dimensions only from GRiB1 messages
+  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
+    call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
+         xaux2in,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
+         yaux1in,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
+         yaux2in,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    xaux1=xaux1in
+    xaux2=xaux2in
+    yaux1=yaux1in
+    yaux2=yaux2in
+    if (xaux1.gt.180.) xaux1=xaux1-360.0
+    if (xaux2.gt.180.) xaux2=xaux2-360.0
+    if (xaux1.lt.-180.) xaux1=xaux1+360.0
+    if (xaux2.lt.-180.) xaux2=xaux2+360.0
+    if (xaux2.lt.xaux1) xaux2=xaux2+360.0
+    xlon0=xaux1
+    ylat0=yaux1
+    dx=(xaux2-xaux1)/real(nxfield-1)
+    dy=(yaux2-yaux1)/real(ny-1)
+    dxconst=180./(dx*r_earth*pi)
+    dyconst=180./(dy*r_earth*pi)
+    gotGrid=1
+    ! Check whether fields are global
+    ! If they contain the poles, specify polar stereographic map
+    ! projections using the stlmbr- and stcm2p-calls
+    !***********************************************************
+
+    xauxa=abs(xaux2+dx-360.-xaux1)
+    if (xauxa.lt.0.001) then
+      nx=nxfield+1                 ! field is cyclic
+      xglobal=.true.
+      if (abs(nxshift).ge.nx) &
+           stop 'nxshift in file par_mod is too large'
+      xlon0=xlon0+real(nxshift)*dx
+    else
+      nx=nxfield
+      xglobal=.false.
+      if (nxshift.ne.0) &
+           stop 'nxshift (par_mod) must be zero for non-global domain'
+    endif
+    nxmin1=nx-1
+    nymin1=ny-1
+    if (xlon0.gt.180.) xlon0=xlon0-360.
+    xauxa=abs(yaux1+90.)
+    if (xglobal.and.xauxa.lt.0.001) then
+      sglobal=.true.               ! field contains south pole
+      ! Enhance the map scale by factor 3 (*2=6) compared to north-south
+      ! map scale
+      sizesouth=6.*(switchsouth+90.)/dy
+      call stlmbr(southpolemap,-90.,0.)
+      call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
+           sizesouth,switchsouth,180.)
+      switchsouthg=(switchsouth-ylat0)/dy
+    else
+      sglobal=.false.
+      switchsouthg=999999.
+    endif
+    xauxa=abs(yaux2-90.)
+    if (xglobal.and.xauxa.lt.0.001) then
+      nglobal=.true.               ! field contains north pole
+      ! Enhance the map scale by factor 3 (*2=6) compared to north-south
+      ! map scale
+      sizenorth=6.*(90.-switchnorth)/dy
+      call stlmbr(northpolemap,90.,0.)
+      call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
+           sizenorth,switchnorth,180.)
+      switchnorthg=(switchnorth-ylat0)/dy
+    else
+      nglobal=.false.
+      switchnorthg=999999.
+    endif
+    if (nxshift.lt.0) &
+         stop 'nxshift (par_mod) must not be negative'
+    if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
+  endif ! gotGrid
+
+  if (nx.gt.nxmax) then
+    write(*,*) 'FLEXPART error: Too many grid points in x direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) nx,nxmax
+    stop
+  endif
+
+  if (ny.gt.nymax) then
+    write(*,*) 'FLEXPART error: Too many grid points in y direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) ny,nymax
+    stop
+  endif
+
+  k=isec1(8)
+  if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
+  if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
+
+  if(isec1(6).eq.129) then
+    do jy=0,ny-1
+      do ix=0,nxfield-1
+        oro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)/ga
+      end do
+    end do
+  endif
+  if(isec1(6).eq.172) then
+    do jy=0,ny-1
+      do ix=0,nxfield-1
+        lsm(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
+      end do
+    end do
+  endif
+  if(isec1(6).eq.160) then
+    do jy=0,ny-1
+      do ix=0,nxfield-1
+        excessoro(ix,jy)=zsec4(nxfield*(ny-jy-1)+ix+1)
+      end do
+    end do
+  endif
+
+  call grib_release(igrib)
+  goto 10                      !! READ NEXT LEVEL OR PARAMETER
+  !
+  ! CLOSING OF INPUT DATA FILE
+  !
+
+30 call grib_close_file(ifile)
+
+  !error message if no fields found with correct first longitude in it
+  if (gotGrid.eq.0) then
+    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
+         'messages'
+    stop
+  endif
+
+  nuvz=iumax
+  nwz =iwmax
+  if(nuvz.eq.nlev_ec) nwz=nlev_ec+1
+
+  if (nuvz+1.gt.nuvzmax) then
+    write(*,*) 'FLEXPART error: Too many u,v grid points in z '// &
+         'direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) nuvz+1,nuvzmax
+    stop
+  endif
+
+  if (nwz.gt.nwzmax) then
+    write(*,*) 'FLEXPART error: Too many w grid points in z '// &
+         'direction.'
+    write(*,*) 'Reduce resolution of wind fields.'
+    write(*,*) 'Or change parameter settings in file par_mod.'
+    write(*,*) nwz,nwzmax
+    stop
+  endif
+
+  ! If desired, shift all grids by nxshift grid cells
+  !**************************************************
+
+  if (xglobal) then
+    call shift_field_0(oro,nxfield,ny)
+    call shift_field_0(lsm,nxfield,ny)
+    call shift_field_0(excessoro,nxfield,ny)
+  endif
+
+  ! Output of grid info
+  !********************
+
+  if (lroot) then
+    write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
+         nuvz+1,nwz
+    write(*,*)
+    write(*,'(a)') ' Mother domain:'
+    write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
+         xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
+    write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
+         ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
+    write(*,*)
+  end if
+
+  ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
+  ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
+
+  numskip=nlev_ec-nuvz  ! number of ecmwf model layers not used
+  ! by trajectory model
+  !do 8940 i=1,244
+  !   write (*,*) 'zsec2:',i,ifield,zsec2(i),numskip
+  !940  continue
+  !   stop
+  ! SEC SEC SEC
+  ! for unknown reason zsec 1 to 10 is filled in this version
+  ! compared to the old one
+  ! SEC SEC SE
+  do i=1,nwz
+    j=numskip+i
+    k=nlev_ec+1+numskip+i
+    akm(nwz-i+1)=zsec2(j)
+    !   write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
+    bkm(nwz-i+1)=zsec2(k)
+  end do
+
+  !
+  ! CALCULATION OF AKZ, BKZ
+  ! AKZ,BKZ: model discretization parameters at the center of each model
+  !     layer
+  !
+  ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
+  ! i.e. ground level
+  !*****************************************************************************
+
+  akz(1)=0.
+  bkz(1)=1.0
+  do i=1,nuvz
+    akz(i+1)=0.5*(akm(i+1)+akm(i))
+    bkz(i+1)=0.5*(bkm(i+1)+bkm(i))
+  end do
+  nuvz=nuvz+1
+
+  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+  ! upon the transformation to z levels. In order to save computer memory, this is
+  ! not done anymore in the standard version. However, this option can still be
+  ! switched on by replacing the following lines with those below, that are
+  ! currently commented out. For this, similar changes are necessary in
+  ! verttransform.f and verttranform_nests.f
+  !*****************************************************************************
+
+  nz=nuvz
+  if (nz.gt.nzmax) stop 'nzmax too small'
+  do i=1,nuvz
+    aknew(i)=akz(i)
+    bknew(i)=bkz(i)
+  end do
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !nz=nuvz+nwz-1
+  !if (nz.gt.nzmax) stop 'nzmax too small'
+  !do 100 i=1,nwz
+  !  aknew(2*(i-1)+1)=akm(i)
+  !00     bknew(2*(i-1)+1)=bkm(i)
+  !do 110 i=2,nuvz
+  !  aknew(2*(i-1))=akz(i)
+  !10     bknew(2*(i-1))=bkz(i)
+  ! End doubled vertical resolution
+
+
+  ! Determine the uppermost level for which the convection scheme shall be applied
+  ! by assuming that there is no convection above 50 hPa (for standard SLP)
+  !*****************************************************************************
+
+  do i=1,nuvz-2
+    pint=akz(i)+bkz(i)*101325.
+    if (pint.lt.5000.) goto 96
+  end do
+96 nconvlev=i
+  if (nconvlev.gt.nconvlevmax-1) then
+    nconvlev=nconvlevmax-1
+    write(*,*) 'Attention, convection only calculated up to ', &
+         akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
+  endif
+
+  return
+
+999 write(*,*)
+  write(*,*) ' ###########################################'// &
+       '###### '
+  write(*,*) '       TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
+  write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
+  write(*,*) ' ###########################################'// &
+       '###### '
+  write(*,*)
+  write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND   !!!'
+  write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE...     !!!'
+  write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
+  write(*,'(a)') '!!! THE "X" KEY...                   !!!'
+  write(*,*)
+  read(*,'(a)') opt
+  if(opt.eq.'X') then
+    stop
+  else
+    goto 5
+  endif
+
+end subroutine gridcheck_ecmwf
+
Index: src/makefile
===================================================================
--- src/makefile	(revision d1a8707797c5c3fabdbd530210b84b85fb18b258)
+++ src/makefile	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -12,29 +12,38 @@
 #       'make -j ecmwf gcc=4.9',
 #    also set environment variable LD_LIBRARY_PATH to point to compiler libraries
+#
+#    Makefile was modified to produce unified executable for both ECMWF and GFS meteo data formats
+#    gributils were included to detect format of meteo data
+#
+#    Cpp directives USE_MPIINPLACE were added to three source files. The effect of these directives 
+#    are to enable the MPI_IN_PLACE option only if compiled with a -DUSE_MPIINPLACE directive. 
+#    Otherwise, a safer option (which requires the allocation of another array) is used by default. 
+#    In makefile added the -x f95-cpp-input flag for compiling of cpp directives.
 # 
 #  USAGE
-#    Compile serial FLEXPART (ECMWF)
-#      make [-j] ecmwf
-#
-#    Compile parallel FLEXPART (ECMWF)
-#      make [-j] ecmwf-mpi
+#    Compile serial FLEXPART 
+#      make [-j] serial
+#
+#    Compile parallel FLEXPART 
+#      make [-j] mpi
 #     
-#    Compile for debugging parallel FLEXPART (ECMWF)
-#      make [-j] ecmwf-mpi-dbg
-#
-#    Compile serial FLEXPART (GFS)
-#      make [-j] gfs
-#
-#    Compile parallel FLEXPART (GFS)
-#      make [-j] gfs-mpi
+#    Compile for debugging parallel FLEXPART
+#      make [-j] mpi-dbg
 #
 ################################################################################
 
 ## PROGRAMS
-FLEXPART-ECMWF-MPI 	= FP_ecmwf_MPI
-FLEXPART-ECMWF-MPI-DBG 	= DBG_FP_ecmwf_MPI
-FLEXPART-ECMWF 		= FP_ecmwf_gfortran
-FLEXPART-GFS 		= FP_gfs_gfortran
-FLEXPART-GFS-MPI 	= FP_gfs_MPI
+# Unified executable names
+# The same executable is used for both ECMWF and GFS metdata
+
+# Parallel processing executable
+FLEXPART-MPI = FLEXPART_MPI
+
+# Parallel processing executable with debugging info
+FLEXPART-MPI-DBG = DBG_FLEXPART_MPI
+
+# Serial processing executable
+FLEXPART-SERIAL = FLEXPART
+
 
 ifeq ($(gcc), 4.9)
@@ -59,4 +68,8 @@
 	LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
 endif
+
+
+# path to gributils used to detect meteodata format
+VPATH = gributils/
 
 
@@ -82,5 +95,6 @@
 xmass_mod.o             flux_mod.o \
 point_mod.o             outg_mod.o \
-mean_mod.o	 	random_mod.o
+mean_mod.o	 	random_mod.o \
+class_gribfile_mod.o
 
 MPI_MODOBJS = \
@@ -99,5 +113,6 @@
 	redist.o		\
 	concoutput_surf.o	concoutput_surf_nest.o	\
-	getfields.o
+	getfields.o \
+        readwind_ecmwf.o
 
 ## For MPI version 
@@ -112,32 +127,6 @@
 	redist_mpi.o		\
 	concoutput_surf_mpi.o	concoutput_surf_nest_mpi.o	\
-	getfields_mpi.o
-
-### WINDFIELDS
-## For ECMWF (serial) version:
-OBJECTS_ECMWF = \
-	calcpar.o          readwind.o \
-	richardson.o       verttransform.o \
-	obukhov.o          gridcheck.o  \
-	convmix.o          calcmatrix.o \
-	ecmwf_mod.o
-
-
-## For ECMWF MPI version:
-OBJECTS_ECMWF_MPI = \
-	gridcheck.o	   readwind_mpi.o \
-	calcpar.o          \
-	richardson.o       verttransform.o \
-	obukhov.o          \
-	convmix.o	   calcmatrix.o \
-	ecmwf_mod.o
-
-## For GFS (serial) version:
-OBJECTS_GFS = \
-	calcpar_gfs.o          readwind_gfs.o \
-	richardson_gfs.o       verttransform_gfs.o \
-	obukhov_gfs.o          gridcheck_gfs.o  \
-	convmix_gfs.o          calcmatrix_gfs.o \
-	gfs_mod.o
+	getfields_mpi.o \
+        readwind_ecmwf_mpi.o 
 
 OBJECTS = \
@@ -201,62 +190,28 @@
 %.o: %.mod
 
-ecmwf: $(FLEXPART-ECMWF)
-ecmwf: FC := $(F90)
-
-ecmwf-mpi: $(FLEXPART-ECMWF-MPI)
-ecmwf-mpi: FC := $(MPIF90)
-
-ecmwf-mpi-dbg: $(FLEXPART-ECMWF-MPI-DBG)
-ecmwf-mpi-dbg: FFLAGS := $(DBGFLAGS) 
-ecmwf-mpi-dbg: LDFLAGS:= $(LDDEBUG) 
-ecmwf-mpi-dbg: FC := $(MPIF90)
-
-gfs: $(FLEXPART-GFS)
-gfs: FC := $(F90)
-
-gfs-mpi: $(FLEXPART-GFS-MPI)
-gfs-mpi: FC := $(MPIF90)
-
-#all: $(FLEXPART-ECMWF)
-#all: $(FLEXPART-ECMWF-MPI)
-
-## This allows for switching between ECMWF/GFS without editing source code.
-wind_mod = ecmwf_mod.o # default wind
-# ifeq ($(MAKECMDGOALS),ecmwf)
-# wind_mod = ecmwf_mod.o
-# endif
-# ifeq ($(MAKECMDGOALS),ecmwf-mpi)
-# wind_mod = ecmwf_mod.o
-# endif
-# ifeq ($(MAKECMDGOALS),ecmwf-mpi-dbg)
-# wind_mod = ecmwf_mod.o
-# endif
-
-ifeq ($(MAKECMDGOALS),gfs)
-wind_mod = gfs_mod.o
-endif
-ifeq ($(MAKECMDGOALS),gfs-mpi)
-wind_mod = gfs_mod.o
-endif
-
-$(FLEXPART-ECMWF): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_ECMWF)
-	+$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_ECMWF) $(LDFLAGS)
-
-$(FLEXPART-ECMWF-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) $(OBJECTS_ECMWF_MPI)
+# serial executable
+serial: $(FLEXPART-SERIAL)
+serial: FC := $(F90)
+
+# parallel processing executable
+mpi: $(FLEXPART-MPI)
+mpi: FC := $(MPIF90)
+
+# parallel processing with debugging info
+mpi-dbg: $(FLEXPART-MPI-DBG)
+mpi-dbg: FFLAGS := $(DBGFLAGS)
+mpi-dbg: LDFLAGS:= $(LDDEBUG)
+mpi-dbg: FC := $(MPIF90)
+
+$(FLEXPART-SERIAL): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL)
+	+$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(LDFLAGS)
+
+$(FLEXPART-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI)
 	+$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \
-	$(OBJECTS_ECMWF_MPI) $(LDFLAGS)
-#	+$(FC) -o $@ *.o $(LDFLAGS)
-
-$(FLEXPART-ECMWF-MPI-DBG): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \
-	$(OBJECTS_ECMWF_MPI)
+	$(LDFLAGS)
+
+$(FLEXPART-MPI-DBG): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI)
 	+$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \
-	$(OBJECTS_ECMWF_MPI) $(LDFLAGS)
-
-$(FLEXPART-GFS): $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_GFS)
-	+$(FC) -o $@ $(MODOBJS) $(OBJECTS) $(OBJECTS_SERIAL) $(OBJECTS_GFS) $(LDFLAGS)
-
-$(FLEXPART-GFS-MPI): $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) $(OBJECTS_GFS)
-	+$(FC) -o $@ $(MODOBJS) $(MPI_MODOBJS) $(OBJECTS) $(OBJECTS_MPI) \
-	$(OBJECTS_GFS) $(LDFLAGS)
+	$(LDFLAGS)
 
 %.o: %.f90
@@ -267,6 +222,6 @@
 
 cleanall:
-	\rm -f *.o *.mod $(FLEXPART-ECMWF-MPI) $(FLEXPART-ECMWF-MPI-DBG) $(FLEXPART-ECMWF) \
-	$(FLEXPART-GFS-MPI) $(FLEXPART-GFS)
+	\rm -f *.o *.mod $(FLEXPART-MPI) $(FLEXPART-MPI-DBG) $(FLEXPART-SERIAL)
+
 
 .SUFFIXES = $(SUFFIXES) .f90
@@ -282,8 +237,6 @@
 	random_mod.o
 calcfluxes.o: com_mod.o flux_mod.o outg_mod.o par_mod.o
-calcmatrix.o: com_mod.o conv_mod.o par_mod.o
-calcmatrix_gfs.o: com_mod.o conv_mod.o par_mod.o
-calcpar.o: com_mod.o par_mod.o
-calcpar_gfs.o: com_mod.o par_mod.o
+calcmatrix.o: com_mod.o conv_mod.o par_mod.o class_gribfile_mod.o
+calcpar.o: com_mod.o par_mod.o class_gribfile_mod.o
 calcpar_nests.o: com_mod.o par_mod.o
 calcpv.o: com_mod.o par_mod.o
@@ -311,7 +264,7 @@
 conv_mod.o: par_mod.o
 convect43c.o: conv_mod.o par_mod.o
-convmix.o: com_mod.o conv_mod.o flux_mod.o par_mod.o
-convmix_gfs.o: com_mod.o conv_mod.o par_mod.o
+convmix.o: com_mod.o conv_mod.o flux_mod.o par_mod.o class_gribfile_mod.o
 coordtrafo.o: com_mod.o par_mod.o point_mod.o
+detectformat.o: com_mod.o par_mod.o class_gribfile_mod.o
 distance.o: par_mod.o
 distance2.o: par_mod.o
@@ -319,11 +272,11 @@
 drydepokernel_nest.o: com_mod.o par_mod.o unc_mod.o
 erf.o: par_mod.o
-FLEXPART.o: com_mod.o conv_mod.o par_mod.o point_mod.o random_mod.o netcdf_output_mod.o
+FLEXPART.o: com_mod.o conv_mod.o par_mod.o point_mod.o random_mod.o netcdf_output_mod.o class_gribfile_mod.o
 FLEXPART_MPI.o: com_mod.o conv_mod.o mpi_mod.o par_mod.o point_mod.o \
-	random_mod.o netcdf_output_mod.o
+	random_mod.o netcdf_output_mod.o class_gribfile_mod.o
 fluxoutput.o: com_mod.o flux_mod.o outg_mod.o par_mod.o
 get_settling.o: com_mod.o par_mod.o
-getfields.o: com_mod.o par_mod.o
-getfields_mpi.o: com_mod.o par_mod.o mpi_mod.o
+getfields.o: com_mod.o par_mod.o class_gribfile_mod.o
+getfields_mpi.o: com_mod.o par_mod.o mpi_mod.o class_gribfile_mod.o
 gethourlyOH.o: com_mod.o oh_mod.o par_mod.o
 getrb.o: par_mod.o
@@ -331,5 +284,5 @@
 getvdep.o: com_mod.o par_mod.o
 getvdep_nests.o: com_mod.o par_mod.o
-gridcheck.o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o
+gridcheck_ecmwf.o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o
 gridcheck_emos.o: com_mod.o conv_mod.o par_mod.o
 gridcheck_fnl.o: cmapf_mod.o com_mod.o conv_mod.o par_mod.o
@@ -374,5 +327,5 @@
 outgrid_init.o: com_mod.o flux_mod.o oh_mod.o outg_mod.o par_mod.o unc_mod.o
 outgrid_init_nest.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
-par_mod.o : $(wind_mod)
+par_mod.o :
 part0.o: par_mod.o
 partdep.o: par_mod.o
@@ -402,9 +355,9 @@
 readreleases.o: com_mod.o par_mod.o point_mod.o xmass_mod.o
 readspecies.o: com_mod.o par_mod.o
-readwind.o: com_mod.o par_mod.o
+readwind_ecmwf.o: com_mod.o par_mod.o
 readwind_emos.o: com_mod.o par_mod.o
 readwind_gfs.o: com_mod.o par_mod.o
 readwind_gfs_emos.o: com_mod.o par_mod.o
-readwind_mpi.o: com_mod.o mpi_mod.o par_mod.o
+readwind_ecmwf_mpi.o: com_mod.o mpi_mod.o par_mod.o
 readwind_nests.o: com_mod.o par_mod.o
 readwind_nests_emos.o: com_mod.o par_mod.o
@@ -415,6 +368,5 @@
 releaseparticles_mpi.o: com_mod.o mpi_mod.o par_mod.o point_mod.o \
 	random_mod.o xmass_mod.o
-richardson.o: par_mod.o
-richardson_gfs.o: par_mod.o
+richardson.o: par_mod.o class_gribfile_mod.o
 scalev.o: par_mod.o
 shift_field.o: par_mod.o
@@ -425,5 +377,5 @@
 	par_mod.o point_mod.o unc_mod.o xmass_mod.o netcdf_output_mod.o	
 unc_mod.o: par_mod.o
-verttransform.o: cmapf_mod.o com_mod.o par_mod.o 
+verttransform_ecmwf.o: cmapf_mod.o com_mod.o par_mod.o 
 verttransform_gfs.o: cmapf_mod.o com_mod.o par_mod.o
 verttransform_nests.o: com_mod.o par_mod.o
Index: src/par_mod.f90
===================================================================
--- src/par_mod.f90	(revision d1a8707797c5c3fabdbd530210b84b85fb18b258)
+++ src/par_mod.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -37,9 +37,4 @@
 
 module par_mod
-
-!************************************************************************
-! wind_mod: is gfs_mod.f90 for target gfs, ecmwf_mod.f90 for target ecmwf
-!************************************************************************
-  use wind_mod 
 
   implicit none
@@ -146,5 +141,32 @@
   !*********************************************
   
-  ! Moved to ecmwf_mod.f90 (for ECMWF) / gfs_mod.f90 (for GFS)
+!  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF new 
+  integer,parameter :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !ECMWF new 
+
+!  INTEGER,PARAMETER :: nxmax=361,nymax=181,nuvzmax=138,nwzmax=138,nzmax=138 !NCEP data
+
+!  !integer,parameter :: nxshift=359  ! for ECMWF
+!  integer,parameter :: nxshift=0     ! for GFS
+!  integer,parameter :: nxmax=15,nymax=15,nuvzmax=140,nwzmax=140,nzmax=140
+  integer,parameter :: nxshift=359  ! for ECMWF
+!  integer,parameter :: nxshift=0     ! for GFS
+
+  !*********************************************
+  ! Maximum dimensions of the nested input grids
+  !*********************************************
+
+  integer,parameter :: maxnests=0,nxmaxn=451,nymaxn=226
+
+  ! nxmax,nymax        maximum dimension of wind fields in x and y
+  !                    direction, respectively
+  ! nuvzmax,nwzmax     maximum dimension of (u,v) and (w) wind fields in z
+  !                    direction (for fields on eta levels)
+  ! nzmax              maximum dimension of wind fields in z direction
+  !                    for the transformed Cartesian coordinates
+  ! nxshift            for global grids (in x), the grid can be shifted by
+  !                    nxshift grid points, in order to accomodate nested
+  !                    grids, and output grids overlapping the domain "boundary"
+  !                    nxshift must not be negative; "normal" setting would be 0
+
   
   integer,parameter :: nconvlevmax = nuvzmax-1
Index: src/readwind_ecmwf.f90
===================================================================
--- src/readwind_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/readwind_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,508 @@
+!**********************************************************************
+! 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/>.   *
+!**********************************************************************
+
+subroutine readwind_ecmwf(indj,n,uuh,vvh,wwh)
+
+!**********************************************************************
+!                                                                     *
+!             TRAJECTORY MODEL SUBROUTINE READWIND                    *
+!                                                                     *
+!**********************************************************************
+!                                                                     *
+!             AUTHOR:      G. WOTAWA                                  *
+!             DATE:        1997-08-05                                 *
+!             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
+!             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
+!                                 ECMWF grib_api                      *
+!             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
+!                                 ECMWF grib_api                      *
+!                                                                     *
+!**********************************************************************
+!  Changes, Bernd C. Krueger, Feb. 2001:
+!   Variables tth and qvh (on eta coordinates) in common block
+!
+!   Unified ECMWF and GFS builds                                      
+!   Marian Harustak, 12.5.2017                                        
+!     - Renamed from readwind to readwind_ecmwf                     
+!
+!**********************************************************************
+!                                                                     *
+! DESCRIPTION:                                                        *
+!                                                                     *
+! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
+! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
+!                                                                     *
+! INPUT:                                                              *
+! indj               indicates number of the wind field to be read in *
+! n                  temporal index for meteorological fields (1 to 3)*
+!                                                                     *
+! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
+!                                                                     *
+! wfname             File name of data to be read in                  *
+! nx,ny,nuvz,nwz     expected field dimensions                        *
+! nlev_ec            number of vertical levels ecmwf model            *
+! uu,vv,ww           wind fields                                      *
+! tt,qv              temperature and specific humidity                *
+! ps                 surface pressure                                 *
+!                                                                     *
+!**********************************************************************
+
+  use grib_api
+  use par_mod
+  use com_mod
+
+  implicit none
+
+!  include 'grib_api.h'
+
+!HSO  parameters for grib_api
+  integer :: ifile
+  integer :: iret
+  integer :: igrib
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
+  integer :: gotGrid
+!HSO  end
+
+  real(sp) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real(sp) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real(sp) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+  integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
+
+! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
+
+! dimension of isec2 at least (22+n), where n is the number of parallels or
+! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
+
+! dimension of zsec2 at least (10+nn), where nn is the number of vertical
+! coordinate parameters
+
+  integer :: isec1(56),isec2(22+nxmax+nymax)
+  real(sp) :: zsec4(jpunp)
+  real(sp) :: xaux,yaux
+  real(dp) :: xauxin,yauxin
+  real(sp),parameter :: eps=1.e-4
+  real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
+  real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
+
+  logical :: hflswitch,strswitch!,readcloud
+
+!HSO  grib api error messages
+  character(len=24) :: gribErrorMsg = 'Error reading grib file'
+  character(len=20) :: gribFunction = 'readwind'
+
+  hflswitch=.false.
+  strswitch=.false.
+!ZHG test the grib fields that have lcwc without using them
+!  readcloud=.false.
+
+  levdiff2=nlev_ec-nwz+1
+  iumax=0
+  iwmax=0
+
+!
+! OPENING OF DATA FILE (GRIB CODE)
+!
+5 call grib_open_file(ifile,path(3)(1:length(3)) &
+       //trim(wfname(indj)),'r',iret)
+  if (iret.ne.GRIB_SUCCESS) then
+    goto 888   ! ERROR DETECTED
+  endif
+!turn on support for multi fields messages */
+!call grib_multi_support_on
+
+  gotGrid=0
+  ifield=0
+10 ifield=ifield+1
+!
+! GET NEXT FIELDS
+!
+  call grib_new_from_file(ifile,igrib,iret)
+  if (iret.eq.GRIB_END_OF_FILE)  then
+    goto 50    ! EOF DETECTED
+  elseif (iret.ne.GRIB_SUCCESS) then
+    goto 888   ! ERROR DETECTED
+  endif
+
+!first see if we read GRIB1 or GRIB2
+  call grib_get_int(igrib,'editionNumber',gribVer,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+
+  if (gribVer.eq.1) then ! GRIB Edition 1
+
+!print*,'GRiB Edition 1'
+!read the grib2 identifiers
+    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',isec1(8),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+!change code for etadot to code for omega
+    if (isec1(6).eq.77) then
+      isec1(6)=135
+    endif
+
+    conversion_factor=1.
+
+  else
+
+!print*,'GRiB Edition 2'
+!read the grib2 identifiers
+    call grib_get_int(igrib,'discipline',discipl,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterCategory',parCat,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterNumber',parNum,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',valSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'paramId',parId,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+!print*,discipl,parCat,parNum,typSurf,valSurf
+
+!convert to grib1 identifiers
+    isec1(6)=-1
+    isec1(7)=-1
+    isec1(8)=-1
+    isec1(8)=valSurf     ! level
+    conversion_factor=1.
+    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
+      isec1(6)=130         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
+      isec1(6)=131         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
+      isec1(6)=132         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
+      isec1(6)=133         ! indicatorOfParameter
+! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field QC 
+    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
+      isec1(6)=246         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
+      isec1(6)=247         ! indicatorOfParameter
+! ESO qc(=clwc+ciwc):
+    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
+      isec1(6)=201031         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
+      isec1(6)=134         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
+      isec1(6)=151         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
+      isec1(6)=165         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
+      isec1(6)=166         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
+      isec1(6)=167         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
+      isec1(6)=168         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
+      isec1(6)=141         ! indicatorOfParameter
+      conversion_factor=1000.
+    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
+      isec1(6)=164         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
+      isec1(6)=142         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
+      isec1(6)=143         ! indicatorOfParameter
+      conversion_factor=1000.
+    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
+      isec1(6)=146         ! indicatorOfParameter
+    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
+      isec1(6)=176         ! indicatorOfParameter
+!    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS --wrong
+    elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
+      isec1(6)=180         ! indicatorOfParameter
+!    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
+    elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
+      isec1(6)=181         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
+      isec1(6)=129         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
+      isec1(6)=160         ! indicatorOfParameter
+    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
+         (typSurf.eq.1)) then ! LSM
+      isec1(6)=172         ! indicatorOfParameter
+    else
+      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
+           parCat,parNum,typSurf
+    endif
+    if(parId .ne. isec1(6) .and. parId .ne. 77) then
+      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
+!    stop
+    endif
+
+  endif
+
+!HSO  get the size and data of the values array
+  if (isec1(6).ne.-1) then
+    call grib_get_real4_array(igrib,'values',zsec4,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+  endif
+
+!HSO  get the required fields from section 2 in a gribex compatible manner
+  if (ifield.eq.1) then
+    call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
+    call grib_check(iret,gribFunction,gribErrorMsg)
+! CHECK GRID SPECIFICATIONS
+    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
+    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
+    if(isec2(12)/2-1.ne.nlev_ec) &
+         stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
+  endif ! ifield
+
+!HSO  get the second part of the grid dimensions only from GRiB1 messages
+  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
+    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
+         xauxin,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
+         yauxin,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    if (xauxin.gt.180.) xauxin=xauxin-360.0
+    if (xauxin.lt.-180.) xauxin=xauxin+360.0
+
+    xaux=xauxin+real(nxshift)*dx
+    yaux=yauxin
+    if (xaux.gt.180.) xaux=xaux-360.0
+    if(abs(xaux-xlon0).gt.eps) &
+         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
+    if(abs(yaux-ylat0).gt.eps) &
+         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
+    gotGrid=1
+  endif ! gotGrid
+
+  do j=0,nymin1
+    do i=0,nxfield-1
+      k=isec1(8)
+      if(isec1(6).eq.130) tth(i,j,nlev_ec-k+2,n)= &!! TEMPERATURE
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.131) uuh(i,j,nlev_ec-k+2)= &!! U VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.132) vvh(i,j,nlev_ec-k+2)= &!! V VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.133) then                      !! SPEC. HUMIDITY
+        qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
+             qvh(i,j,nlev_ec-k+2,n) = 0.
+!        this is necessary because the gridded data may contain
+!        spurious negative values
+      endif
+      if(isec1(6).eq.134) ps(i,j,1,n)= &!! SURF. PRESS.
+           zsec4(nxfield*(ny-j-1)+i+1)
+
+      if(isec1(6).eq.135) wwh(i,j,nlev_ec-k+1)= &!! W VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
+           zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
+      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.164) tcc(i,j,1,n)= &!! CLOUD COVER
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.165) u10(i,j,1,n)= &!! 10 M U VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.166) v10(i,j,1,n)= &!! 10 M V VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.167) tt2(i,j,1,n)= &!! 2 M TEMPERATURE
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.168) td2(i,j,1,n)= &!! 2 M DEW POINT
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.142) then                      !! LARGE SCALE PREC.
+        lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
+      endif
+      if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
+        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
+        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
+      endif
+      if(isec1(6).eq.146) sshf(i,j,1,n)= &!! SENS. HEAT FLUX
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if((isec1(6).eq.146).and.(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) &
+           hflswitch=.true.    ! Heat flux available
+      if(isec1(6).eq.176) then                      !! SOLAR RADIATION
+        ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
+      endif
+      if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
+           (zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
+!sec        strswitch=.true.
+      if(isec1(6).eq.129) oro(i,j)= &!! ECMWF OROGRAPHY
+           zsec4(nxfield*(ny-j-1)+i+1)/ga
+      if(isec1(6).eq.160) excessoro(i,j)= &!! STANDARD DEVIATION OF OROGRAPHY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.172) lsm(i,j)= &!! ECMWF LAND SEA MASK
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
+      if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
+!ZHG READING CLOUD FIELDS ASWELL
+! ESO TODO: add check for if one of clwc/ciwc missing (error),
+! also if all 3 cw fields present, use qc and disregard the others
+      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
+        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        readclouds=.true.
+        sumclouds=.false.
+      endif
+      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
+        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+      endif
+!ZHG end
+!ESO read qc (=clwc+ciwc)
+      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
+        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        readclouds=.true.
+        sumclouds=.true.
+      endif
+
+    end do
+  end do
+
+  call grib_release(igrib)
+  goto 10                      !! READ NEXT LEVEL OR PARAMETER
+!
+! CLOSING OF INPUT DATA FILE
+!
+
+50 call grib_close_file(ifile)
+
+!error message if no fields found with correct first longitude in it
+  if (gotGrid.eq.0) then
+    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
+         'messages'
+    stop
+  endif
+
+  if(levdiff2.eq.0) then
+    iwmax=nlev_ec+1
+    do i=0,nxmin1
+      do j=0,nymin1
+        wwh(i,j,nlev_ec+1)=0.
+      end do
+    end do
+  endif
+
+! For global fields, assign the leftmost data column also to the rightmost
+! data column; if required, shift whole grid by nxshift grid points
+!*************************************************************************
+
+  if (xglobal) then
+    call shift_field_0(ewss,nxfield,ny)
+    call shift_field_0(nsss,nxfield,ny)
+    call shift_field_0(oro,nxfield,ny)
+    call shift_field_0(excessoro,nxfield,ny)
+    call shift_field_0(lsm,nxfield,ny)
+    call shift_field(ps,nxfield,ny,1,1,2,n)
+    call shift_field(sd,nxfield,ny,1,1,2,n)
+    call shift_field(msl,nxfield,ny,1,1,2,n)
+    call shift_field(tcc,nxfield,ny,1,1,2,n)
+    call shift_field(u10,nxfield,ny,1,1,2,n)
+    call shift_field(v10,nxfield,ny,1,1,2,n)
+    call shift_field(tt2,nxfield,ny,1,1,2,n)
+    call shift_field(td2,nxfield,ny,1,1,2,n)
+    call shift_field(lsprec,nxfield,ny,1,1,2,n)
+    call shift_field(convprec,nxfield,ny,1,1,2,n)
+    call shift_field(sshf,nxfield,ny,1,1,2,n)
+    call shift_field(ssr,nxfield,ny,1,1,2,n)
+    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
+    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
+    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
+    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
+    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
+!ZHG
+    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
+    if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
+!ZHG end
+
+  endif
+
+  do i=0,nxmin1
+    do j=0,nymin1
+      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
+    end do
+  end do
+
+  if ((.not.hflswitch).or.(.not.strswitch)) then
+    write(*,*) 'WARNING: No flux data contained in GRIB file ', &
+         wfname(indj)
+
+! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
+! As ECMWF has increased the model resolution, such that now the first model
+! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
+! (3rd model level in FLEXPART) for the profile method
+!***************************************************************************
+
+    do i=0,nxmin1
+      do j=0,nymin1
+        plev1=akz(3)+bkz(3)*ps(i,j,1,n)
+        pmean=0.5*(ps(i,j,1,n)+plev1)
+        tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
+        fu=-r_air*tv/ga/pmean
+        hlev1=fu*(plev1-ps(i,j,1,n))   ! HEIGTH OF FIRST MODEL LAYER
+        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
+        fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
+        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
+             tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
+             surfstr(i,j,1,n),sshf(i,j,1,n))
+        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
+        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
+      end do
+    end do
+  endif
+
+
+! Assign 10 m wind to model level at eta=1.0 to have one additional model
+! level at the ground
+! Specific humidity is taken the same as at one level above
+! Temperature is taken as 2 m temperature
+!**************************************************************************
+
+  do i=0,nxmin1
+    do j=0,nymin1
+      uuh(i,j,1)=u10(i,j,1,n)
+      vvh(i,j,1)=v10(i,j,1,n)
+      qvh(i,j,1,n)=qvh(i,j,2,n)
+      tth(i,j,1,n)=tt2(i,j,1,n)
+    end do
+  end do
+
+  if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
+  if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
+
+  return
+
+888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
+  write(*,*) ' #### ',wfname(indj),'                    #### '
+  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
+  stop 'Execution terminated'
+
+end subroutine readwind_ecmwf
+
Index: src/readwind_ecmwf_mpi.f90
===================================================================
--- src/readwind_ecmwf_mpi.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/readwind_ecmwf_mpi.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,533 @@
+!**********************************************************************
+! 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/>.   *
+!**********************************************************************
+
+subroutine readwind_ecmwf(indj,n,uuh,vvh,wwh)
+
+!**********************************************************************
+!                                                                     *
+!             TRAJECTORY MODEL SUBROUTINE READWIND                    *
+!                                                                     *
+!**********************************************************************
+!                                                                     *
+!             AUTHOR:      G. WOTAWA                                  *
+!             DATE:        1997-08-05                                 *
+!             LAST UPDATE: 2000-10-17, Andreas Stohl                  *
+!             CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
+!                                 ECMWF grib_api                      *
+!             CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
+!                                 ECMWF grib_api                      *
+!                                                                     *
+!**********************************************************************
+!  Changes, Bernd C. Krueger, Feb. 2001:
+!   Variables tth and qvh (on eta coordinates) in common block
+!  eso 2014:
+!    This version for reading windfields in advance by a separate
+!    MPI process.
+!    TODO: can be merged with readwind.f90 if using numwfmem in
+!          shift_field
+!
+!   Unified ECMWF and GFS builds             
+!   Marian Harustak, 12.5.2017                
+!     - Renamed from reawind to readwind_ecmwf 
+!
+!**********************************************************************
+!                                                                     *
+! DESCRIPTION:                                                        *
+!                                                                     *
+! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE   *
+! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE          *
+!                                                                     *
+! INPUT:                                                              *
+! indj               indicates number of the wind field to be read in *
+! n                  temporal index for meteorological fields (1 to 3)*
+!                                                                     *
+! IMPORTANT VARIABLES FROM COMMON BLOCK:                              *
+!                                                                     *
+! wfname             File name of data to be read in                  *
+! nx,ny,nuvz,nwz     expected field dimensions                        *
+! nlev_ec            number of vertical levels ecmwf model            *
+! uu,vv,ww           wind fields                                      *
+! tt,qv              temperature and specific humidity                *
+! ps                 surface pressure                                 *
+!                                                                     *
+!**********************************************************************
+
+  use GRIB_API
+  use par_mod
+  use com_mod
+  use mpi_mod
+
+  implicit none
+
+!  include 'grib_api.h'
+
+!HSO  parameters for grib_api
+  integer :: ifile
+  integer :: iret
+  integer :: igrib
+  integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
+  integer :: gotGrid
+!HSO  end
+
+  real(sp), intent(inout) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real(sp), intent(inout) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real(sp), intent(inout) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+  integer, intent(in) :: indj,n
+  integer :: i,j,k,levdiff2,ifield,iumax,iwmax
+
+! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
+
+! dimension of isec2 at least (22+n), where n is the number of parallels or
+! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
+
+! dimension of zsec2 at least (10+nn), where nn is the number of vertical
+! coordinate parameters
+
+  integer :: isec1(56),isec2(22+nxmax+nymax)
+  real(kind=4) :: zsec4(jpunp)
+  real(kind=4) :: xaux,yaux
+  real(kind=8) :: xauxin,yauxin
+  real,parameter :: eps=1.e-4
+  real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
+  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
+
+  logical :: hflswitch,strswitch !,readcloud
+
+!HSO  grib api error messages
+  character(len=24) :: gribErrorMsg = 'Error reading grib file'
+  character(len=20) :: gribFunction = 'readwind'
+
+! measure time
+  if (mp_measure_time) call mpif_mtime('readwind',0)
+
+  hflswitch=.false.
+  strswitch=.false.
+!ZHG test the grib fields that have lcwc without using them
+!  readcloud=.false.
+!hg end
+  levdiff2=nlev_ec-nwz+1
+  iumax=0
+  iwmax=0
+
+!
+! OPENING OF DATA FILE (GRIB CODE)
+!
+5 call grib_open_file(ifile,path(3)(1:length(3)) &
+       //trim(wfname(indj)),'r',iret)
+  if (iret.ne.GRIB_SUCCESS) then
+    goto 888   ! ERROR DETECTED
+  endif
+!turn on support for multi fields messages */
+!call grib_multi_support_on
+
+  gotGrid=0
+  ifield=0
+10 ifield=ifield+1
+!
+! GET NEXT FIELDS
+!
+  call grib_new_from_file(ifile,igrib,iret)
+  if (iret.eq.GRIB_END_OF_FILE)  then
+    goto 50    ! EOF DETECTED
+  elseif (iret.ne.GRIB_SUCCESS) then
+    goto 888   ! ERROR DETECTED
+  endif
+
+!first see if we read GRIB1 or GRIB2
+  call grib_get_int(igrib,'editionNumber',gribVer,iret)
+  call grib_check(iret,gribFunction,gribErrorMsg)
+
+  if (gribVer.eq.1) then ! GRIB Edition 1
+
+!print*,'GRiB Edition 1'
+!read the grib2 identifiers
+    call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',isec1(8),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+!change code for etadot to code for omega
+    if (isec1(6).eq.77) then
+      isec1(6)=135
+    endif
+
+    conversion_factor=1.
+
+  else
+
+!print*,'GRiB Edition 2'
+!read the grib2 identifiers
+    call grib_get_int(igrib,'discipline',discipl,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterCategory',parCat,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'parameterNumber',parNum,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'level',valSurf,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'paramId',parId,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+
+!print*,discipl,parCat,parNum,typSurf,valSurf
+
+!convert to grib1 identifiers
+    isec1(6)=-1
+    isec1(7)=-1
+    isec1(8)=-1
+    isec1(8)=valSurf     ! level
+    conversion_factor=1.
+    if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
+      isec1(6)=130         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
+      isec1(6)=131         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
+      isec1(6)=132         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
+      isec1(6)=133         ! indicatorOfParameter
+! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field qc 
+    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
+      isec1(6)=246         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
+      isec1(6)=247         ! indicatorOfParameter
+! ESO qc(=clwc+ciwc):
+    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
+      isec1(6)=201031         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
+      isec1(6)=134         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
+      isec1(6)=135         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
+      isec1(6)=151         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
+      isec1(6)=165         ! indicatorOfParameter
+    elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
+      isec1(6)=166         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
+      isec1(6)=167         ! indicatorOfParameter
+    elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
+      isec1(6)=168         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
+      isec1(6)=141         ! indicatorOfParameter
+      conversion_factor=1000.
+    elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
+      isec1(6)=164         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
+      isec1(6)=142         ! indicatorOfParameter
+    elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
+      isec1(6)=143         ! indicatorOfParameter
+      conversion_factor=1000.
+    elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
+      isec1(6)=146         ! indicatorOfParameter
+    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
+      isec1(6)=176         ! indicatorOfParameter
+!    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS --wrong
+    elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
+      isec1(6)=180         ! indicatorOfParameter
+!    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
+    elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
+      isec1(6)=181         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
+      isec1(6)=129         ! indicatorOfParameter
+    elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
+      isec1(6)=160         ! indicatorOfParameter
+    elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
+         (typSurf.eq.1)) then ! LSM
+      isec1(6)=172         ! indicatorOfParameter
+    else
+      print*,'***WARNING: undefined GRiB2 message found!',discipl, &
+           parCat,parNum,typSurf
+    endif
+    if(parId .ne. isec1(6) .and. parId .ne. 77) then
+      write(*,*) 'parId',parId, 'isec1(6)',isec1(6)
+!    stop
+    endif
+
+  endif
+
+!HSO  get the size and data of the values array
+  if (isec1(6).ne.-1) then
+    call grib_get_real4_array(igrib,'values',zsec4,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+  endif
+
+!HSO  get the required fields from section 2 in a gribex compatible manner
+  if (ifield.eq.1) then
+    call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
+    call grib_check(iret,gribFunction,gribErrorMsg)
+! CHECK GRID SPECIFICATIONS
+    if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
+    if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
+    if(isec2(12)/2-1.ne.nlev_ec) &
+         stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
+  endif ! ifield
+
+!HSO  get the second part of the grid dimensions only from GRiB1 messages
+  if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
+    call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
+         xauxin,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
+         yauxin,iret)
+    call grib_check(iret,gribFunction,gribErrorMsg)
+    if (xauxin.gt.180.) xauxin=xauxin-360.0
+    if (xauxin.lt.-180.) xauxin=xauxin+360.0
+
+    xaux=xauxin+real(nxshift)*dx
+    yaux=yauxin
+    if (xaux.gt.180.) xaux=xaux-360.0
+    if(abs(xaux-xlon0).gt.eps) &
+         stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
+    if(abs(yaux-ylat0).gt.eps) &
+         stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
+    gotGrid=1
+  endif ! gotGrid
+
+  do j=0,nymin1
+    do i=0,nxfield-1
+      k=isec1(8)
+      if(isec1(6).eq.130) tth(i,j,nlev_ec-k+2,n)= &!! TEMPERATURE
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.131) uuh(i,j,nlev_ec-k+2)= &!! U VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.132) vvh(i,j,nlev_ec-k+2)= &!! V VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.133) then                      !! SPEC. HUMIDITY
+        qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
+             qvh(i,j,nlev_ec-k+2,n) = 0.
+!        this is necessary because the gridded data may contain
+!        spurious negative values
+      endif
+      if(isec1(6).eq.134) ps(i,j,1,n)= &!! SURF. PRESS.
+           zsec4(nxfield*(ny-j-1)+i+1)
+
+      if(isec1(6).eq.135) wwh(i,j,nlev_ec-k+1)= &!! W VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
+           zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
+      if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.164) tcc(i,j,1,n)= &!! CLOUD COVER
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.165) u10(i,j,1,n)= &!! 10 M U VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.166) v10(i,j,1,n)= &!! 10 M V VELOCITY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.167) tt2(i,j,1,n)= &!! 2 M TEMPERATURE
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.168) td2(i,j,1,n)= &!! 2 M DEW POINT
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.142) then                      !! LARGE SCALE PREC.
+        lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
+      endif
+      if(isec1(6).eq.143) then                      !! CONVECTIVE PREC.
+        convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
+        if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
+      endif
+      if(isec1(6).eq.146) sshf(i,j,1,n)= &!! SENS. HEAT FLUX
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if((isec1(6).eq.146).and.(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) &
+           hflswitch=.true.    ! Heat flux available
+      if(isec1(6).eq.176) then                      !! SOLAR RADIATION
+        ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
+      endif
+      if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
+           (zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
+!sec        strswitch=.true.
+      if(isec1(6).eq.129) oro(i,j)= &!! ECMWF OROGRAPHY
+           zsec4(nxfield*(ny-j-1)+i+1)/ga
+      if(isec1(6).eq.160) excessoro(i,j)= &!! STANDARD DEVIATION OF OROGRAPHY
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.172) lsm(i,j)= &!! ECMWF LAND SEA MASK
+           zsec4(nxfield*(ny-j-1)+i+1)
+      if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
+      if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
+!ZHG READING CLOUD FIELDS ASWELL
+      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
+        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        readclouds = .true.
+        sumclouds=.false.
+      endif
+      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
+        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+      endif
+!ZHG end
+!ESO read qc (=clwc+ciwc)
+      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
+        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
+        readclouds=.true.
+        sumclouds=.true.
+      endif
+
+    end do
+  end do
+
+  call grib_release(igrib)
+  goto 10                      !! READ NEXT LEVEL OR PARAMETER
+!
+! CLOSING OF INPUT DATA FILE
+!
+
+50 call grib_close_file(ifile)
+
+!error message if no fields found with correct first longitude in it
+  if (gotGrid.eq.0) then
+    print*,'***ERROR: input file needs to contain GRiB1 formatted'// &
+         'messages'
+    stop
+  endif
+
+  if(levdiff2.eq.0) then
+    iwmax=nlev_ec+1
+    do i=0,nxmin1
+      do j=0,nymin1
+        wwh(i,j,nlev_ec+1)=0.
+      end do
+    end do
+  endif
+
+! For global fields, assign the leftmost data column also to the rightmost
+! data column; if required, shift whole grid by nxshift grid points
+!*************************************************************************
+
+  if (xglobal) then
+    call shift_field_0(ewss,nxfield,ny)
+    call shift_field_0(nsss,nxfield,ny)
+    call shift_field_0(oro,nxfield,ny)
+    call shift_field_0(excessoro,nxfield,ny)
+    call shift_field_0(lsm,nxfield,ny)
+    ! call shift_field(ps,nxfield,ny,1,1,2,n)
+    ! call shift_field(sd,nxfield,ny,1,1,2,n)
+    ! call shift_field(msl,nxfield,ny,1,1,2,n)
+    ! call shift_field(tcc,nxfield,ny,1,1,2,n)
+    ! call shift_field(u10,nxfield,ny,1,1,2,n)
+    ! call shift_field(v10,nxfield,ny,1,1,2,n)
+    ! call shift_field(tt2,nxfield,ny,1,1,2,n)
+    ! call shift_field(td2,nxfield,ny,1,1,2,n)
+    ! call shift_field(lsprec,nxfield,ny,1,1,2,n)
+    ! call shift_field(convprec,nxfield,ny,1,1,2,n)
+    ! call shift_field(sshf,nxfield,ny,1,1,2,n)
+    ! call shift_field(ssr,nxfield,ny,1,1,2,n)
+    ! call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
+    ! call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
+    call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
+    call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
+    call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
+! numwfmem
+    call shift_field(ps,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(sd,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(msl,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(tcc,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(u10,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(v10,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(tt2,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(td2,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(lsprec,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(convprec,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(sshf,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(ssr,nxfield,ny,1,1,numwfmem,n)
+    call shift_field(tth,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
+    call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
+!hg
+    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
+    if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
+!hg end
+
+  endif
+
+  do i=0,nxmin1
+    do j=0,nymin1
+      surfstr(i,j,1,n)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
+    end do
+  end do
+
+  if (mp_measure_time) call mpif_mtime('readwind',1)
+
+  if ((.not.hflswitch).or.(.not.strswitch)) then
+    write(*,*) 'WARNING: No flux data contained in GRIB file ', &
+         wfname(indj)
+
+! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
+! As ECMWF has increased the model resolution, such that now the first model
+! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
+! (3rd model level in FLEXPART) for the profile method
+!***************************************************************************
+
+    do i=0,nxmin1
+      do j=0,nymin1
+        plev1=akz(3)+bkz(3)*ps(i,j,1,n)
+        pmean=0.5*(ps(i,j,1,n)+plev1)
+        tv=tth(i,j,3,n)*(1.+0.61*qvh(i,j,3,n))
+        fu=-r_air*tv/ga/pmean
+        hlev1=fu*(plev1-ps(i,j,1,n))   ! HEIGTH OF FIRST MODEL LAYER
+        ff10m= sqrt(u10(i,j,1,n)**2+v10(i,j,1,n)**2)
+        fflev1=sqrt(uuh(i,j,3)**2+vvh(i,j,3)**2)
+        call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
+             tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
+             surfstr(i,j,1,n),sshf(i,j,1,n))
+        if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
+        if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
+      end do
+    end do
+  endif
+
+
+! Assign 10 m wind to model level at eta=1.0 to have one additional model
+! level at the ground
+! Specific humidity is taken the same as at one level above
+! Temperature is taken as 2 m temperature
+!**************************************************************************
+
+  do i=0,nxmin1
+    do j=0,nymin1
+      uuh(i,j,1)=u10(i,j,1,n)
+      vvh(i,j,1)=v10(i,j,1,n)
+      qvh(i,j,1,n)=qvh(i,j,2,n)
+      tth(i,j,1,n)=tt2(i,j,1,n)
+    end do
+  end do
+
+  if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
+  if(iwmax.ne.nwz)    stop 'READWIND: NWZ NOT CONSISTENT'
+
+  return
+
+888 write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
+  write(*,*) ' #### ',wfname(indj),'                    #### '
+  write(*,*) ' #### IS NOT GRIB FORMAT !!!                  #### '
+  stop 'Execution terminated'
+
+end subroutine readwind_ecmwf
+
Index: src/verttransform_ecmwf.f90
===================================================================
--- src/verttransform_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
+++ src/verttransform_ecmwf.f90	(revision 61e07bac31955a83b2f343eb477aac0514ec1b8d)
@@ -0,0 +1,866 @@
+!**********************************************************************
+! 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/>.   *
+!**********************************************************************
+
+subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
+!                         i  i   i   i   i
+!*****************************************************************************
+!                                                                            *
+!     This subroutine transforms temperature, dew point temperature and      *
+!     wind components from eta to meter coordinates.                         *
+!     The vertical wind component is transformed from Pa/s to m/s using      *
+!     the conversion factor pinmconv.                                        *
+!     In addition, this routine calculates vertical density gradients        *
+!     needed for the parameterization of the turbulent velocities.           *
+!                                                                            *
+!     Author: A. Stohl, G. Wotawa                                            *
+!                                                                            *
+!     12 August 1996                                                         *
+!     Update: 16 January 1998                                                *
+!                                                                            *
+!     Major update: 17 February 1999                                         *
+!     by G. Wotawa                                                           *
+!                                                                            *
+!     - Vertical levels for u, v and w are put together                      *
+!     - Slope correction for vertical velocity: Modification of calculation  *
+!       procedure                                                            *
+!                                                                            *
+!*****************************************************************************
+!  Changes, Bernd C. Krueger, Feb. 2001:
+!   Variables tth and qvh (on eta coordinates) from common block
+!
+! Sabine Eckhardt, March 2007
+! added the variable cloud for use with scavenging - descr. in com_mod
+!
+! Unified ECMWF and GFS builds
+! Marian Harustak, 12.5.2017 
+!     - Renamed from verttransform to verttransform_ecmwf
+!
+!*****************************************************************************
+!                                                                            *
+! Variables:                                                                 *
+! nx,ny,nz                        field dimensions in x,y and z direction    *
+! clouds(0:nxmax,0:nymax,0:nzmax,numwfmem) cloud field for wet deposition    *
+! uu(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in x-direction [m/s]*
+! vv(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in y-direction [m/s]*
+! ww(0:nxmax,0:nymax,nzmax,numwfmem)     wind components in z-direction      *
+!                                          [deltaeta/s]                      *
+! tt(0:nxmax,0:nymax,nzmax,numwfmem)     temperature [K]                     *
+! pv(0:nxmax,0:nymax,nzmax,numwfmem)     potential voriticity (pvu)          *
+! ps(0:nxmax,0:nymax,numwfmem)           surface pressure [Pa]               *
+!                                                                            *
+!*****************************************************************************
+
+  use par_mod
+  use com_mod
+  use cmapf_mod, only: cc2gll
+!  use mpi_mod
+
+  implicit none
+
+  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: uuh,vvh,pvh
+  real,intent(in),dimension(0:nxmax-1,0:nymax-1,nwzmax) :: wwh
+
+  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
+  real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
+  real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
+  real,dimension(0:nymax-1) :: cosf
+
+  integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx
+
+  integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixm,jym,kz_inv
+  real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
+  real :: ew,dz1,dz2,dz
+  real :: xlon,ylat,xlonr,dzdx,dzdy
+  real :: dzdx1,dzdx2,dzdy1,dzdy2
+  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
+  real,parameter :: const=r_air/ga
+  real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
+
+  logical :: init = .true.
+
+  !ZHG SEP 2014 tests  
+  ! integer :: cloud_ver,cloud_min, cloud_max 
+  ! integer ::teller(5), convpteller=0, lspteller=0
+  ! real :: cloud_col_wat, cloud_water
+  !ZHG 2015 temporary variables for testing
+  ! real :: rcw(0:nxmax-1,0:nymax-1)
+  ! real :: rpc(0:nxmax-1,0:nymax-1)
+  ! character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
+  ! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
+  ! CHARACTER(LEN=3)  :: aspec
+  ! integer :: virr=0
+  real :: tot_cloud_h
+  real :: dbg_height(nzmax) 
+!ZHG
+
+!*************************************************************************
+! If verttransform is called the first time, initialize heights of the   *
+! z levels in meter. The heights are the heights of model levels, where  *
+! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
+! the vertical resolution in the z system is doubled. As reference point,*
+! the lower left corner of the grid is used.                             *
+! Unlike in the eta system, no difference between heights for u,v and    *
+! heights for w exists.                                                  *
+!*************************************************************************
+
+
+!eso measure CPU time
+!  call mpif_mtime('verttransform',0)
+
+  if (init) then
+
+! Search for a point with high surface pressure (i.e. not above significant topography)
+! Then, use this point to construct a reference z profile, to be used at all times
+!*****************************************************************************
+
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        if (ps(ix,jy,1,n).gt.100000.) then
+          ixm=ix
+          jym=jy
+          goto 3
+        endif
+      end do
+    end do
+3   continue
+
+
+    tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew*(td2(ixm,jym,1,n))/ &
+         ps(ixm,jym,1,n))
+    pold(ixm,jym)=ps(ixm,jym,1,n)
+    height(1)=0.
+
+    do kz=2,nuvz
+      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
+      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
+
+      if (abs(tv(ixm,jym)-tvold(ixm,jym)).gt.0.2) then
+        height(kz)= &
+             height(kz-1)+const*log(pold(ixm,jym)/pint(ixm,jym))* &
+             (tv(ixm,jym)-tvold(ixm,jym))/log(tv(ixm,jym)/tvold(ixm,jym))
+      else
+        height(kz)=height(kz-1)+ &
+             const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym)
+      endif
+
+      tvold(ixm,jym)=tv(ixm,jym)
+      pold(ixm,jym)=pint(ixm,jym)
+    end do
+
+
+! Determine highest levels that can be within PBL
+!************************************************
+
+    do kz=1,nz
+      if (height(kz).gt.hmixmax) then
+        nmixz=kz
+        goto 9
+      endif
+    end do
+9   continue
+
+! Do not repeat initialization of the Cartesian z grid
+!*****************************************************
+
+    init=.false.
+
+    dbg_height = height
+
+  endif
+
+
+! Loop over the whole grid
+!*************************
+
+
+  do jy=0,nymin1
+    do ix=0,nxmin1
+      tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew*(td2(ix,jy,1,n))/ &
+           ps(ix,jy,1,n))
+    enddo
+  enddo
+  pold=ps(:,:,1,n)
+  uvzlev(:,:,1)=0.
+  wzlev(:,:,1)=0.
+  rhoh(:,:,1)=pold/(r_air*tvold)
+
+
+! Compute heights of eta levels
+!******************************
+
+  do kz=2,nuvz
+    pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
+    tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
+    rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
+
+    where (abs(tv-tvold).gt.0.2) 
+      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)* &
+           (tv-tvold)/log(tv/tvold)
+    elsewhere
+      uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)*tv
+    endwhere
+
+    tvold=tv
+    pold=pint
+  end do
+
+
+  do kz=2,nwz-1
+    wzlev(:,:,kz)=(uvzlev(:,:,kz+1)+uvzlev(:,:,kz))/2.
+  end do
+  wzlev(:,:,nwz)=wzlev(:,:,nwz-1)+ &
+       uvzlev(:,:,nuvz)-uvzlev(:,:,nuvz-1)
+
+! pinmconv=(h2-h1)/(p2-p1)
+
+  pinmconv(:,:,1)=(uvzlev(:,:,2))/ &
+       ((aknew(2)+bknew(2)*ps(:,:,1,n))- &
+       (aknew(1)+bknew(1)*ps(:,:,1,n)))
+  do kz=2,nz-1
+    pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
+         ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
+         (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
+  end do
+  pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
+       ((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
+       (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
+
+! Levels, where u,v,t and q are given
+!************************************
+
+
+  uu(:,:,1,n)=uuh(:,:,1)
+  vv(:,:,1,n)=vvh(:,:,1)
+  tt(:,:,1,n)=tth(:,:,1,n)
+  qv(:,:,1,n)=qvh(:,:,1,n)
+!hg adding the cloud water 
+  if (readclouds) then
+    clwc(:,:,1,n)=clwch(:,:,1,n)
+    if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
+  end if
+!hg 
+  pv(:,:,1,n)=pvh(:,:,1)
+  rho(:,:,1,n)=rhoh(:,:,1)
+
+  uu(:,:,nz,n)=uuh(:,:,nuvz)
+  vv(:,:,nz,n)=vvh(:,:,nuvz)
+  tt(:,:,nz,n)=tth(:,:,nuvz,n)
+  qv(:,:,nz,n)=qvh(:,:,nuvz,n)
+!hg adding the cloud water
+  if (readclouds) then
+    clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
+    if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
+  end if
+!hg
+  pv(:,:,nz,n)=pvh(:,:,nuvz)
+  rho(:,:,nz,n)=rhoh(:,:,nuvz)
+
+
+  kmin=2
+  idx=kmin
+  do iz=2,nz-1
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        if(height(iz).gt.uvzlev(ix,jy,nuvz)) then
+          uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
+          vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
+          tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
+          qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
+!hg adding the cloud water
+          if (readclouds) then
+            clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
+            if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
+          end if
+!hg
+          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
+          rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
+        else
+          innuvz: do kz=idx(ix,jy),nuvz
+            if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
+                 (height(iz).le.uvzlev(ix,jy,kz))) then
+              idx(ix,jy)=kz
+              exit innuvz
+            endif
+          enddo innuvz
+        endif
+      enddo
+    enddo
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        if(height(iz).le.uvzlev(ix,jy,nuvz)) then
+          kz=idx(ix,jy)
+          dz1=height(iz)-uvzlev(ix,jy,kz-1)
+          dz2=uvzlev(ix,jy,kz)-height(iz)
+          dz=dz1+dz2
+          uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+          vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+          tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+               +tth(ix,jy,kz,n)*dz1)/dz
+          qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+               +qvh(ix,jy,kz,n)*dz1)/dz
+!hg adding the cloud water
+          if (readclouds) then
+            clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
+            if (.not.sumclouds) &
+                 &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
+          end if
+!hg
+          pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+          rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
+        endif
+      enddo
+    enddo
+  enddo
+
+
+! Levels, where w is given
+!*************************
+
+  ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
+  ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
+  kmin=2
+  idx=kmin
+  do iz=2,nz
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        inn:         do kz=idx(ix,jy),nwz
+          if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. &
+               height(iz).le.wzlev(ix,jy,kz)) then
+            idx(ix,jy)=kz
+            exit inn
+          endif
+        enddo inn
+      enddo
+    enddo
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        kz=idx(ix,jy)
+        dz1=height(iz)-wzlev(ix,jy,kz-1)
+        dz2=wzlev(ix,jy,kz)-height(iz)
+        dz=dz1+dz2
+        ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(ix,jy,kz-1)*dz2 &
+             +wwh(ix,jy,kz)*pinmconv(ix,jy,kz)*dz1)/dz
+      enddo
+    enddo
+  enddo
+
+! Compute density gradients at intermediate levels
+!*************************************************
+
+  drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/ &
+       (height(2)-height(1))
+  do kz=2,nz-1
+    drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
+         (height(kz+1)-height(kz-1))
+  end do
+  drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
+
+!    end do
+!  end do
+
+
+!****************************************************************
+! Compute slope of eta levels in windward direction and resulting
+! vertical wind correction
+!****************************************************************
+
+  do jy=1,ny-2
+    cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180)
+  enddo
+
+  kmin=2
+  idx=kmin
+  do iz=2,nz-1
+    do jy=1,ny-2
+      do ix=1,nx-2
+
+        inneta: do kz=idx(ix,jy),nz
+          if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. &
+               (height(iz).le.uvzlev(ix,jy,kz))) then
+            idx(ix,jy)=kz
+            exit inneta
+          endif
+        enddo inneta
+      enddo
+    enddo
+
+    do jy=1,ny-2
+      do ix=1,nx-2
+        kz=idx(ix,jy)
+        dz1=height(iz)-uvzlev(ix,jy,kz-1)
+        dz2=uvzlev(ix,jy,kz)-height(iz)
+        dz=dz1+dz2
+        ix1=ix-1
+        jy1=jy-1
+        ixp=ix+1
+        jyp=jy+1
+
+        dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
+        dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
+        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
+
+        dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
+        dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
+        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
+
+        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy)+dzdy*vv(ix,jy,iz,n)*dyconst)
+
+      end do
+
+    end do
+  end do
+
+! If north pole is in the domain, calculate wind velocities in polar
+! stereographic coordinates
+!*******************************************************************
+
+  if (nglobal) then
+    do iz=1,nz
+      do jy=int(switchnorthg)-2,nymin1
+        ylat=ylat0+real(jy)*dy
+        do ix=0,nxmin1
+          xlon=xlon0+real(ix)*dx
+          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
+        end do
+      end do
+    end do
+
+
+    do iz=1,nz
+
+! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+!
+!   AMSnauffer Nov 18 2004 Added check for case vv=0
+!
+      xlon=xlon0+real(nx/2-1)*dx
+      xlonr=xlon*pi/180.
+      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+           vv(nx/2-1,nymin1,iz,n)**2)
+      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
+      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
+        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
+      else
+        ddpol=pi/2-xlonr
+      endif
+      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+      xlon=180.0
+      xlonr=xlon*pi/180.
+      ylat=90.0
+      uuaux=-ffpol*sin(xlonr+ddpol)
+      vvaux=-ffpol*cos(xlonr+ddpol)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
+      jy=nymin1
+      do ix=0,nxmin1
+        uupol(ix,jy,iz,n)=uupolaux
+        vvpol(ix,jy,iz,n)=vvpolaux
+      end do
+    end do
+
+
+! Fix: Set W at pole to the zonally averaged W of the next equator-
+! ward parallel of latitude
+
+    do iz=1,nz
+      wdummy=0.
+      jy=ny-2
+      do ix=0,nxmin1
+        wdummy=wdummy+ww(ix,jy,iz,n)
+      end do
+      wdummy=wdummy/real(nx)
+      jy=nymin1
+      do ix=0,nxmin1
+        ww(ix,jy,iz,n)=wdummy
+      end do
+    end do
+
+  endif
+
+
+! If south pole is in the domain, calculate wind velocities in polar
+! stereographic coordinates
+!*******************************************************************
+
+  if (sglobal) then
+    do iz=1,nz
+      do jy=0,int(switchsouthg)+3
+        ylat=ylat0+real(jy)*dy
+        do ix=0,nxmin1
+          xlon=xlon0+real(ix)*dx
+          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
+        end do
+      end do
+    end do
+
+    do iz=1,nz
+
+! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+!
+!   AMSnauffer Nov 18 2004 Added check for case vv=0
+!
+      xlon=xlon0+real(nx/2-1)*dx
+      xlonr=xlon*pi/180.
+      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+           vv(nx/2-1,0,iz,n)**2)
+      if (vv(nx/2-1,0,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
+      else if (vv(nx/2-1,0,iz,n).gt.0.) then
+        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
+      else
+        ddpol=pi/2-xlonr
+      endif
+      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+      xlon=180.0
+      xlonr=xlon*pi/180.
+      ylat=-90.0
+      uuaux=+ffpol*sin(xlonr-ddpol)
+      vvaux=-ffpol*cos(xlonr-ddpol)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
+      jy=0
+      do ix=0,nxmin1
+        uupol(ix,jy,iz,n)=uupolaux
+        vvpol(ix,jy,iz,n)=vvpolaux
+      end do
+    end do
+
+
+! Fix: Set W at pole to the zonally averaged W of the next equator-
+! ward parallel of latitude
+
+    do iz=1,nz
+      wdummy=0.
+      jy=1
+      do ix=0,nxmin1
+        wdummy=wdummy+ww(ix,jy,iz,n)
+      end do
+      wdummy=wdummy/real(nx)
+      jy=0
+      do ix=0,nxmin1
+        ww(ix,jy,iz,n)=wdummy
+      end do
+    end do
+  endif
+
+
+!***********************************************************************************  
+  if (readclouds) then !HG METHOD
+! The method is loops all grids vertically and constructs the 3D matrix for clouds
+! Cloud top and cloud bottom gid cells are assigned as well as the total column 
+! cloud water. For precipitating grids, the type and whether it is in or below 
+! cloud scavenging are assigned with numbers 2-5 (following the old metod).
+! Distinction is done for lsp and convp though they are treated the same in regards
+! to scavenging. Also clouds that are not precipitating are defined which may be 
+! to include future cloud processing by non-precipitating-clouds. 
+!***********************************************************************************
+    write(*,*) 'Global ECMWF fields: using cloud water'
+    clw(:,:,:,n)=0.0
+!    icloud_stats(:,:,:,n)=0.0
+    ctwc(:,:,n)=0.0
+    clouds(:,:,:,n)=0
+! If water/ice are read separately into clwc and ciwc, store sum in clwc
+    if (.not.sumclouds) then 
+      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
+    end if
+    do jy=0,nymin1
+      do ix=0,nxmin1
+        lsp=lsprec(ix,jy,1,n)
+        convp=convprec(ix,jy,1,n)
+        prec=lsp+convp
+        tot_cloud_h=0
+! Find clouds in the vertical
+        do kz=1, nz-1 !go from top to bottom
+          if (clwc(ix,jy,kz,n).gt.0) then      
+! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3 
+            clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)*rho(ix,jy,kz,n))*(height(kz+1)-height(kz))
+            tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz)) 
+            
+!            icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n)          ! Column cloud water [m3/m3]
+            ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
+!            icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz))                     ! Cloud BOT height stats      [m]
+            cloudh_min=min(height(kz+1),height(kz))
+!ZHG 2015 extra for testing
+!            clh(ix,jy,kz,n)=height(kz+1)-height(kz)
+!            icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
+!            icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz))          ! Cloud TOP height            [m]
+!ZHG
+          endif
+        end do
+
+! If Precipitation. Define removal type in the vertical
+        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+
+          do kz=nz,1,-1 !go Bottom up!
+            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
+              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1) 
+              clouds(ix,jy,kz,n)=1                               ! is a cloud
+              if (lsp.ge.convp) then
+                clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
+              else
+                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
+              endif                                              ! convective or large scale
+            elseif((clw(ix,jy,kz,n).le.0) .and. (cloudh_min.ge.height(kz))) then ! is below cloud
+              if (lsp.ge.convp) then
+                clouds(ix,jy,kz,n)=5                             ! lsp dominated washout
+              else
+                clouds(ix,jy,kz,n)=4                             ! convp dominated washout
+              endif                                              ! convective or large scale 
+            endif
+
+            if (height(kz).ge. 19000) then                        ! set a max height for removal
+              clouds(ix,jy,kz,n)=0
+            endif !clw>0
+          end do !nz
+        endif ! precipitation
+      end do
+    end do
+
+! eso: copy the relevant data to clw4 to reduce amount of communicated data for MPI
+!    ctwc(:,:,n) = icloud_stats(:,:,4,n)
+
+!**************************************************************************
+  else       ! use old definitions
+!**************************************************************************
+!   create a cloud and rainout/washout field, clouds occur where rh>80%
+!   total cloudheight is stored at level 0
+    write(*,*) 'Global fields: using cloud water from Parameterization'
+    do jy=0,nymin1
+      do ix=0,nxmin1
+! OLD METHOD
+        rain_cloud_above(ix,jy)=0
+        lsp=lsprec(ix,jy,1,n)
+        convp=convprec(ix,jy,1,n)
+        cloudsh(ix,jy,n)=0
+        do kz_inv=1,nz-1
+          kz=nz-kz_inv+1
+          pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+          rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+          clouds(ix,jy,kz,n)=0
+          if (rh.gt.0.8) then ! in cloud
+            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+              rain_cloud_above(ix,jy)=1
+              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+                   height(kz)-height(kz-1)
+              if (lsp.ge.convp) then
+                clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+              else
+                clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+              endif
+            else ! no precipitation
+              clouds(ix,jy,kz,n)=1 ! cloud
+            endif
+          else ! no cloud
+            if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
+              if (lsp.ge.convp) then
+                clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+              else
+                clouds(ix,jy,kz,n)=4 ! convp dominated washout
+              endif
+            endif
+          endif
+        end do
+!END OLD METHOD
+      end do
+    end do
+  endif !readclouds
+
+
+     !********* TEST ***************
+     ! WRITE OUT SOME TEST VARIABLES
+     !********* TEST ************'**
+!teller(:)=0
+!virr=virr+1
+!WRITE(aspec, '(i3.3)'), virr
+
+!if (readclouds) then
+!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement.txt'
+!else
+!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement_old.txt'
+!endif
+!
+!OPEN(UNIT=118, FILE=fnameH,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!do kz_inv=1,nz-1
+!  kz=nz-kz_inv+1
+!  !kz=91
+!  do jy=0,nymin1
+!     do ix=0,nxmin1
+!          if (clouds(ix,jy,kz,n).eq.1) teller(1)=teller(1)+1 ! no precipitation cloud
+!          if (clouds(ix,jy,kz,n).eq.2) teller(2)=teller(2)+1 ! convp dominated rainout
+!          if (clouds(ix,jy,kz,n).eq.3) teller(3)=teller(3)+1 ! lsp dominated rainout
+!          if (clouds(ix,jy,kz,n).eq.4) teller(4)=teller(4)+1 ! convp dominated washout
+!          if (clouds(ix,jy,kz,n).eq.5) teller(5)=teller(5)+1 ! lsp dominated washout
+!          
+!        !  write(*,*) height(kz),teller
+!     end do
+!  end do
+!  write(118,*) height(kz),teller
+!  teller(:)=0
+!end do
+!teller(:)=0
+!write(*,*) teller 
+!write(*,*) aspec
+!
+!fnameA=trim(zhgpath)//trim(aspec)//'cloudV.txt'
+!fnameB=trim(zhgpath)//trim(aspec)//'cloudT.txt'
+!fnameC=trim(zhgpath)//trim(aspec)//'cloudB.txt'
+!fnameD=trim(zhgpath)//trim(aspec)//'cloudW.txt'
+!fnameE=trim(zhgpath)//trim(aspec)//'old_cloudV.txt'
+!fnameF=trim(zhgpath)//trim(aspec)//'lsp.txt'
+!fnameG=trim(zhgpath)//trim(aspec)//'convp.txt'
+!if (readclouds) then
+!OPEN(UNIT=111, FILE=fnameA,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=112, FILE=fnameB,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=113, FILE=fnameC,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=114, FILE=fnameD,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!else
+!OPEN(UNIT=115, FILE=fnameE,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=116, FILE=fnameF,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!OPEN(UNIT=117, FILE=fnameG,FORM='FORMATTED',STATUS = 'UNKNOWN')
+!endif
+!
+!do ix=0,nxmin1
+!if (readclouds) then
+!write(111,*) (icloud_stats(ix,jy,1,n),jy=0,nymin1)
+!write(112,*) (icloud_stats(ix,jy,2,n),jy=0,nymin1)
+!write(113,*) (icloud_stats(ix,jy,3,n),jy=0,nymin1)
+!write(114,*) (icloud_stats(ix,jy,4,n),jy=0,nymin1)
+!else
+!write(115,*) (cloudsh(ix,jy,n),jy=0,nymin1)    !integer
+!write(116,*) (lsprec(ix,jy,1,n),jy=0,nymin1)   !7.83691406E-02 
+!write(117,*) (convprec(ix,jy,1,n),jy=0,nymin1) !5.38330078E-02
+!endif
+!end do
+!
+!if (readclouds) then
+!CLOSE(111)
+!CLOSE(112)
+!CLOSE(113)
+!CLOSE(114)
+!else
+!CLOSE(115)
+!CLOSE(116)
+!CLOSE(117)
+!endif
+!
+!END ********* TEST *************** END
+! WRITE OUT SOME TEST VARIABLES
+!END ********* TEST *************** END
+
+
+! PS 2012
+!      lsp=lsprec(ix,jy,1,n)
+!      convp=convprec(ix,jy,1,n)
+!          prec=lsp+convp
+!          if (lsp.gt.convp) then !  prectype='lsp'
+!            lconvectprec = .false.
+!          else ! prectype='cp'
+!            lconvectprec = .true.
+!           endif
+!      else ! windfields does not contain cloud data 
+!          rhmin = 0.90 ! standard condition for presence of clouds
+!PS       note that original by Sabine Eckhart was 80%
+!PS       however, for T<-20 C we consider saturation over ice
+!PS       so I think 90% should be enough          
+!          icloudbot(ix,jy,n)=icmv
+!          icloudtop=icmv ! this is just a local variable
+!98        do kz=1,nz
+!            pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+!            rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+!ps            if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
+!            if (rh .gt. rhmin) then
+!              if (icloudbot(ix,jy,n) .eq. icmv) then
+!                icloudbot(ix,jy,n)=nint(height(kz))
+!              endif
+!              icloudtop=nint(height(kz)) ! use int to save memory
+!            endif
+    ! end do
+!PS try to get a cloud thicker than 50 m 
+!PS if there is at least .01 mm/h  - changed to 0.002 and put into
+!PS parameter precpmin        
+!          if ((icloudbot(ix,jy,n) .eq. icmv .or. &
+!              icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
+!              prec .gt. precmin) then
+!            rhmin = rhmin - 0.05
+!            if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
+!   end if
+
+!PS is based on looking at a limited set of comparison data
+!          if (lconvectprec .and. icloudtop .lt. 6000 .and. &
+!             prec .gt. precmin) then 
+!
+!            if (convp .lt. 0.1) then
+!              icloudbot(ix,jy,n) = 500
+!              icloudtop =         8000
+!            else
+!              icloudbot(ix,jy,n) = 0
+!              icloudtop =      10000
+!            endif
+!          endif
+!          if (icloudtop .ne. icmv) then
+!            icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
+!          else
+!            icloudthck(ix,jy,n) = icmv
+!          endif
+!PS  get rid of too thin clouds      
+!          if (icloudthck(ix,jy,n) .lt. 50) then
+!            icloudbot(ix,jy,n)=icmv
+!            icloudthck(ix,jy,n)=icmv
+!          endif
+!hg__________________________________
+!           rcw(ix,jy)=2E-7*prec**0.36
+!           rpc(ix,jy)=prec
+!hg end______________________________
+
+!      endif !hg read clouds
+
+
+
+
+!eso measure CPU time
+!  call mpif_mtime('verttransform',1)
+
+!eso print out the same measure as in Leo's routine
+    ! write(*,*) 'verttransform: ', &
+    !      sum(tt(:,:,:,n)*tt(:,:,:,n)), &
+    !      sum(uu(:,:,:,n)*uu(:,:,:,n)),sum(vv(:,:,:,n)*vv(:,:,:,n)),&
+    !      sum(qv(:,:,:,n)*qv(:,:,:,n)),sum(pv(:,:,:,n)*pv(:,:,:,n)),&
+    !      sum(rho(:,:,:,n)*rho(:,:,:,n)),sum(drhodz(:,:,:,n)*drhodz(:,:,:,n)),&
+    !      sum(ww(:,:,:,n)*ww(:,:,:,n)), &
+    !      sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
+end subroutine verttransform_ecmwf
+
