mode_read_grib.F90

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #####################
MODULE MODE_READ_GRIB
!     #####################
!-------------------------------------------------------------------
!
USE MODI_ABOR1_SFX
USE GRIB_API
!
USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
USE PARKIND1  ,ONLY : JPRB
!
CONTAINS
!-------------------------------------------------------------------
!     ####################
      SUBROUTINE MAKE_GRIB_INDEX(HGRIB,HINMODEL)
!     ####################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_GRID_GRIB, ONLY : CGRIB_FILE, NIDX, NGRIB_VERSION
!
IMPLICIT NONE
!
 CHARACTER(LEN=*), INTENT(IN) :: HGRIB
 CHARACTER(LEN=6), INTENT(IN)    :: HINMODEL  ! Grib originating model
!
INTEGER(KIND=kindOfInt) :: IRET
REAL(KIND=JPRB) :: ZHOOK_HANDLE
INTEGER, DIMENSION(:), ALLOCATABLE  :: IVERSION
INTEGER :: ILEN
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:MAKE_GRIB_INDEX',0,ZHOOK_HANDLE)
!
IF (CGRIB_FILE==HGRIB) CALL DR_HOOK('MODE_READ_GRIB:MAKE_GRIB_INDEX',1,ZHOOK_HANDLE)
IF (CGRIB_FILE==HGRIB) RETURN
!
CGRIB_FILE=HGRIB
!
 CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'editionNumber',IRET)
IF (IRET/=0) THEN
    CALL ABOR1_SFX("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while searching edition number")
ELSE
  call GRIB_INDEX_GET_SIZE(NIDX,'editionNumber',ILEN)
  ALLOCATE(IVERSION(ILEN))
  CALL GRIB_INDEX_GET(NIDX, 'editionNumber', IVERSION , IRET)
  NGRIB_VERSION=IVERSION(1)

  CALL GRIB_INDEX_RELEASE(NIDX,IRET)
  IF (IRET/=0) CALL ABOR1_SFX("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while deleting the grib index")
ENDIF

!test version
IF (NGRIB_VERSION==1) THEN
 CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'indicatorOfParameter',IRET)
ELSEIF (NGRIB_VERSION==2) THEN
 IF(HINMODEL=='ARPEGE') THEN
   print*,"CALL GRIB_INDEX_CREATE"
   CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'parameterNumber',IRET)
print*,IRET
 ELSE
   CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'paramId',IRET)        
 ENDIF
ENDIF
IF (IRET/=0) CALL ABOR1_SFX("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while creating the grib index")
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:MAKE_GRIB_INDEX',1,ZHOOK_HANDLE)
!
END SUBROUTINE MAKE_GRIB_INDEX
!-------------------------------------------------------------------
!     ####################
      SUBROUTINE CLEAR_GRIB_INDEX
!     ####################
!
USE MODD_GRID_GRIB, ONLY : CGRIB_FILE, NIDX
!
IMPLICIT NONE
!
INTEGER(KIND=kindOfInt) :: IRET
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:CLEAR_GRIB_INDEX',0,ZHOOK_HANDLE)
!
IF (CGRIB_FILE.NE."") THEN
  CGRIB_FILE=""
  CALL GRIB_INDEX_RELEASE(NIDX,IRET)
  IF (IRET/=0) CALL ABOR1_SFX("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while deleting the grib index")
ENDIF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:CLEAR_GRIB_INDEX',1,ZHOOK_HANDLE)
!
END SUBROUTINE CLEAR_GRIB_INDEX
!-------------------------------------------------------------------
!     ####################
      SUBROUTINE GET_GRIB_MESSAGE(KLUOUT,KLTYPE,KLEV1,KLEV2,KGRIB,KFOUND,HTYPELEVEL,PLEV1,PLEV2)
!     ####################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_GRID_GRIB, ONLY : NIDX
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KLUOUT
INTEGER, INTENT(INOUT)  :: KLTYPE
INTEGER, INTENT(INOUT)  :: KLEV1
INTEGER, INTENT(INOUT)  :: KLEV2
INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KGRIB
INTEGER, INTENT(OUT) :: KFOUND
CHARACTER(LEN=*), INTENT(INOUT), OPTIONAL  :: HTYPELEVEL      ! TypeOfLevel JPMODIF
REAL, INTENT(INOUT), OPTIONAL  :: PLEV1      ! top level of soil
REAL, INTENT(INOUT), OPTIONAL  :: PLEV2      ! Bottom level of soil
!

INTEGER :: ILTYPE
INTEGER :: ILEV1
INTEGER :: ILEV2
CHARACTER(LEN=50) :: YTYPELEVEL      ! TypeOfLevel JPMODIF
REAL :: ZLEV1,ZLEV2
INTEGER(KIND=kindOfInt) :: IRET
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:GET_GRIB_MESSAGE',0,ZHOOK_HANDLE)
!
IRET = 0
KFOUND=0
!
DO WHILE (IRET /= GRIB_END_OF_INDEX .AND. KFOUND/=3)
print*,"===============new message=============="
  !
  IRET = 0
  KFOUND=0
  !
  IF (KLTYPE/=-2) THEN
    CALL GRIB_GET(KGRIB,'indicatorOfTypeOfLevel',ILTYPE,IRET)
print*,IRET
    IF(IRET/=0) THEN
      CALL GRIB_GET(KGRIB,'typeOfFirstFixedSurface',ILTYPE,IRET)
    ENDIF
print*,IRET
    CALL TEST_IRET(KLUOUT,ILTYPE,KLTYPE,IRET)
    print*,"ILTYPE,KLTYPE,IRET",ILTYPE,KLTYPE,IRET
  ELSE
     IF (PRESENT(HTYPELEVEL)) THEN
        CALL GRIB_GET(KGRIB,'typeOfLevel',YTYPELEVEL,IRET)
        CALL TEST_IRET_STR(KLUOUT,TRIM(YTYPELEVEL),HTYPELEVEL,IRET)
     ENDIF
  ENDIF
  !
  IF (IRET.EQ.0) THEN
    !
    KFOUND = KFOUND + 1
    !
    IF (KLEV1/=-2) THEN
      CALL GRIB_GET(KGRIB,'topLevel',ILEV1,IRET)
      CALL TEST_IRET(KLUOUT,ILEV1,KLEV1,IRET)
    ELSE IF (PRESENT(PLEV1)) THEN  !JP
      CALL GRIB_GET(KGRIB,'topLevel',ZLEV1,IRET)
      CALL TEST_IRET_FLOAT(KLUOUT,ZLEV1,PLEV1,IRET)

    ENDIF
    !
print*,KFOUND
    IF (IRET.EQ.0) THEN
      !
      KFOUND = KFOUND + 1
      !
      IF (KLEV2/=-2) THEN
        CALL GRIB_GET(KGRIB,'bottomLevel',ILEV2,IRET)
        CALL TEST_IRET(KLUOUT,ILEV2,KLEV2,IRET)
      ELSE IF (PRESENT(PLEV2)) THEN  !JP
        CALL GRIB_GET(KGRIB,'bottomLevel',ZLEV2,IRET)
        CALL TEST_IRET_FLOAT(KLUOUT,ZLEV2,PLEV2,IRET)
      ENDIF
      !
      IF (IRET.EQ.0) KFOUND = KFOUND + 1
      !
print*,KFOUND
    ENDIF
    !
  ENDIF
  !
  IF (KFOUND.NE.3) THEN
    CALL GRIB_RELEASE(KGRIB)
    CALL GRIB_NEW_FROM_INDEX(NIDX,KGRIB,IRET)
  ENDIF
  !
ENDDO
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:GET_GRIB_MESSAGE',1,ZHOOK_HANDLE)
!
CONTAINS
!
!       ##############
        SUBROUTINE TEST_IRET(KLUOUT,VAL1,VAL0,KRET)
!       ##############
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
INTEGER, INTENT(IN) :: VAL1
INTEGER, INTENT(INOUT) :: VAL0
INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KRET   ! number of the message researched
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',0,ZHOOK_HANDLE)
!
IF (KRET > 0) THEN
  WRITE (KLUOUT,'(A)')' | Error encountered in the Grib file, skipping field'
ELSE IF (KRET == -6) THEN
  WRITE (KLUOUT,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
ELSEIF (VAL1 /= VAL0) THEN
  IF (VAL0 == -1) THEN
    VAL0 = VAL1
  ELSE
    KRET=1
  ENDIF
ENDIF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',1,ZHOOK_HANDLE)
END SUBROUTINE TEST_IRET

!       ##############
        SUBROUTINE TEST_IRET_STR(KLUOUT,VAL1,VAL0,KRET)
!       ##############
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
CHARACTER(LEN=*), INTENT(IN) :: VAL1
CHARACTER(LEN=*), INTENT(INOUT) :: VAL0
INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KRET   ! number of the message researched
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',0,ZHOOK_HANDLE)
!
IF (KRET > 0) THEN
  WRITE (KLUOUT,'(A)')' | Error encountered in the Grib file, skipping field'
ELSE IF (KRET == -6) THEN
  WRITE (KLUOUT,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
ELSEIF (VAL1 /= VAL0) THEN
  IF (VAL0 == '') THEN
    VAL0 = VAL1
  ELSE
    KRET=1
  ENDIF
ENDIF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',1,ZHOOK_HANDLE)
END SUBROUTINE TEST_IRET_STR

!       ##############
        SUBROUTINE TEST_IRET_FLOAT(KLUOUT,VAL1,VAL0,KRET)
!       ##############
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
REAL, INTENT(IN) :: VAL1
REAL, INTENT(INOUT) :: VAL0
INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KRET   ! number of the message researched
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',0,ZHOOK_HANDLE)
!
IF (KRET > 0) THEN
  WRITE (KLUOUT,'(A)')' | Error encountered in the Grib file, skipping field'
ELSE IF (KRET == -6) THEN
  WRITE (KLUOUT,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
ELSEIF (VAL1 /= VAL0) THEN
  IF (VAL0 == -1.0) THEN
    VAL0 = VAL1
  ELSE
    KRET=1
  ENDIF
ENDIF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',1,ZHOOK_HANDLE)
END SUBROUTINE TEST_IRET_FLOAT

!
END SUBROUTINE GET_GRIB_MESSAGE
!-------------------------------------------------------------------
!     ####################
      SUBROUTINE READ_GRIB(HGRIB,HINMODEL,KLUOUT,KPARAM,KRET,PFIELD,KLTYPE,KLEV1,KLEV2,KPARAM2, &
         KLTYPE2,HTYPELEVEL,PLEV1,PLEV2)
!     ####################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_GRID_GRIB, ONLY : NIDX, NGRIB_VERSION
!
IMPLICIT NONE
!
 CHARACTER(LEN=*), INTENT(IN) :: HGRIB      ! name of the GRIB file
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
INTEGER, INTENT(IN) :: KLUOUT
INTEGER,INTENT(IN)    :: KPARAM ! Parameter to read
INTEGER(KIND=kindOfInt), INTENT(OUT)  :: KRET
REAL, DIMENSION(:), POINTER :: PFIELD
INTEGER,INTENT(INOUT), OPTIONAL :: KLTYPE ! Level type
INTEGER,INTENT(INOUT), OPTIONAL :: KLEV1  ! Level parameter 1
INTEGER,INTENT(INOUT), OPTIONAL :: KLEV2  ! Level parameter 2
INTEGER, INTENT(INOUT), OPTIONAL :: KPARAM2
INTEGER,INTENT(INOUT), OPTIONAL :: KLTYPE2 ! Level type
CHARACTER(LEN=*), INTENT(INOUT), OPTIONAL :: HTYPELEVEL
!
REAL, INTENT(INOUT), OPTIONAL :: PLEV1,PLEV2
!
INTEGER :: ILTYPE, ILEV1, ILEV2,ILTYPE2
INTEGER(KIND=kindOfInt) :: IGRIB
INTEGER :: ISIZE, IFOUND
REAL(KIND=JPRB) :: ZHOOK_HANDLE
REAL :: ZLEV1,ZLEV2
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB',0,ZHOOK_HANDLE)
!
ILTYPE=-2
IF (PRESENT(KLTYPE)) ILTYPE=KLTYPE
ILEV1=-2
IF (PRESENT(KLEV1)) ILEV1=KLEV1
ILEV2=-2
IF (PRESENT(KLEV2)) ILEV2=KLEV2
ZLEV1=-2.0
IF (PRESENT(PLEV1)) ZLEV1=PLEV1
ZLEV2=-2.0
IF (PRESENT(PLEV2)) ZLEV2=PLEV2
!
 CALL MAKE_GRIB_INDEX(HGRIB,HINMODEL)
!
IFOUND=0
KRET=0
!
IF (NGRIB_VERSION == 1) THEN
 CALL GRIB_INDEX_SELECT(NIDX,'indicatorOfParameter',KPARAM,KRET)
ELSEIF (NGRIB_VERSION == 2) THEN
 IF (HINMODEL=='ARPEGE') THEN
  print*,"GRIB_INDEX_SELECT :",KPARAM
   CALL GRIB_INDEX_SELECT(NIDX,'parameterNumber',KPARAM,KRET)
print*,KRET
 ELSE
   CALL GRIB_INDEX_SELECT(NIDX,'paramId',KPARAM,KRET)         
 ENDIF
END IF
 CALL GRIB_NEW_FROM_INDEX(NIDX,IGRIB,KRET)


WRITE (KLUOUT,*) 'READ_GRIB GRIB_NEW_FROM_INDEX ',KPARAM,IGRIB,KRET
IF (KRET.EQ.0) THEN
   IF (PRESENT(HTYPELEVEL)) THEN
      CALL GET_GRIB_MESSAGE(KLUOUT,ILTYPE,ILEV1,ILEV2,IGRIB,IFOUND,HTYPELEVEL,ZLEV1,ZLEV2)
   ELSE
print*,"CALL GET_GRIB_MESSAGE"
      CALL GET_GRIB_MESSAGE(KLUOUT,ILTYPE,ILEV1,ILEV2,IGRIB,IFOUND)
   ENDIF
ENDIF
!

WRITE (KLUOUT,*) 'READ_GRIB GRIB_NEW_FROM_INDEX ',KPARAM,IGRIB,KRET,IFOUND

IF (PRESENT(KPARAM2)) THEN
  IF (IFOUND/=3) THEN
    IF (HINMODEL=='ARPEGE' .AND.NGRIB_VERSION == 2) THEN
      CALL GRIB_INDEX_SELECT(NIDX,'parameterNumber',KPARAM2,KRET)
    ELSE
      CALL GRIB_INDEX_SELECT(NIDX,'indicatorOfParameter',KPARAM2,KRET)
    ENDIF
    CALL GRIB_NEW_FROM_INDEX(NIDX,IGRIB,KRET)
    IF (KRET.EQ.0) THEN
      ILTYPE=-2
      IF (PRESENT(KLTYPE)) ILTYPE=KLTYPE
      IF (PRESENT(KLTYPE2)) ILTYPE=KLTYPE2
      CALL GET_GRIB_MESSAGE(KLUOUT,ILTYPE,ILEV1,ILEV2,IGRIB,IFOUND)
    ENDIF
  ELSE
    KPARAM2 = KPARAM
  ENDIF
ENDIF
!
print*,"IFOUND=",IFOUND

IF (IFOUND==3) THEN
  !
  IF (PRESENT(KLTYPE)) KLTYPE = ILTYPE
  IF (PRESENT(KLEV1))  KLEV1  = ILEV1
  IF (PRESENT(KLEV2))  KLEV2  = ILEV2
  !

  IF (.NOT.ASSOCIATED(PFIELD)) THEN
    CALL GRIB_GET_SIZE(IGRIB,'values',ISIZE,KRET)
    IF (KRET.NE.0) CALL ABOR1_SFX("MODE_READ_GRIB:READ_GRIB: Problem getting size of values")
    ALLOCATE(PFIELD(ISIZE))
  ENDIF
  !
  CALL GRIB_GET(IGRIB,'values',PFIELD,KRET)
  IF (KRET.NE.0) CALL ABOR1_SFX("MODE_READ_GRIB:READ_GRIB: Problem getting values")
  CALL GRIB_RELEASE(IGRIB,KRET)
  IF (KRET.NE.0) CALL ABOR1_SFX("MODE_READ_GRIB:READ_GRIB: Problem releasing memory")
  !
ELSE
  !
  KRET = 1
  !
ENDIF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB
!
!-------------------------------------------------------------------
!-------------------------------------------------------------------
!     ####################
      SUBROUTINE READ_GRIB_LAND_MASK(HGRIB,KLUOUT,HINMODEL,PMASK)
!     ####################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_GRID_GRIB, ONLY : NGRIB_VERSION
!
IMPLICIT NONE
!
CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), POINTER       :: PMASK     ! Land mask
!
INTEGER(KIND=kindOfInt)                 :: IRET      ! return code
INTEGER                           :: ILTYPE,ILTYPE2    ! leveltype
INTEGER                           :: ILEV      ! level
INTEGER :: INUM_ZS,ISIZE,ICOUNT,JLOOP,IPARAM,IGRIB,IPARAM2
CHARACTER(LEN=24) :: YLTYPELU
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
print*,"HINMODEL=",HINMODEL
print*,"NGRIB_VERSION=",NGRIB_VERSION
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_LAND_MASK',0,ZHOOK_HANDLE)
WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_LAND_MASK: | Reading land mask from ',HINMODEL
!
SELECT CASE (HINMODEL)
  CASE ('ECMWF ')
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,172,IRET,PMASK) 
  CASE ('NCEP  ')
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,172,IRET,PMASK) 
  CASE ('ARPEGE','ALADIN','MOCAGE')
    print*,"NGRIB_VERSION=",NGRIB_VERSION
    IF(HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
      ILTYPE=1
      CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,0,IRET,PMASK,KLTYPE=ILTYPE)   
    ELSE
      CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,81,IRET,PMASK)          
    ENDIF
print*,"NGRIB_VERSION=",NGRIB_VERSION
  CASE ('HIRLAM')        
    ILTYPE=105
    ILEV  =0    
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,81,IRET,PMASK,KLTYPE=ILTYPE,KLEV1=ILEV)            
  CASE DEFAULT
    CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_LAND_MASK: OPTION NOT SUPPORTED '//HINMODEL)
END SELECT
!
IF (IRET /= 0) THEN
  CALL ABOR1_SFX('MODE_READ_GRIB: LAND SEA MASK MISSING (READ_GRIB_LAND_MASK)')
END IF
!
WHERE (PMASK>0.5)
  PMASK = 1.
ELSEWHERE
  PMASK = 0.
END WHERE
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_LAND_MASK',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_LAND_MASK
!-------------------------------------------------------------------
!     ############################
      SUBROUTINE READ_GRIB_ZS(HGRIB,KLUOUT,HINMODEL,PZS)
!     ############################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_CSTS,       ONLY : XG
USE MODD_GRID_GRIB, ONLY : NGRIB_VERSION
!
IMPLICIT NONE
!
 CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), POINTER       :: PZS       ! 
!
INTEGER(KIND=kindOfInt)                           :: IRET      ! return code

INTEGER :: ISIZE,JLOOP,ICOUNT
CHARACTER(LEN=24) :: YLTYPELU

REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
!* Read orography
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_ZS',0,ZHOOK_HANDLE)
WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_ZS: | Reading orography from ',HINMODEL
!
SELECT CASE (HINMODEL)
  CASE ('ECMWF ')
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,129,IRET,PZS) 
  CASE ('NCEP  ')
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228002,IRET,PZS)               
  CASE ('ARPEGE','MOCAGE')
    IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
      CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,5,IRET,PZS)
    ELSE
      CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,8,IRET,PZS)               
    ENDIF
  CASE ('HIRLAM','ALADIN')
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,6,IRET,PZS)                  
  CASE DEFAULT
    CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_ZS:OPTION NOT SUPPORTED '//HINMODEL)
END SELECT
!
IF (IRET /= 0) THEN
  CALL ABOR1_SFX('MODE_READ_GRIB: OROGRAPHY MISSING (READ_GRIB_ZS_LAND)')
END IF
!
! Datas given in archives are multiplied by the gravity acceleration
PZS = PZS / XG
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_ZS',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_ZS
!-------------------------------------------------------------------
!     ############################
      SUBROUTINE READ_GRIB_ZS_LAND(HGRIB,KLUOUT,HINMODEL,PMASK,PZSL)
!     ############################
!
IMPLICIT NONE
!
 CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), INTENT(IN)    :: PMASK     ! grib land mask
REAL, DIMENSION(:), POINTER       :: PZSL      ! 
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_ZS_LAND',0,ZHOOK_HANDLE)
!
 CALL READ_GRIB_ZS(HGRIB,KLUOUT,HINMODEL,PZSL)
!
IF (SIZE(PMASK)==SIZE(PZSL)) &
  WHERE (PMASK(:)/=1.) PZSL = 0.
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_ZS_LAND',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_ZS_LAND
!-------------------------------------------------------------------
!     ############################
      SUBROUTINE READ_GRIB_ZS_SEA(HGRIB,KLUOUT,HINMODEL,PMASK,PZSS)
!     ############################
!
IMPLICIT NONE
!
 CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), INTENT(IN)    :: PMASK     ! grib land mask
REAL, DIMENSION(:), POINTER       :: PZSS      ! 
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_ZS_SEA',0,ZHOOK_HANDLE)
!
  CALL READ_GRIB_ZS(HGRIB,KLUOUT,HINMODEL,PZSS)
!
IF (SIZE(PMASK)==SIZE(PZSS)) &
  WHERE (PMASK(:)/=0.) PZSS = 0.
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_ZS_SEA',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_ZS_SEA
!-------------------------------------------------------------------
!     ###########################
      SUBROUTINE READ_GRIB_T(HGRIB,KLUOUT,HINMODEL,PT)
!     ###########################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_GRID_GRIB, ONLY : NGRIB_VERSION

IMPLICIT NONE
!
 CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), POINTER       :: PT        ! 
!
INTEGER(KIND=kindOfInt)                           :: IRET      ! return code
INTEGER                           :: ILTYPE    ! type of level (Grib code table 3)
INTEGER                           :: ILEV      ! level definition
INTEGER                           :: IPARAM    ! number of ParamId or IndicatorOfParameter or ParameterNumber 
REAL(KIND=JPRB) :: ZHOOK_HANDLE
CHARACTER(LEN=7)   :: YTYPELEVEL  ! Type of searched level 

!-------------------------------------------------------------------
!* Read surface temperature
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_T',0,ZHOOK_HANDLE)
WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_T: | Reading surface temperature'
!
SELECT CASE (HINMODEL)
  CASE ('ECMWF ')
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,139,IRET,PT)    

  CASE ('ARPEGE','ALADIN','MOCAGE')
    ILEV=0
    IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
      IPARAM=18
      ILTYPE=1
    ELSE
      IPARAM=11
      ILTYPE=111
    ENDIF
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV) 
    IF (IRET /= 0) THEN
      ILTYPE=1
      CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,PT,KLTYPE=ILTYPE) 
      IF (IRET /= 0) THEN
         ILTYPE=105
        CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV)   
      ENDIF        
    END IF

  CASE ('HIRLAM ')
    WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_T: | Reading surface temperature tile 4'
     ILTYPE=105
     ILEV=904
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,11,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV)

  CASE ('NCEP  ')
    YTYPELEVEL='surface'
    ILEV=0
    ILTYPE=-2
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,130,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV,KLEV2=ILEV,HTYPELEVEL=YTYPELEVEL) 

  CASE DEFAULT
    CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_T:OPTION NOT SUPPORTED '//HINMODEL)
END SELECT
!
IF (IRET /= 0) THEN
  CALL ABOR1_SFX('MODE_READ_GRIB: SURFACE TEMPERATURE MISSING (READ_GRIB_T)')
END IF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_T',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_T
!-------------------------------------------------------------------
!     ###########################
      SUBROUTINE READ_GRIB_TS(HGRIB,KLUOUT,HINMODEL,PMASK,PTS)
!     ###########################
!
USE MODD_SURF_PAR,   ONLY : XUNDEF
!
IMPLICIT NONE
!
 CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), INTENT(IN)    :: PMASK     ! grib land mask
REAL, DIMENSION(:), POINTER       :: PTS       ! 
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TS',0,ZHOOK_HANDLE)
!
 CALL READ_GRIB_T(HGRIB,KLUOUT,HINMODEL,PTS)
!
IF (SIZE(PMASK)==SIZE(PTS)) WHERE (PMASK(:)/=1.) PTS = XUNDEF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TS',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_TS
!-------------------------------------------------------------------
!     ###########################
      SUBROUTINE READ_GRIB_SST(HGRIB,KLUOUT,HINMODEL,PMASK,PSST)
!     ###########################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_SURF_PAR,   ONLY : XUNDEF
!
IMPLICIT NONE
!
 CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), INTENT(IN)    :: PMASK     ! grib land mask
REAL, DIMENSION(:), POINTER       :: PSST      ! 
!
INTEGER :: IRET
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_SST',0,ZHOOK_HANDLE)
!
SELECT CASE (HINMODEL)
  CASE ('ECMWF ')
     CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,34,IRET,PSST)
     IF (IRET /= 0) CALL READ_GRIB_T(HGRIB,KLUOUT,HINMODEL,PSST)
  CASE ('ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP  ')
    CALL READ_GRIB_T(HGRIB,KLUOUT,HINMODEL,PSST)
   CASE DEFAULT
     CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_SST:OPTION NOT SUPPORTED '//HINMODEL)    
END SELECT
!
IF (SIZE(PMASK)==SIZE(PSST)) WHERE (PMASK(:)/=0.) PSST = XUNDEF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_SST',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_SST
!-------------------------------------------------------------------
!     ###########################
      SUBROUTINE READ_GRIB_TSWATER(HGRIB,KLUOUT,HINMODEL,PMASK,PTS)
!     ###########################
!
!    MODIFICATIONS
!    Julien Pergaud(via J.Escobar) 10/2018 : correction for GFS undef=9999. temperature values
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_SURF_PAR,   ONLY : XUNDEF
!
IMPLICIT NONE
!
 CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), INTENT(IN)    :: PMASK     ! grib land mask
REAL, DIMENSION(:), POINTER       :: PTS     ! 
!
INTEGER :: IRET
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TSWATER',0,ZHOOK_HANDLE)
!
SELECT CASE (HINMODEL)
  CASE ('ECMWF ')
     CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,3080,IRET,PTS)
     IF (IRET /= 0) CALL READ_GRIB_T2(HGRIB,KLUOUT,HINMODEL,PMASK,PTS)
  CASE ('ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP  ')
    CALL READ_GRIB_T2(HGRIB,KLUOUT,HINMODEL,PMASK,PTS)
   CASE DEFAULT
     CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_TSWATER:OPTION NOT SUPPORTED '//HINMODEL)    
END SELECT
!
IF (SIZE(PMASK)==SIZE(PTS)) WHERE ((PMASK(:)/=1.) .OR. ((PMASK(:)==1.) .AND.(PTS(:)==9999.))) PTS = XUNDEF

!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TSWATER',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_TSWATER
!-------------------------------------------------------------------
!     ###########################
      SUBROUTINE READ_GRIB_T2(HGRIB,KLUOUT,HINMODEL,PMASK,PT2)
!     ###########################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
USE MODD_SURF_PAR,   ONLY : XUNDEF
USE MODD_GRID_GRIB, ONLY : NGRIB_VERSION
!
IMPLICIT NONE
 !
 CHARACTER(LEN=*),   INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), INTENT(IN)    :: PMASK     ! grib land mask
REAL, DIMENSION(:), POINTER       :: PT2       ! 
!
INTEGER(KIND=kindOfInt)                           :: IRET
INTEGER                           :: ILTYPE, ILEV    ! type of level (Grib code table 3)
INTEGER                           :: IPARAM    ! number of ParamId or IndicatorOfParameter or ParameterNumber
REAL(KIND=JPRB) :: ZHOOK_HANDLE
CHARACTER(LEN=19)   :: YTYPELEVEL  ! Type of searched level 
REAL :: ZLEV1,ZLEV2
!-------------------------------------------------------------------
!* Read deep soil temperature
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_T2',0,ZHOOK_HANDLE)
WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_T2: | Reading deep soil temperature'
!
SELECT CASE (HINMODEL)
  CASE ('ECMWF ')
     CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,170,IRET,PT2)   
  CASE ('ARPEGE','ALADIN','MOCAGE')
     IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
      IPARAM=18
      ILTYPE=106
    ELSE
      IPARAM=11
      ILTYPE=111
    ENDIF
     CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,PT2,KLTYPE=ILTYPE) 
    IF (IRET /= 0) THEN
       ILTYPE=105
      CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,PT2,KLTYPE=ILTYPE)   
    ENDIF
  CASE ('NCEP  ') ! toto
     YTYPELEVEL='depthBelowLandLayer'
     ILTYPE=-2
     ZLEV1=0.1
     ZLEV2=0.4
     CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228139,IRET,PT2,KLTYPE=ILTYPE,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)  
  CASE ('HIRLAM ')
     ILTYPE=105
     ILEV=954
    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,11,IRET,PT2,KLTYPE=ILTYPE,KLEV1=ILEV)      
   CASE DEFAULT
     CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_T2:OPTION NOT SUPPORTED '//HINMODEL)
END SELECT
!
IF (IRET /= 0) THEN
  CALL ABOR1_SFX('MODE_READ_GRIB: DEEP SOIL TEMPERATURE MISSING (READ_GRIB_T2)')
END IF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_T2',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_T2
!-------------------------------------------------------------------
!-------------------------------------------------------------------
!     ###########################
      SUBROUTINE READ_GRIB_T2_LAND(HGRIB,KLUOUT,HINMODEL,PMASK,ZFIELD)
!     ###########################
!
!    MODIFICATIONS
!    Julien Pergaud(via J.Escobar) 10/2018 : correction for GFS undef=9999. temperature values
!
USE MODD_SURF_PAR,   ONLY : XUNDEF
!
IMPLICIT NONE
!
 CHARACTER(LEN=*),   INTENT(IN)   :: HGRIB     ! Grib file name
INTEGER,            INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),   INTENT(IN)   :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:), INTENT(IN)    :: PMASK     ! grib land mask
REAL, DIMENSION(:), POINTER       :: ZFIELD    ! 
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
!* Read deep soil temperature
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_T2_LAND',0,ZHOOK_HANDLE)
!
 CALL READ_GRIB_T2(HGRIB,KLUOUT,HINMODEL,PMASK,ZFIELD)
!
IF (SIZE(PMASK)==SIZE(ZFIELD)) WHERE ((PMASK(:)/=1.) .OR. ((PMASK(:)==1.) .AND.(ZFIELD(:)==9999.))) ZFIELD = XUNDEF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_T2_LAND',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_T2_LAND
!-------------------------------------------------------------------
!-------------------------------------------------------------------
SUBROUTINE PUT_LAYER_DEPTH(KLUOUT,KLEV,HROUT,KLTYPE,KLEV1,KLEV2,KNLAYERDEEP,PV4,PV,PD)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KLUOUT
INTEGER, INTENT(IN) :: KLEV
 CHARACTER(LEN=*), INTENT(IN) :: HROUT
INTEGER, INTENT(INOUT) :: KLTYPE
INTEGER, INTENT(IN) :: KLEV1
INTEGER, INTENT(IN) :: KLEV2
INTEGER, INTENT(IN) :: KNLAYERDEEP
REAL, INTENT(IN) :: PV4
REAL, INTENT(IN) :: PV
REAL, INTENT(OUT) :: PD
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:PUT_LAYER_DEPTH',0,ZHOOK_HANDLE)
!
IF (KLEV2 == -1) KLTYPE = 0
IF (KLTYPE==112) THEN
  PD = (KLEV2 - KLEV1) / 100.
ELSE
  IF (KNLAYERDEEP == 4) THEN
    PD = PV4            
  ELSE
    PD = PV
  END IF
  WRITE (KLUOUT,'(A,i1,A,f5.2,A)') 'MODE_READ_GRIB:'//TRIM(HROUT)//': | Level ',&
                                    KLEV,' height not available, assuming ',PD,' m'
END IF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:PUT_LAYER_DEPTH',1,ZHOOK_HANDLE)
END SUBROUTINE PUT_LAYER_DEPTH
!-------------------------------------------------------------------
!     #######################
      SUBROUTINE FILL_PFIELD(KLUOUT,HROUT,KNLAYERDEEP,PDIN,PFIELDIN,PMASK,PFIELDOUT,PDOUT)
!     #######################
!
!    MODIFICATIONS
!    Julien Pergaud(via J.Escobar) 10/2018 : correction for GFS undef=9999. temperature values
!
USE MODD_SURF_PAR,   ONLY : XUNDEF
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KLUOUT
 CHARACTER(LEN=*), INTENT(IN) :: HROUT
INTEGER, INTENT(IN) :: KNLAYERDEEP
REAL, DIMENSION(:), INTENT(IN) :: PDIN
REAL, DIMENSION(:,:), INTENT(IN) :: PFIELDIN
REAL, DIMENSION(:), INTENT(IN) :: PMASK
REAL, DIMENSION(:,:), POINTER :: PFIELDOUT
REAL, DIMENSION(:,:), POINTER :: PDOUT
!
 CHARACTER(LEN=20) :: FMT0
INTEGER :: JL
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:FILL_PFIELD',0,ZHOOK_HANDLE)
!--------------------------------------------------------------------------------
! 1.  Display the number of layer found
!     -----------------------
WRITE(FMT0,FMT='(A8,I1,A11)') '(A,I1,A,',KNLAYERDEEP,'(F5.2,","))'
WRITE (KLUOUT,FMT=FMT0) 'MODE_READ_GRIB:'//TRIM(HROUT)//': | ',KNLAYERDEEP,&
                        ' deep layers, heights are : ',PDIN(1:KNLAYERDEEP)
!--------------------------------------------------------------------------------
! 2.  Set temperature profile and layer thicknesses
!     -----------------------------------------------
ALLOCATE(PFIELDOUT(SIZE(PFIELDIN,1),KNLAYERDEEP))
ALLOCATE(PDOUT(SIZE(PFIELDIN,1),KNLAYERDEEP))
!
DO JL=1,KNLAYERDEEP
  PDOUT(:,JL)=SUM(PDIN(1:JL))
  PFIELDOUT(:,JL)=PFIELDIN(:,JL)
  IF (SIZE(PMASK)==SIZE(PFIELDOUT,1)) &
    WHERE ((PMASK(:)/=1.) .OR. ((PMASK(:)==1.) .AND.(PFIELDOUT(:,JL)==9999.))) PFIELDOUT(:,JL) = XUNDEF

ENDDO
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:FILL_PFIELD',1,ZHOOK_HANDLE)
END SUBROUTINE FILL_PFIELD
!-------------------------------------------------------------------
!     #######################
      SUBROUTINE READ_GRIB_TG_ECMWF(HGRIB,KLUOUT,HINMODEL,PMASK,PTG,PD)
!     #######################
!    MODIFICATIONS
!    Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 
!
IMPLICIT NONE
!
!* dummy arguments
!  ---------------
 CHARACTER(LEN=*),     INTENT(IN)    :: HGRIB     ! Grib file name
INTEGER,              INTENT(IN)    :: KLUOUT    ! logical unit of output listing
 CHARACTER(LEN=6),     INTENT(IN)    :: HINMODEL  ! Grib originating model
REAL, DIMENSION(:),   INTENT(IN)    :: PMASK     ! grib land mask
REAL, DIMENSION(:,:), POINTER       :: PTG       ! field to initialize
REAL, DIMENSION(:,:), POINTER       :: PD        ! thickness of each layer
!
!* local variables
!  ---------------
INTEGER(KIND=kindOfInt)                           :: IRET      ! return code
INTEGER                           :: ILTYPE    ! type of level (Grib code table 3)
INTEGER                           :: ILEV1     ! level definition
INTEGER                           :: ILEV2     ! level definition
INTEGER                           :: JL         ! layer loop counter
INTEGER                           :: INLAYERDEEP! number of deep moisture layers
REAL,    DIMENSION(:), POINTER    :: ZFIELD => NULL()  ! first layer temperature
REAL,  DIMENSION(:,:), ALLOCATABLE:: ZTG      ! first layer temperature
REAL, DIMENSION(:)   , ALLOCATABLE:: ZD


REAL(KIND=JPRB) :: ZHOOK_HANDLE
!--------------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TG_ECMWF',0,ZHOOK_HANDLE)
WRITE  (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_TG_ECMWF: | Reading soil temperature'
!
ALLOCATE(ZD(5))
!
! 1.  Search and read level 1 (and its depth)
!     --------------------------------------
ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
 CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
!
IF (IRET== 0) THEN
  CALL PUT_LAYER_DEPTH(KLUOUT,1,'READ_GRIB_TG_ECMWF',ILTYPE,ILEV1,ILEV2,4,0.07,0.07,ZD(1))
  ALLOCATE(ZTG(SIZE(ZFIELD),5))
  ZTG(:,1)=ZFIELD
ELSE
  CALL ABOR1_SFX('MODE_READ_GRIB: SOIL TEMPERATURE LEVEL 1 MISSING (READ_GRIB_TG_ECMWF)')
ENDIF
!
! 2.  Search and read level 4 (and its depth) This level is optionnal
!     ---------------------------------------------------------------
ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
 CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,236,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
!
IF (IRET == 0) THEN
  INLAYERDEEP = 4
  CALL PUT_LAYER_DEPTH(KLUOUT,4,'READ_GRIB_TG_ECMWF',ILTYPE,ILEV1,ILEV2,INLAYERDEEP,1.89,1.89,ZD(4))
  ZTG(:,4)=ZFIELD
ELSE
  INLAYERDEEP = 3
  ZD(4) = 0.
ENDIF