diff --git a/src/MNH/read_all_data_grib_case.f90 b/src/MNH/read_all_data_grib_case.f90
index 1dc6e59950ba47f1b176781458cfcc3d94e56c35..8f4e4835a34e939884aa517bc353b000a844383b 100644
--- a/src/MNH/read_all_data_grib_case.f90
+++ b/src/MNH/read_all_data_grib_case.f90
@@ -127,6 +127,7 @@ END MODULE MODI_READ_ALL_DATA_GRIB_CASE
!! 08/03/2018 (P.Wautelet) replace ADD_FORECAST_TO_DATE by DATETIME_CORRECTDATE
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!! Pergaud : 2018 add GFS
+!! 01/2019 (G.Delautier via Q.Rodier) for GRIB2 ARPEGE and AROME from EPYGRAM
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -323,6 +324,7 @@ REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZR_DUM
INTEGER :: IMI
TYPE(TFILEDATA),POINTER :: TZFILE
INTEGER, DIMENSION(JP_GFS) :: IP_GFS ! list of pressure levels for GFS model
+INTEGER :: IVERSION,ILEVTYPE
!---------------------------------------------------------------------------------------
IP_GFS=(/1000,975,950,925,900,850,800,750,700,650,600,550,500,450,400,350,300,&
250,200,150,100,70,50,30,20,10/)!
@@ -448,8 +450,13 @@ SELECT CASE (ICENTER)
CASE (85)
SELECT CASE (HGRID)
CASE('lambert')
- WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Aladin france model'
- IMODEL = 1
+ WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arome france model'
+ CALL GRIB_GET(IGRIB(1),'editionNumber',IVERSION,IRET_GRIB)
+ IF (IVERSION==2) THEN
+ IMODEL = 6 ! GRIB2 since summer 2018 (epygram)
+ ELSE
+ IMODEL = 1 ! GRIB1 befor summer 2018 (lfi2mv)
+ ENDIF
ALLOCATE(ZPARAM(10))
CASE('mercator')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Aladin reunion model'
@@ -458,8 +465,13 @@ SELECT CASE (ICENTER)
CASE('unknown_PLPresent','reduced_stretched_rotated_gg')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arpege model'
- IMODEL = 3
ALLOCATE(ZPARAM(10))
+ CALL GRIB_GET(IGRIB(1),'editionNumber',IVERSION,IRET_GRIB)
+ IF (IVERSION==2) THEN
+ IMODEL = 7 ! GRIB2 since summer 2018 (epygram)
+ ELSE
+ IMODEL = 3 ! GRIB1 befor summer 2018 (lfi2mv)
+ ENDIF
CASE('regular_gg')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arpege model'
@@ -488,29 +500,31 @@ END IF
WRITE (ILUOUT0,'(A)') ' | Searching orography'
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
- CALL SEARCH_FIELD(129,109,1,0,IGRIB,INUM_ZS)
+ CALL SEARCH_FIELD(IGRIB,INUM_ZS,KPARAM=129)
IF(INUM_ZS < 0) THEN
- CALL SEARCH_FIELD(129,1,0,0,IGRIB,INUM_ZS)
- IF(INUM_ZS < 0) THEN
- WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
- END IF
- END IF
+ WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
+ END IF
CASE(3,4,5) ! arpege et mocage
- CALL SEARCH_FIELD(8,-1,-1,-1,IGRIB,INUM_ZS)
+ CALL SEARCH_FIELD(IGRIB,INUM_ZS,KPARAM=8)
IF(INUM_ZS < 0) THEN
WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
ENDIF
CASE(1,2) ! aladin et aladin reunion
- CALL SEARCH_FIELD(6,-1,-1,-1,IGRIB,INUM_ZS)
+ CALL SEARCH_FIELD(IGRIB,INUM_ZS,KPARAM=6)
+ IF(INUM_ZS < 0) THEN
+ WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
+ ENDIF
+ CASE(6,7) ! arpege and arome GRIB2
+ CALL SEARCH_FIELD(IGRIB,INUM_ZS,KDIS=0,KCAT=3,KNUMBER=5)
IF(INUM_ZS < 0) THEN
WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
ENDIF
CASE(10) ! NCEP
- DO IVAR=0,222
- CALL SEARCH_FIELD(IVAR,1,0,0,IGRIB,INUM_ZS)
- IF(INUM_ZS < 0) THEN
- WRITE (ILUOUT0,'(A)')'Orography is missing'
- ENDIF
+ DO IVAR=0,222
+ CALL SEARCH_FIELD(IGRIB,INUM_ZS,KPARAM=IVAR)
+ IF(INUM_ZS < 0) THEN
+ WRITE (ILUOUT0,'(A)')'Orography is missing'
+ ENDIF
END DO
INUM_ZS=218
WRITE (ILUOUT0,*) 'lsm ',IGRIB(350)
@@ -561,11 +575,13 @@ WRITE (ILUOUT0,'(A)') ' | Searching pressure'
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
- CALL SEARCH_FIELD(152,-1,-1,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=152)
CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
- CALL SEARCH_FIELD(1,-1,-1,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=1)
+ CASE(6,7) ! NEW AROME,ARPEGE
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=3,KNUMBER=25)
CASE(10) ! NCEP
- CALL SEARCH_FIELD(134,1,0,0,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=134)
END SELECT
IF(INUM < 0) THEN
WRITE (ILUOUT0,'(A)')'Surface pressure is missing - abort'
@@ -578,7 +594,7 @@ ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
! determination des tableaux ZPS_G et ZLNPS_G
SELECT CASE (IMODEL)
- CASE(0) ! ECMWF
+ CASE(0,6,7) ! ECMWF
ALLOCATE (ZPS_G (ISIZE))
ALLOCATE (ZLNPS_G(ISIZE))
ZLNPS_G(:) = ZVALUE(1:ISIZE)
@@ -645,64 +661,55 @@ DEALLOCATE (ZLNPS_G)
!
WRITE (ILUOUT0,'(A)') ' | Reading T and Q fields'
!
-SELECT CASE (IMODEL)
- CASE(0) ! ECMWF
+IF (IMODEL/=10) THEN
+ SELECT CASE (IMODEL)
+ CASE(0) ! ECMWF
ISTARTLEVEL=1
IT=130
IQ=133
- CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
- IF(INUM < 0) THEN
- ISTARTLEVEL=0
- CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
- ENDIF
- IF(INUM < 0) THEN
- WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
- CALL ABORT
- STOP
- ENDIF
- CALL SEARCH_FIELD(IQ,109,ISTARTLEVEL,-1,IGRIB,INUM)
- IF(INUM < 0) THEN
- WRITE (ILUOUT0,'(A)')'Atmospheric specific humidity is missing - abort'
- CALL ABORT
- STOP
- ENDIF
- CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
+ CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
IT=11
IQ=51
ISTARTLEVEL=1
- CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
- IF(INUM < 0) THEN
- ISTARTLEVEL=0
- CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
- ENDIF
- IF(INUM < 0) THEN
- WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
- CALL ABORT
- STOP
- ENDIF
- CALL SEARCH_FIELD(IQ,109,ISTARTLEVEL,-1,IGRIB,INUM)
- IF(INUM < 0) THEN
- WRITE (ILUOUT0,'(A)')'Atmospheric specific humidity is missing - abort'
- CALL ABORT
- STOP
- ENDIF
- CASE(10) ! NCEP
- ISTARTLEVEL=10
- IT=130
- IQ=157
- CALL SEARCH_FIELD(IT,100,ISTARTLEVEL,-1,IGRIB,INUM)
- IF(INUM < 0) THEN
- WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
- CALL ABORT
- STOP
- ENDIF
- CALL SEARCH_FIELD(IQ,100,ISTARTLEVEL,-1,IGRIB,INUM)
- IF(INUM < 0) THEN
- WRITE (ILUOUT0,'(A)')'Atmospheric relative humidity is missing - abort'
- CALL ABORT
- STOP
- ENDIF
-END SELECT
+ CASE(6,7) !GRIB2 AROME AND ARPEGE
+ IT=130
+ IQ=133
+ ISTARTLEVEL=1
+ END SELECT
+
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ISTARTLEVEL)
+ IF(INUM < 0) THEN
+ ISTARTLEVEL=0
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ISTARTLEVEL)
+ ENDIF
+ IF(INUM < 0) THEN
+ WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
+ CALL ABORT
+ STOP
+ ENDIF
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ISTARTLEVEL)
+ IF(INUM < 0) THEN
+ WRITE (ILUOUT0,'(A)')'Atmospheric specific humidity is missing - abort'
+ CALL ABORT
+ STOP
+ ENDIF
+ELSE ! NCEP
+ ISTARTLEVEL=10
+ IT=130
+ IQ=157
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ISTARTLEVEL)
+ IF(INUM < 0) THEN
+ WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
+ CALL ABORT
+ STOP
+ ENDIF
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ISTARTLEVEL)
+ IF(INUM < 0) THEN
+ WRITE (ILUOUT0,'(A)')'Atmospheric relative humidity is missing - abort'
+ CALL ABORT
+ STOP
+ ENDIF
+ENDIF
!
IF (IMODEL/=10) THEN ! others than NCEP
CALL GRIB_GET(IGRIB(INUM),'NV',INLEVEL)
@@ -718,14 +725,14 @@ ALLOCATE (ZQ_G(ISIZE,INLEVEL))
IF (IMODEL/=10) THEN ! others than NCEP
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1-1+ISTARTLEVEL
- CALL SEARCH_FIELD(IQ,109,ILEV1,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ILEV1)
IF (INUM< 0) THEN
!callabortstop
WRITE(YMSG,*) 'atmospheric humidity level ',JLOOP1,' is missing'
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
CALL GRIB_GET(IGRIB(INUM),'values',ZQ_G(:,INLEVEL-JLOOP1+1))
- CALL SEARCH_FIELD(IT,109,ILEV1,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ILEV1)
IF (INUM< 0) THEN
!callabortstop
WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
@@ -737,7 +744,7 @@ IF (IMODEL/=10) THEN ! others than NCEP
ELSE ! NCEP
DO JLOOP1=1, INLEVEL
ILEV1 = IP_GFS(JLOOP1)
- CALL SEARCH_FIELD(IQ,100,ILEV1,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ILEV1)
IF (INUM< 0) THEN
!callabortstop
WRITE(YMSG,*) 'atmospheric humidity level ',JLOOP1,' is missing'
@@ -745,7 +752,7 @@ ELSE ! NCEP
END IF
CALL GRIB_GET(IGRIB(INUM),'values',ZQ_G(:,JLOOP1),IRET_GRIB)
WRITE (ILUOUT0,*) 'Q ',ILEV1,IRET_GRIB
- CALL SEARCH_FIELD(IT,100,ILEV1,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ILEV1)
IF (INUM< 0) THEN
!callabortstop
WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
@@ -790,7 +797,7 @@ IF (IMODEL/=10) THEN ! others than NCEP
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','there is no PV value in this message')
ENDIF
SELECT CASE (IMODEL)
- CASE (0,3,4)
+ CASE (0,3,4,6,7)
DO JLOOP1 = 1, INLEVEL
XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+JLOOP1) / XP00_LS
XB_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
@@ -860,10 +867,10 @@ IF (IMODEL==1) THEN
! search cloud_water in Arome case (forecast)
ISTARTLEVEL = 1
IPAR=246
- CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
IF (INUM < 0) THEN
ISTARTLEVEL = 0
- CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
END IF
IF (INUM > 0) THEN
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arome model (forecast)'
@@ -873,10 +880,33 @@ IF (IMODEL==1) THEN
! search also turbulent kinetic energy
ISTARTLEVEL = 1
IPAR=251
- CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
IF (INUM < 0) THEN
ISTARTLEVEL = 0
- CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
+ END IF
+ IF (INUM > 0) CTURB='TKEL'
+END IF
+
+IF (IMODEL==6) THEN ! GRIB2 AROME
+! search cloud_water in Arome case (forecast)
+ ISTARTLEVEL = 1
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=6,KNUMBER=6,KLEV1=ISTARTLEVEL)
+ IF (INUM < 0) THEN
+ ISTARTLEVEL = 0
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=6,KNUMBER=6,KLEV1=ISTARTLEVEL)
+ END IF
+ IF (INUM > 0) THEN
+ WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arome model (forecast)'
+ LCPL_AROME=.TRUE.
+ NRR=6
+ END IF
+ ! search also turbulent kinetic energy
+ ISTARTLEVEL = 1
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=19,KNUMBER=11,KLEV1=ISTARTLEVEL)
+ IF (INUM < 0) THEN
+ ISTARTLEVEL = 0
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=19,KNUMBER=11,KLEV1=ISTARTLEVEL)
END IF
IF (INUM > 0) CTURB='TKEL'
END IF
@@ -1034,12 +1064,35 @@ DEALLOCATE (ZEXNM_LS)
!* 2.5.8 Read the other specific ratios (if Arome model)
!
IF (NRR >1) THEN
- WRITE (ILUOUT0,'(A)') ' | Reading Q fields (except humidity)'
- DO JLOOP2=1,NRR-1
- IPAR=246+JLOOP2-1
+ IF (IMODEL==1) THEN
+ WRITE (ILUOUT0,'(A)') ' | Reading Q fields (except humidity)'
+ DO JLOOP2=1,NRR-1
+ IPAR=246+JLOOP2-1
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = JLOOP1-1+ISTARTLEVEL
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ILEV1)
+
+ IF (INUM < 0) THEN
+ !callabortstop
+ WRITE(YMSG,*) 'Specific ratio ',IPAR,' at level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
+ ALLOCATE(ZVALUE(ISIZE))
+ CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
+ ALLOCATE(ZOUT(INO))
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+ CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,1+JLOOP2))
+ DEALLOCATE (ZVALUE)
+ DEALLOCATE (ZOUT)
+ END DO
+ END DO
+ ELSE ! GRIB2 AROME IMODEL =6
+ WRITE (ILUOUT0,'(A)') ' | Reading Q fields (except humidity)'
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1-1+ISTARTLEVEL
- CALL SEARCH_FIELD(IPAR,109,ILEV1,-1,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=6,KNUMBER=6,KLEV1=ILEV1)
IF (INUM < 0) THEN
!callabortstop
@@ -1052,19 +1105,102 @@ IF (NRR >1) THEN
ALLOCATE(ZOUT(INO))
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
- CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,1+JLOOP2))
+ CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,2))
DEALLOCATE (ZVALUE)
DEALLOCATE (ZOUT)
END DO
- END DO
+
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = JLOOP1-1+ISTARTLEVEL
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=1,KNUMBER=85,KLEV1=ILEV1)
+
+ IF (INUM < 0) THEN
+ !callabortstop
+ WRITE(YMSG,*) 'Specific ratio for rain at level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
+ ALLOCATE(ZVALUE(ISIZE))
+ CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
+ ALLOCATE(ZOUT(INO))
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+ CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,3))
+ DEALLOCATE (ZVALUE)
+ DEALLOCATE (ZOUT)
+ END DO
+
+
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = JLOOP1-1+ISTARTLEVEL
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=1,KNUMBER=82,KLEV1=ILEV1)
+ IF (INUM < 0) THEN
+ !callabortstop
+ WRITE(YMSG,*) 'Specific ratio for ICE at level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
+ ALLOCATE(ZVALUE(ISIZE))
+ CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
+ ALLOCATE(ZOUT(INO))
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+ CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,4))
+ DEALLOCATE (ZVALUE)
+ DEALLOCATE (ZOUT)
+ END DO
+
+
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = JLOOP1-1+ISTARTLEVEL
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=1,KNUMBER=86,KLEV1=ILEV1)
+ IF (INUM < 0) THEN
+ !callabortstop
+ WRITE(YMSG,*) 'Specific ratio ',IPAR,' at level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
+ ALLOCATE(ZVALUE(ISIZE))
+ CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
+ ALLOCATE(ZOUT(INO))
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+ CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,5))
+ DEALLOCATE (ZVALUE)
+ DEALLOCATE (ZOUT)
+ END DO
+
+
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = JLOOP1-1+ISTARTLEVEL
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=1,KNUMBER=32,KLEV1=ILEV1)
+ IF (INUM < 0) THEN
+ !callabortstop
+ WRITE(YMSG,*) 'Specific ratio ',IPAR,' at level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
+ ALLOCATE(ZVALUE(ISIZE))
+ CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
+ ALLOCATE(ZOUT(INO))
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+ CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,6))
+ DEALLOCATE (ZVALUE)
+ DEALLOCATE (ZOUT)
+ END DO
+ END IF
END IF
!
IF (CTURB=='TKEL') THEN
WRITE (ILUOUT0,'(A)') ' | Reading TKE field'
- IPAR=251
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1-1+ISTARTLEVEL
- CALL SEARCH_FIELD(IPAR,109,ILEV1,-1,IGRIB,INUM)
+ IF (IMODEL==1) THEN
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=251,KLEV1=ILEV1)
+ ELSE ! case 6 new arome
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=19,KNUMBER=11,KLEV1=ILEV1)
+ END IF
IF (INUM < 0) THEN
!callabortstop
WRITE(YMSG,*) 'TKE at level ',JLOOP1,' is missing'
@@ -1109,9 +1245,7 @@ IF (IMODEL==5) THEN
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','Mocage model: Bad input argument in read_all_data_grib_case')
END IF
XSV_LS(:,:,:,:) = 0.
- ITYP=109
ILEV1=-1
- ILEV2=-1
!
WRITE (ILUOUT0,'(A,A4,A)') ' | Reading Mocage species (ppp) from ',HFILE,' file'
!
@@ -1166,7 +1300,7 @@ IF (IMODEL==5) THEN
IF (INACT .NE. 0) THEN
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1
- CALL SEARCH_FIELD(INUMGRIB(JN),ITYP,ILEV1,ILEV2,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=INUMGRIB(JN),KLEV1=ILEV1)
IF (INUM < 0) THEN
!callabortstop
WRITE(YMSG,*) 'Atmospheric ',INUMGRIB(JN),' grib chemical species level ',JLOOP1,' is missing'
@@ -1233,32 +1367,29 @@ END IF
!
IF (HFILE(1:3)=='ATM') THEN
IF (IMODEL/=10) THEN ! others than NCEP
- ITYP = 109
ISTARTLEVEL = 1
ELSE
- ITYP = 100
ISTARTLEVEL = 10
END IF
-ILEV2 = -1
ALLOCATE (XU_LS(IIU,IJU,INLEVEL))
ALLOCATE (XV_LS(IIU,IJU,INLEVEL))
ALLOCATE (ZTU_LS(INO))
ALLOCATE (ZTV_LS(INO))
!
SELECT CASE (IMODEL)
- CASE (0)
+ CASE (0,6,7)
IPAR = 131
- CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
IF (INUM< 0) THEN
ISTARTLEVEL = 0
- CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
END IF
CASE (1,2,3)
IPAR = 33
- CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
IF (INUM < 0) THEN
ISTARTLEVEL = 0
- CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
END IF
CASE (10)
IPAR = 131
@@ -1272,7 +1403,7 @@ DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
ILEV1 = IP_GFS(JLOOP1)
END IF
! read component u
- CALL SEARCH_FIELD(IPAR,ITYP,ILEV1,ILEV2,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ILEV1)
IF (INUM < 0) THEN
!callabortstop
WRITE(YMSG,*) 'wind vector component "u" at level ',JLOOP1,' is missing'
@@ -1281,7 +1412,7 @@ DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
- IF (IMODEL==3) THEN
+ IF (IMODEL==3.OR.(IMODEL==7)) THEN
ALLOCATE(ZTU0_LS(INI))
ZTU0_LS(:) = ZVALUE(:)
ELSE
@@ -1304,7 +1435,7 @@ DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
ELSE
ILEV1 = IP_GFS(JLOOP1)
END IF
- CALL SEARCH_FIELD(IPAR+1,ITYP,ILEV1,ILEV2,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR+1,KLEV1=ILEV1)
IF (INUM < 0) THEN
!callabortstop
WRITE(YMSG,*) 'wind vector component "v" at level ',JLOOP1,' is missing'
@@ -1313,7 +1444,7 @@ DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
- IF ((IMODEL==3)) THEN
+ IF ((IMODEL==3).OR.(IMODEL==7)) THEN
CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
ALLOCATE(IINLO(INJ))
CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
@@ -1332,8 +1463,8 @@ DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
DEALLOCATE(ZOUT)
END IF
DEALLOCATE (ZVALUE)
- ! interpolations for arpege grid
- IF ((IMODEL==3)) THEN
+ ! interpolations for arpege grid
+ IF ((IMODEL==3).OR.(IMODEL==7)) THEN
! Comes back to real winds instead of stretched winds
! (but still with components according to Arpege grid axes)
ZLATPOLE = ZPARAM(7) * XPI/180.
@@ -1413,7 +1544,7 @@ DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
ZTU_LS(JLOOP4)*SIN(ZALPHA) + ZTV_LS(JLOOP4)*COS(ZALPHA)
ENDDO
ENDDO
- IF ((IMODEL==3)) THEN ! deallocation of Arpege arrays
+ IF ((IMODEL==3).OR.(IMODEL==7)) THEN ! deallocation of Arpege arrays
DEALLOCATE (ZTU0_LS)
DEALLOCATE (ZTV0_LS)
END IF
@@ -1618,7 +1749,7 @@ IF (ODUMMY_REAL) THEN
!
DO JI = 1, IMOC
WRITE(ILUOUT0,'(A,4(I3,1X))') CDUMMY_2D(JI),INUMGRIB(JI),INUMLEV(JI),INUMLEV1(JI),INUMLEV2(JI)
- CALL SEARCH_FIELD(IPAR,ITYP,ILEV1,ILEV2,IGRIB,INUM)
+ CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ILEV1)
IF (INUM < 0) THEN
WRITE (ILUOUT0,'(A,I3,A,I2,A)') ' -> 2D field ',INUMGRIB(JI),' is missing - abort'
!callabortstop
@@ -1715,7 +1846,7 @@ END SUBROUTINE ARRAY_1D_TO_2D
!---------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------
!#################################################################################
-SUBROUTINE SEARCH_FIELD(KPARAM,KLTYPE,KLEV1,KLEV2,KGRIB,KNUM)
+SUBROUTINE SEARCH_FIELD(KGRIB,KNUM,KPARAM,KDIS,KCAT,KNUMBER,KLEV1)
!#################################################################################
! search the grib message corresponding to KPARAM,KLTYPE,KLEV1,KLEV2 in all
! the KGIRB messages
@@ -1728,40 +1859,37 @@ USE MODE_IO_ll
IMPLICIT NONE
!
!
-INTEGER,INTENT(IN) :: KPARAM ! Parameter to read
-INTEGER,INTENT(IN) :: KLTYPE ! Level type
-INTEGER,INTENT(IN) :: KLEV1 ! Level parameter 1
-INTEGER,INTENT(IN) :: KLEV2 ! Level parameter 2
INTEGER(KIND=kindOfInt),DIMENSION(:),INTENT(IN) :: KGRIB ! number of grib messages
INTEGER,INTENT(OUT) :: KNUM ! number of the message researched
+INTEGER,INTENT(IN),OPTIONAL :: KPARAM ! INdicator of parameter/paramId
+INTEGER,INTENT(IN),OPTIONAL :: KDIS ! Discipline (GRIB2)
+INTEGER,INTENT(IN),OPTIONAL :: KCAT ! Catégorie (GRIB2)
+INTEGER,INTENT(IN),OPTIONAL :: KNUMBER ! parameterNumber (GRIB2)
+INTEGER,INTENT(IN),OPTIONAL :: KLEV1 ! Level
!
! Declaration of local variables
!
INTEGER :: IFOUND ! Number of correct parameters
+INTEGER :: ISEARCH ! Number of correct parameters to find
INTEGER :: IRET ! error code
-INTEGER :: IPARAM ! Parameter read
+INTEGER :: IPARAM,IDIS,ICAT,INUMBER
INTEGER :: ILEV1 ! Level parameter 1
-INTEGER :: ILEV2 ! Level parameter 2
INTEGER :: JLOOP ! Dummy counter
INTEGER :: IVERSION
-CHARACTER(LEN=24) :: YLTYPELU
-CHARACTER(LEN=20) :: YLTYPE
-!
! Variables used to display messages
INTEGER :: ILUOUT0 ! Logical unit number of the listing
!
ILUOUT0 = TLUOUT0%NLU
!
-SELECT CASE (KLTYPE)
-CASE(100)
- YLTYPE='isobaricInhPa'
-CASE(109)
- YLTYPE='hybrid'
-CASE(1)
- YLTYPE='surface'
-CASE DEFAULT
- YLTYPE='unknown'
-END SELECT
+ISEARCH=0
+IF (PRESENT(KPARAM)) ISEARCH=ISEARCH+1
+IF (PRESENT(KDIS)) ISEARCH=ISEARCH+1
+IF (PRESENT(KCAT)) ISEARCH=ISEARCH+1
+IF (PRESENT(KNUMBER)) ISEARCH=ISEARCH+1
+IF (PRESENT(KLEV1)) ISEARCH=ISEARCH+1
+
+
+
DO JLOOP=1,SIZE(KGRIB)
IFOUND = 0
!
@@ -1773,50 +1901,89 @@ DO JLOOP=1,SIZE(KGRIB)
WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
CYCLE
ENDIF
- IF (IVERSION == 2) THEN
- CALL GRIB_GET(KGRIB(JLOOP),'paramId',IPARAM,IRET_GRIB)
- ELSE
- CALL GRIB_GET(KGRIB(JLOOP),'indicatorOfParameter',IPARAM,IRET_GRIB)
- ENDIF
- IF (IRET_GRIB > 0) THEN
- WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
- CYCLE
- ELSE IF (IRET_GRIB == -6) THEN
- WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
- CYCLE
+ IF (PRESENT(KPARAM)) THEN
+ IF (IVERSION == 2) THEN
+ CALL GRIB_GET(KGRIB(JLOOP),'paramId',IPARAM,IRET_GRIB)
+ ELSE
+ CALL GRIB_GET(KGRIB(JLOOP),'indicatorOfParameter',IPARAM,IRET_GRIB)
+ ENDIF
+ IF (IRET_GRIB > 0) THEN
+ WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
+ CYCLE
+ ELSE IF (IRET_GRIB == -6) THEN
+ WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
+ CYCLE
+ ENDIF
+ IF (IPARAM==KPARAM) THEN
+ IFOUND = IFOUND + 1
+ ELSE
+ CYCLE
+ ENDIF
ENDIF
!
- CALL GRIB_GET(KGRIB(JLOOP),'typeOfLevel',YLTYPELU,IRET_GRIB)
- IF (IRET_GRIB > 0) THEN
- WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
- CYCLE
- ELSE IF (IRET_GRIB == -6) THEN
- WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
- CYCLE
+ IF (PRESENT(KDIS)) THEN
+ CALL GRIB_GET(KGRIB(JLOOP),'discipline',IDIS,IRET_GRIB)
+ IF (IRET_GRIB > 0) THEN
+ WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
+ CYCLE
+ ELSE IF (IRET_GRIB == -6) THEN
+ WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
+ CYCLE
+ ENDIF
+ IF (IDIS==KDIS) THEN
+ IFOUND = IFOUND + 1
+ ELSE
+ CYCLE
+ ENDIF
+ ENDIF
+ IF (PRESENT(KCAT)) THEN
+ CALL GRIB_GET(KGRIB(JLOOP),'parameterCategory',ICAT,IRET_GRIB)
+ IF (IRET_GRIB > 0) THEN
+ WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
+ CYCLE
+ ELSE IF (IRET_GRIB == -6) THEN
+ WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
+ CYCLE
+ ENDIF
+ IF (ICAT==KCAT) THEN
+ IFOUND = IFOUND + 1
+ ELSE
+ CYCLE
+ ENDIF
+ ENDIF
+ IF (PRESENT(KNUMBER)) THEN
+ CALL GRIB_GET(KGRIB(JLOOP),'parameterNumber',INUMBER,IRET_GRIB)
+ IF (IRET_GRIB > 0) THEN
+ WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
+ CYCLE
+ ELSE IF (IRET_GRIB == -6) THEN
+ WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
+ CYCLE
+ ENDIF
+ IF (INUMBER==KNUMBER) THEN
+ IFOUND = IFOUND + 1
+ ELSE
+ CYCLE
+ ENDIF
ENDIF
!
- CALL GRIB_GET(KGRIB(JLOOP),'topLevel',ILEV1,IRET_GRIB)
- IF (IRET_GRIB > 0) THEN
- WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
- CYCLE
- ELSE IF (IRET_GRIB == -6) THEN
- WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
- CYCLE
+ IF(PRESENT(KLEV1)) THEN
+ CALL GRIB_GET(KGRIB(JLOOP),'topLevel',ILEV1,IRET_GRIB)
+ IF (IRET_GRIB > 0) THEN
+ WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
+ CYCLE
+ ELSE IF (IRET_GRIB == -6) THEN
+ WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
+ CYCLE
+ ENDIF
+ IF (ILEV1==KLEV1) THEN
+ IFOUND = IFOUND + 1
+ ELSE
+ CYCLE
+ ENDIF
ENDIF
!
- CALL GRIB_GET(KGRIB(JLOOP),'bottomLevel',ILEV2,IRET_GRIB)
- IF (IRET_GRIB > 0) THEN
- WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
- CYCLE
- ELSE IF (IRET_GRIB == -6) THEN
- WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
- CYCLE
- ENDIF
- IF (IPARAM==KPARAM) IFOUND = 1
- IF (YLTYPELU==YLTYPE .OR. -1==KLTYPE) IFOUND = IFOUND + 1
- IF (ILEV1== KLEV1 .OR. -1== KLEV1) IFOUND = IFOUND + 1
- IF (ILEV2== KLEV2 .OR. -1== KLEV2) IFOUND = IFOUND + 1
- IF (IFOUND == 4) THEN
+ IF (IFOUND == ISEARCH) THEN
KNUM=JLOOP
EXIT
ELSE ! field not found
@@ -1826,7 +1993,6 @@ DO JLOOP=1,SIZE(KGRIB)
END DO
END SUBROUTINE SEARCH_FIELD
-
!#################################################################################
SUBROUTINE COORDINATE_CONVERSION(KMODEL,KGRIB,KNOLON,KNOLARG,&
PLONOUT,PLATOUT,PLXOUT,PLYOUT,KNI,PPARAM,KINLO)
@@ -1961,7 +2127,7 @@ CASE(0,5) ! CEP/MOCAGE
PLYOUT=PLATOUT
ENDIF
!
-CASE(1) ! ALADIN
+CASE(1,6) ! ALADIN
!
CALL GRIB_GET(KGRIB,'Nj',IINLA,IRET_GRIB)
CALL GRIB_GET(KGRIB,'Ni',INLO_GRIB(1),IRET_GRIB)
@@ -2084,8 +2250,14 @@ CASE(2) ! ALADIN REUNION
DEALLOCATE (ZXM)
END IF
!
-CASE(3,4) ! ARPEGE
+CASE(3,4,7) ! ARPEGE
!
+!print*,"=========COORDINATE CONVERSION CASE ARPEGE ============="
+! PROBLEME AVEC LES GRIB d'EPYGRAM
+! dans longitudeOfLastGridPointInDegrees on la la longitude du dernier point du
+! tableau (donc au pole sud)
+! dans les GRIB1 ont avait la valeur max du tableau des longitude (donc Ã
+! l'equateur)
CALL GRIB_GET(KGRIB,'latitudeOfFirstGridPointInDegrees',ZILA1)
CALL GRIB_GET(KGRIB,'longitudeOfFirstGridPointInDegrees',ZILO1)
CALL GRIB_GET(KGRIB,'latitudeOfLastGridPointInDegrees',ZILA2)
@@ -2111,6 +2283,7 @@ CASE(3,4) ! ARPEGE
DO JLOOP1=1 ,IINLA
KNI = KNI + INLO_GRIB(JLOOP1)
ENDDO
+ ZILO2=360.-360./(MAXVAL(INLO_GRIB))
GREADY= (PPARAM(1)==INLO_GRIB(1) .AND.&
PPARAM(3)==ZILA1 .AND. PPARAM(4)==ZILO1 .AND.&
PPARAM(5)==ZILA2 .AND. PPARAM(6)==ZILO2 .AND.&
diff --git a/src/SURFEX/mode_read_grib.F90 b/src/SURFEX/mode_read_grib.F90
index 501ae26da178a98e7ee04f3330aaaaad92f24211..ddd50d15ea3a3397107ac71dea73282a5edd4b19 100644
--- a/src/SURFEX/mode_read_grib.F90
+++ b/src/SURFEX/mode_read_grib.F90
@@ -16,14 +16,17 @@ USE PARKIND1 ,ONLY : JPRB
CONTAINS
!-------------------------------------------------------------------
! ####################
- SUBROUTINE MAKE_GRIB_INDEX(HGRIB)
+ 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
@@ -54,7 +57,13 @@ ENDIF
IF (NGRIB_VERSION==1) THEN
CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'indicatorOfParameter',IRET)
ELSEIF (NGRIB_VERSION==2) THEN
- CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'paramId',IRET)
+ 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")
!
@@ -89,6 +98,8 @@ 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
!
@@ -120,13 +131,20 @@ 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)
@@ -147,6 +165,7 @@ DO WHILE (IRET /= GRIB_END_OF_INDEX .AND. KFOUND/=3)
ENDIF
!
+print*,KFOUND
IF (IRET.EQ.0) THEN
!
KFOUND = KFOUND + 1
@@ -161,6 +180,7 @@ DO WHILE (IRET /= GRIB_END_OF_INDEX .AND. KFOUND/=3)
!
IF (IRET.EQ.0) KFOUND = KFOUND + 1
!
+print*,KFOUND
ENDIF
!
ENDIF
@@ -270,15 +290,18 @@ END SUBROUTINE TEST_IRET_FLOAT
END SUBROUTINE GET_GRIB_MESSAGE
!-------------------------------------------------------------------
! ####################
- SUBROUTINE READ_GRIB(HGRIB,KLUOUT,KPARAM,KRET,PFIELD,KLTYPE,KLEV1,KLEV2,KPARAM2, &
- HTYPELEVEL,PLEV1,PLEV2)
+ 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
@@ -287,11 +310,12 @@ 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
+INTEGER :: ILTYPE, ILEV1, ILEV2,ILTYPE2
INTEGER(KIND=kindOfInt) :: IGRIB
INTEGER :: ISIZE, IFOUND
REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -310,7 +334,7 @@ IF (PRESENT(PLEV1)) ZLEV1=PLEV1
ZLEV2=-2.0
IF (PRESENT(PLEV2)) ZLEV2=PLEV2
!
- CALL MAKE_GRIB_INDEX(HGRIB)
+ CALL MAKE_GRIB_INDEX(HGRIB,HINMODEL)
!
IFOUND=0
KRET=0
@@ -318,30 +342,42 @@ KRET=0
IF (NGRIB_VERSION == 1) THEN
CALL GRIB_INDEX_SELECT(NIDX,'indicatorOfParameter',KPARAM,KRET)
ELSEIF (NGRIB_VERSION == 2) THEN
- CALL GRIB_INDEX_SELECT(NIDX,'paramId',KPARAM,KRET)
+ 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
+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
+WRITE (KLUOUT,*) 'READ_GRIB GRIB_NEW_FROM_INDEX ',KPARAM,IGRIB,KRET,IFOUND
IF (PRESENT(KPARAM2)) THEN
IF (IFOUND/=3) THEN
- CALL GRIB_INDEX_SELECT(NIDX,'indicatorOfParameter',KPARAM2,KRET)
+ 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
@@ -349,7 +385,7 @@ IF (PRESENT(KPARAM2)) THEN
ENDIF
ENDIF
!
-
+print*,"IFOUND=",IFOUND
IF (IFOUND==3) THEN
!
@@ -383,6 +419,10 @@ 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
!
@@ -392,26 +432,35 @@ CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
REAL, DIMENSION(:), POINTER :: PMASK ! Land mask
!
INTEGER(KIND=kindOfInt) :: IRET ! return code
-INTEGER :: ILTYPE ! leveltype
+INTEGER :: ILTYPE,ILTYPE2 ! leveltype
INTEGER :: ILEV ! level
-INTEGER :: INUM_ZS,ISIZE,ICOUNT,JLOOP,IPARAM,IGRIB
+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,KLUOUT,172,IRET,PMASK)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,172,IRET,PMASK)
CASE ('NCEP ')
- CALL READ_GRIB(HGRIB,KLUOUT,172,IRET,PMASK)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,172,IRET,PMASK)
CASE ('ARPEGE','ALADIN','MOCAGE')
- CALL READ_GRIB(HGRIB,KLUOUT,81,IRET,PMASK)
+ 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,KLUOUT,81,IRET,PMASK,KLTYPE=ILTYPE,KLEV1=ILEV)
+ 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
@@ -432,8 +481,11 @@ 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
!
@@ -455,13 +507,17 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_ZS: | Reading orography from ',HI
!
SELECT CASE (HINMODEL)
CASE ('ECMWF ')
- CALL READ_GRIB(HGRIB,KLUOUT,129,IRET,PZS)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,129,IRET,PZS)
CASE ('NCEP ')
- CALL READ_GRIB(HGRIB,KLUOUT,228002,IRET,PZS)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228002,IRET,PZS)
CASE ('ARPEGE','MOCAGE')
- CALL READ_GRIB(HGRIB,KLUOUT,8,IRET,PZS)
+ 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,KLUOUT,6,IRET,PZS)
+ 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
@@ -527,7 +583,11 @@ 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
@@ -538,6 +598,7 @@ 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
@@ -548,18 +609,24 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_T: | Reading surface temperature'
!
SELECT CASE (HINMODEL)
CASE ('ECMWF ')
- CALL READ_GRIB(HGRIB,KLUOUT,139,IRET,PT)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,139,IRET,PT)
CASE ('ARPEGE','ALADIN','MOCAGE')
ILEV=0
- ILTYPE=111
- CALL READ_GRIB(HGRIB,KLUOUT,11,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV)
+ 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,KLUOUT,11,IRET,PT,KLTYPE=ILTYPE)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,PT,KLTYPE=ILTYPE)
IF (IRET /= 0) THEN
ILTYPE=105
- CALL READ_GRIB(HGRIB,KLUOUT,11,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV)
ENDIF
END IF
@@ -567,13 +634,13 @@ SELECT CASE (HINMODEL)
WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_T: | Reading surface temperature tile 4'
ILTYPE=105
ILEV=904
- CALL READ_GRIB(HGRIB,KLUOUT,11,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV)
+ 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,KLUOUT,130,IRET,PT,KLTYPE=ILTYPE,KLEV1=ILEV,KLEV2=ILEV,HTYPELEVEL=YTYPELEVEL)
+ 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)
@@ -614,6 +681,8 @@ 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
!
@@ -632,7 +701,7 @@ IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_SST',0,ZHOOK_HANDLE)
!
SELECT CASE (HINMODEL)
CASE ('ECMWF ')
- CALL READ_GRIB(HGRIB,KLUOUT,34,IRET,PSST)
+ 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)
@@ -651,6 +720,7 @@ END SUBROUTINE READ_GRIB_SST
!
! 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
!
@@ -669,7 +739,7 @@ IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TSWATER',0,ZHOOK_HANDLE)
!
SELECT CASE (HINMODEL)
CASE ('ECMWF ')
- CALL READ_GRIB(HGRIB,KLUOUT,3080,IRET,PTS)
+ 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)
@@ -678,6 +748,7 @@ SELECT CASE (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
@@ -685,8 +756,11 @@ 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
!
@@ -698,6 +772,7 @@ 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
@@ -708,24 +783,30 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_T2: | Reading deep soil temperatu
!
SELECT CASE (HINMODEL)
CASE ('ECMWF ')
- CALL READ_GRIB(HGRIB,KLUOUT,170,IRET,PT2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,170,IRET,PT2)
CASE ('ARPEGE','ALADIN','MOCAGE')
- ILTYPE=111
- CALL READ_GRIB(HGRIB,KLUOUT,11,IRET,PT2,KLTYPE=ILTYPE)
+ 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,KLUOUT,11,IRET,PT2,KLTYPE=ILTYPE)
+ 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,KLUOUT,228139,IRET,PT2,KLTYPE=ILTYPE,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
+ 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,KLUOUT,11,IRET,PT2,KLTYPE=ILTYPE,KLEV1=ILEV)
+ 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
@@ -844,15 +925,17 @@ DO JL=1,KNLAYERDEEP
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
!
@@ -890,7 +973,7 @@ ALLOCATE(ZD(5))
ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
- CALL READ_GRIB(HGRIB,KLUOUT,139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ 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))
@@ -905,7 +988,7 @@ ENDIF
ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
- CALL READ_GRIB(HGRIB,KLUOUT,236,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,236,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
!
IF (IRET == 0) THEN
INLAYERDEEP = 4
@@ -921,7 +1004,7 @@ ENDIF
ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
- CALL READ_GRIB(HGRIB,KLUOUT,183,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,183,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
!
IF (IRET == 0) THEN
CALL PUT_LAYER_DEPTH(KLUOUT,3,'READ_GRIB_TG_ECMWF',ILTYPE,ILEV1,ILEV2,INLAYERDEEP,0.72,0.42,ZD(3))
@@ -936,7 +1019,7 @@ ENDIF
ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
- CALL READ_GRIB(HGRIB,KLUOUT,170,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,170,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
!
IF (IRET== 0) THEN
CALL PUT_LAYER_DEPTH(KLUOUT,2,'READ_GRIB_TG_ECMWF',ILTYPE,ILEV1,ILEV2,INLAYERDEEP,0.21,0.42,ZD(2))
@@ -993,6 +1076,9 @@ END SUBROUTINE READ_GRIB_TG_ECMWF
! given. If they are present, they are used, if not the default value is taken and
! a warning message is issued.
!
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
+!
IMPLICIT NONE
!
!* dummy arguments
@@ -1028,7 +1114,7 @@ ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
IPAR=39
- CALL READ_GRIB(HGRIB,KLUOUT,140,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,KPARAM2=IPAR)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,140,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,KPARAM2=IPAR)
!
IF (IRET == 0) THEN
CALL PUT_LAYER_DEPTH(KLUOUT,1,'READ_GRIB_WG_ECMWF_1',ILTYPE,ILEV1,ILEV2,4,0.07,0.07,ZD(1))
@@ -1046,7 +1132,7 @@ ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
IPAR=42
- CALL READ_GRIB(HGRIB,KLUOUT,237,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,KPARAM2=IPAR)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,237,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,KPARAM2=IPAR)
!
IF (IRET == 0) THEN
INLAYERDEEP = 4
@@ -1065,7 +1151,7 @@ ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
IPAR=41
- CALL READ_GRIB(HGRIB,KLUOUT,184,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,KPARAM2=IPAR)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,184,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,KPARAM2=IPAR)
!
IF (IRET == 0) THEN
CALL PUT_LAYER_DEPTH(KLUOUT,3,'READ_GRIB_WG_ECMWF_1',ILTYPE,ILEV1,ILEV2,INLAYERDEEP,0.72,0.42,ZD(3))
@@ -1083,7 +1169,7 @@ ILTYPE= -1
ILEV1 = -1
ILEV2 = -1
IPAR=40
- CALL READ_GRIB(HGRIB,KLUOUT,171,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,KPARAM2=IPAR)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,171,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,KPARAM2=IPAR)
!
IF (IRET == 0) THEN
CALL PUT_LAYER_DEPTH(KLUOUT,2,'READ_GRIB_WG_ECMWF_1',ILTYPE,ILEV1,ILEV2,INLAYERDEEP,0.21,0.42,ZD(2))
@@ -1228,6 +1314,8 @@ END SUBROUTINE HARMONIZE_GRIB_WG_WGI_ECMWF
! #######################
SUBROUTINE READ_GRIB_WG_ECMWF(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
! #######################
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
USE MODD_GRID_GRIB, ONLY : NNI
USE MODD_SURF_PAR, ONLY : XUNDEF
@@ -1257,7 +1345,7 @@ IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_WG_ECMWF',0,ZHOOK_HANDLE)
!
! 1. Get soil type to compute SWI
! ----------------------------
- CALL READ_GRIB(HGRIB,KLUOUT,43,IRET,ZSLT)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,43,IRET,ZSLT)
!--------------------------------------------------------------------------------
ALLOCATE (ZWFC(SIZE(PFIELD,1)))
ALLOCATE (ZWWILT(SIZE(PFIELD,1)))
@@ -1321,6 +1409,8 @@ END SUBROUTINE READ_GRIB_WG_ECMWF
! #######################
SUBROUTINE READ_GRIB_WGI_ECMWF(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
! #######################
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
USE MODD_GRID_GRIB, ONLY : NNI
USE MODD_SURF_PAR, ONLY : XUNDEF
@@ -1349,7 +1439,7 @@ IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_WGI_ECMWF',0,ZHOOK_HANDLE)
!
! 1. Get soil type to compute WSAT
! ----------------------------
- CALL READ_GRIB(HGRIB,KLUOUT,43,IRET,ZSLT)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,43,IRET,ZSLT)
!--------------------------------------------------------------------------------
ALLOCATE (ZWSAT(SIZE(PFIELD,1)))
ZWSAT(:)=0.
@@ -1458,6 +1548,8 @@ END SUBROUTINE READ_GRIB_TG_METEO_FRANCE
! #######################
SUBROUTINE READ_GRIB_SAND_CLAY_METEO_FRANCE(HGRIB,KLUOUT,HINMODEL,PSAND,PCLAY,GISBA)
! ######################
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
IMPLICIT NONE
!
@@ -1481,7 +1573,7 @@ IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_SAND_CLAY_METEO_FRANCE',0,ZHOO
!
IF (HINMODEL == 'ARPEGE' .OR. HINMODEL == 'MOCAGE') IPAR=171
IF (HINMODEL == 'ALADIN') IPAR=128
- CALL READ_GRIB(HGRIB,KLUOUT,IPAR,IRET,PCLAY)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPAR,IRET,PCLAY)
!
! if not available, the model is not ISBA (IWMODE=1)
IF (IRET /= 0) THEN
@@ -1496,7 +1588,7 @@ END IF
! ------------------------------------------
IF (HINMODEL == 'ARPEGE' .OR. HINMODEL == 'MOCAGE') IPAR=172
IF (HINMODEL == 'ALADIN') IPAR=129
- CALL READ_GRIB(HGRIB,KLUOUT,IPAR,IRET,PSAND)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPAR,IRET,PSAND)
!
! if not available, the model is not ISBA (IWMODE=1)
IF (GISBA) THEN
@@ -1513,8 +1605,10 @@ END SUBROUTINE READ_GRIB_SAND_CLAY_METEO_FRANCE
! #######################
SUBROUTINE READ_GRIB_WG_METEO_FRANCE(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
! #######################
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
-USE MODD_GRID_GRIB, ONLY : NNI
+USE MODD_GRID_GRIB, ONLY : NNI,NGRIB_VERSION
USE MODD_SURF_PAR, ONLY : XUNDEF
!
IMPLICIT NONE
@@ -1535,6 +1629,7 @@ INTEGER(KIND=kindOfInt) :: IRET ! return code
INTEGER :: ILTYPE ! type of level (Grib code table 3)
INTEGER :: ILEV1 ! level definition
INTEGER :: ILEV2 ! level definition
+INTEGER :: IPARAM ! number of ParamId or IndicatorOfParameter or ParameterNumber
REAL, DIMENSION(:), POINTER :: ZCLAY => NULL() ! clay fraction
REAL, DIMENSION(:), POINTER :: ZSAND => NULL() ! sand fraction
REAL, DIMENSION(:), POINTER :: ZFIELD => NULL()
@@ -1565,13 +1660,19 @@ PD(:,2) = 0.20
! ---------------------
ILEV1 = 0
IF (HINMODEL == 'ARPEGE' .OR. HINMODEL=='MOCAGE') THEN
- ILTYPE = 112
- ILEV2 = 1
- CALL READ_GRIB(HGRIB,KLUOUT,153,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
+ IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
+ IPARAM=20
+ ILTYPE=1
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,ZFIELD,KLTYPE=ILTYPE)
+ ELSE
+ IPARAM=153
+ ILEV2 = 1
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
+ ENDIF
ELSE
ILTYPE = 105
ILEV2 = 0
- CALL READ_GRIB(HGRIB,KLUOUT,86,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,86,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
ENDIF
IF (IRET /= 0) THEN
CALL ABOR1_SFX('MODE_READ_GRIB: SOIL MOISTURE LEVEL 1 MISSING (READ_GRIB_WG_METEO_FRANCE)')
@@ -1582,15 +1683,21 @@ PFIELD(:,1) = ZFIELD(:)
! 3. Read layer 2 moisture
! ---------------------
IF (HINMODEL == 'ARPEGE' .OR. HINMODEL=='MOCAGE') THEN
- ILTYPE = 112
- ILEV1 = 0
- ILEV2 = 250
- CALL READ_GRIB(HGRIB,KLUOUT,153,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
+ IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
+ IPARAM=20
+ ILTYPE=106
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,ZFIELD,KLTYPE=ILTYPE)
+ ELSE
+ IPARAM=153
+ ILEV1 = 0
+ ILEV2 = 250
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
+ ENDIF
ELSE
ILTYPE = 111
ILEV1 = -1
ILEV2 = -1
- CALL READ_GRIB(HGRIB,KLUOUT,86,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,86,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
ENDIF
IF (IRET /= 0) THEN
CALL ABOR1_SFX('MODE_READ_GRIB: SOIL MOISTURE LEVEL 2 MISSING (READ_GRIB_WG_METEO_FRANCE)')
@@ -1603,9 +1710,9 @@ PFIELD(:,2) = ZFIELD(:)
!* note that soil water reservoir is considered uniform between 0.2m and GRIB soil depth
IF (GISBA) THEN
IF (HINMODEL == 'ARPEGE' .OR. HINMODEL == 'MOCAGE') THEN
- CALL READ_GRIB(HGRIB,KLUOUT,173,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,173,IRET,ZFIELD)
ELSE
- CALL READ_GRIB(HGRIB,KLUOUT,130,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,130,IRET,ZFIELD)
ENDIF
IF (IRET /= 0) THEN
CALL ABOR1_SFX('MODE_READ_GRIB: SOIL MOISTURE LEVEL 2 DEPTH MISSING (READ_GRIB_WG_METEO_FRANCE)')
@@ -1665,8 +1772,10 @@ END SUBROUTINE READ_GRIB_WG_METEO_FRANCE
! #######################
SUBROUTINE READ_GRIB_WGI_METEO_FRANCE(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
! #######################
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
-USE MODD_GRID_GRIB, ONLY : NNI
+USE MODD_GRID_GRIB, ONLY : NNI,NGRIB_VERSION
USE MODD_SURF_PAR, ONLY : XUNDEF
!
IMPLICIT NONE
@@ -1685,6 +1794,7 @@ REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
LOGICAL :: GISBA ! T: surface scheme in file is ISBA
INTEGER(KIND=kindOfInt) :: IRET ! return code
INTEGER :: ILTYPE ! type of level (Grib code table 3)
+INTEGER :: IPARAM ! number of ParamId or IndicatorOfParameter or ParameterNumber
INTEGER :: ILEV1 ! level definition
INTEGER :: ILEV2 ! level definition
REAL, DIMENSION(:), POINTER :: ZCLAY => NULL() ! clay fraction
@@ -1714,13 +1824,20 @@ PD(:,1) = 0.01
! ---------------------
ILEV1 = 0
IF (HINMODEL == 'ARPEGE' .OR. HINMODEL=='MOCAGE') THEN
- ILTYPE = 112
- ILEV2 = 1
- CALL READ_GRIB(HGRIB,KLUOUT,152,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
+ IPARAM = 22
+ ILTYPE = 1
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,ZFIELD,KLTYPE=ILTYPE)
+ ELSE
+ IPARAM = 152
+ ILTYPE = 112
+ ILEV2 = 1
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ ENDIF
ELSE
ILTYPE = 105
ILEV2 = 0
- CALL READ_GRIB(HGRIB,KLUOUT,139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
END IF
!
IF (IRET == 0) THEN
@@ -1733,15 +1850,22 @@ END IF
! 3. Read layer 2 soil ice
! ---------------------
IF (HINMODEL == 'ARPEGE' .OR. HINMODEL=='MOCAGE') THEN
- ILTYPE = 112
- ILEV1 = 0
- ILEV2 = 250
- CALL READ_GRIB(HGRIB,KLUOUT,152,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
+ IPARAM = 22
+ ILTYPE = 106
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,ZFIELD,KLTYPE=ILTYPE)
+ ELSE
+ IPARAM = 152
+ ILTYPE = 112
+ ILEV1 = 0
+ ILEV2 = 250
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,IPARAM,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ ENDIF
ELSE
ILTYPE = 111
ILEV1 = -1
ILEV2 = -1
- CALL READ_GRIB(HGRIB,KLUOUT,139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
END IF
!
IF (IRET == 0) THEN
@@ -1755,9 +1879,9 @@ END IF
! ------------------------------
IF (GISBA) THEN
IF (HINMODEL == 'ARPEGE' .OR. HINMODEL=='MOCAGE') THEN
- CALL READ_GRIB(HGRIB,KLUOUT,173,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,173,IRET,ZFIELD)
ELSE
- CALL READ_GRIB(HGRIB,KLUOUT,130,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,130,IRET,ZFIELD)
ENDIF
IF (IRET /= 0) THEN
CALL ABOR1_SFX('MODE_READ_GRIB: SOIL ICE LEVEL 2 MISSING (READ_GRIB_WGI_METEO_FRANCE)')
@@ -1806,6 +1930,8 @@ END SUBROUTINE READ_GRIB_WGI_METEO_FRANCE
! #######################
SUBROUTINE READ_GRIB_TG_HIRLAM(HGRIB,KLUOUT,HINMODEL,PMASK,PTG,PDT)
! #######################
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
USE MODD_GRID_GRIB, ONLY : NNI
USE MODD_SURF_PAR, ONLY : XUNDEF
@@ -1837,7 +1963,7 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_TG_HIRLAM: | Reading soil temper
! -----------------------
ILEV1 = 904
ILEV2 = -1
- CALL READ_GRIB(HGRIB,KLUOUT,11,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,11,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
IF (IRET /= 0 ) THEN
CALL ABOR1_SFX('MODE_READ_GRIB: SOIL TEMPERATURE LEVEL 1 MISSING (READ_GRIB_TG_HIRLAM)')
END IF
@@ -1853,7 +1979,7 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_TG_HIRLAM: | Reading deep soil te
!
ILEV1 = 954
ILEV2 = -1
- CALL READ_GRIB(HGRIB,KLUOUT,11,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,11,IRET,ZFIELD,KLEV1=ILEV1,KLEV2=ILEV2)
IF (IRET /= 0) THEN
CALL ABOR1_SFX('MODE_READ_GRIB: DEEP SOIL TEMPERATURE MISSING (READ_GRIB_TG_HIRLAM)')
END IF
@@ -1881,6 +2007,8 @@ END SUBROUTINE READ_GRIB_TG_HIRLAM
! #######################
SUBROUTINE READ_GRIB_WG_HIRLAM(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
! #######################
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
USE MODD_GRID_GRIB, ONLY : NNI
USE MODD_SURF_PAR, ONLY : XUNDEF
@@ -1923,7 +2051,7 @@ ALLOCATE(ZD(2))
ILTYPE=105
ILEV1=904
ILEV2=-1
- CALL READ_GRIB(HGRIB,KLUOUT,86,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,86,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
IF (IRET /= 0 ) THEN
CALL ABOR1_SFX('MODE_READ_GRIB: SOIL MOISTURE LEVEL 1 MISSING (READ_GRIB_WG_HIRLAM)')
END IF
@@ -1941,7 +2069,7 @@ ZWG(:,1) = ZWG(:,1) / ZD(1) ! convert units to m3/m3
ILTYPE=105
ILEV1=954
ILEV2=-1
- CALL READ_GRIB(HGRIB,KLUOUT,86,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,86,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2)
IF (IRET /= 0 ) THEN
CALL ABOR1_SFX('MODE_READ_GRIB: SOIL MOISTURE LEVEL 2 MISSING (READ_GRIB_WG_HIRLAM)')
END IF
@@ -1978,6 +2106,8 @@ END SUBROUTINE READ_GRIB_WG_HIRLAM
! #######################
SUBROUTINE READ_GRIB_WGI_HIRLAM(HGRIB,KLUOUT,PFIELD,PD)
! #######################
+! MODIFICATIONS
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
USE MODD_GRID_GRIB, ONLY : NNI
!
@@ -2014,10 +2144,11 @@ END SUBROUTINE READ_GRIB_WGI_HIRLAM
!! MODIFICATIONS
!! -------------
!! C Ardilouze 07/2013 : possibility to read snow density (ERAI-land)
+!! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!!
!-------------------------------------------------------------------------------
!
-USE MODD_GRID_GRIB, ONLY : NNI
+USE MODD_GRID_GRIB, ONLY : NNI,NGRIB_VERSION
USE MODD_SURF_PAR, ONLY : XUNDEF
USE MODD_SNOW_PAR, ONLY : XRHOSMAX
!
@@ -2045,11 +2176,15 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_SNOW_VEG_AND_DEPTH: | Reading sn
!
SELECT CASE(HINMODEL)
CASE('ECMWF ')
- CALL READ_GRIB(HGRIB,KLUOUT,141,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,141,IRET,ZFIELD)
CASE('ARPEGE','ALADIN','MOCAGE','HIRLAM')
- CALL READ_GRIB(HGRIB,KLUOUT,66,IRET,ZFIELD)
+ IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,11,IRET,ZFIELD)
+ ELSE
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,66,IRET,ZFIELD)
+ ENDIF
CASE('NCEP ')
- CALL READ_GRIB(HGRIB,KLUOUT,3066,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,3066,IRET,ZFIELD)
CASE DEFAULT
CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_SNOW_VEG_AND_DEPTH: OPTION NOT SUPPORTED '//HINMODEL)
END SELECT
@@ -2091,6 +2226,7 @@ END SUBROUTINE READ_GRIB_SNOW_VEG_AND_DEPTH
!! AUTHOR
!! -------------
!! C Ardilouze 07/2013 : possibility to read snow albedo (ERAI-land)
+!! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!!
!-------------------------------------------------------------------------------
!
@@ -2120,7 +2256,7 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_SNOW_ALB: | Reading snow albedo'
ALLOCATE(PSNVA(NNI))
PSNVA(:) = 0.5 * ( XANSMIN + XANSMAX )
IF (HINMODEL == 'ECMWF') THEN
- CALL READ_GRIB(HGRIB,KLUOUT,32,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,32,IRET,ZFIELD)
IF (IRET == 0 ) THEN
DEALLOCATE(PSNVA)
ALLOCATE(PSNVA(SIZE(ZFIELD)))
@@ -2139,6 +2275,7 @@ END SUBROUTINE READ_GRIB_SNOW_ALB
!! AUTHOR
!! -------------
!! C Ardilouze 08/2013 : possibility to read snow density (ERAI-land)
+!! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!!
!-------------------------------------------------------------------------------
!
@@ -2168,7 +2305,7 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_SNOW_DEN: | Reading snow density
ALLOCATE(PSNV(NNI))
PSNV(:) = XRHOSMAX
IF (HINMODEL == 'ECMWF') THEN
- CALL READ_GRIB(HGRIB,KLUOUT,33,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,33,IRET,ZFIELD)
IF (IRET == 0 ) THEN
DEALLOCATE(PSNV)
ALLOCATE(PSNV(SIZE(ZFIELD)))
@@ -2298,6 +2435,7 @@ END SUBROUTINE READ_GRIB_TF_TEB
! #######################
SUBROUTINE READ_GRIB_TG_NCEP(HGRIB,KLUOUT,HINMODEL,PMASK,PTG,PD)
! #######################
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
IMPLICIT NONE
!
@@ -2313,7 +2451,7 @@ REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
!* local variables
! ---------------
INTEGER(KIND=kindOfInt) :: IRET ! return code
-INTEGER :: ILTYPE,ITER ! type of level (Grib code table 3)
+INTEGER :: ILTYPE ! type of level (Grib code table 3)
INTEGER :: ILEV1 ! level definition
INTEGER :: ILEV2 ! level definition
INTEGER :: JL ! layer loop counter
@@ -2341,10 +2479,9 @@ ILEV1 = -2
ILEV2 = -2
ZLEV1= 0.0
ZLEV2= 0.1
-CALL READ_GRIB(HGRIB,KLUOUT,228139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
+CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,&
+ HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
!
-
-
IF (IRET== 0) THEN
ZD(1)=(ZLEV2 - ZLEV1)
@@ -2364,8 +2501,8 @@ ZLEV1= 1.0
ZLEV2= 2.0
!
-CALL READ_GRIB(HGRIB,KLUOUT,228139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
-
+CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,&
+ HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
IF (IRET == 0) THEN
INLAYERDEEP = 4
@@ -2386,9 +2523,9 @@ ILEV2 = -2
ZLEV1= 0.4
ZLEV2= 1.0
- CALL READ_GRIB(HGRIB,KLUOUT,228139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,&
+ HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
!
-
IF (IRET == 0) THEN
ZD(3)=(ZLEV2 - ZLEV1)
ZTG(:,3)=ZFIELD
@@ -2404,10 +2541,9 @@ ILEV1 = -2
ILEV2 = -2
ZLEV1= 0.1
ZLEV2= 0.4
- CALL READ_GRIB(HGRIB,KLUOUT,228139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228139,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,&
+ HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
!
-
-
IF (IRET== 0) THEN
ZD(2)=(ZLEV2 - ZLEV1)
ZTG(:,2)=ZFIELD
@@ -2431,21 +2567,17 @@ IF(SUM(ZD(1:INLAYERDEEP)) < 3.) THEN
ZTG(:,INLAYERDEEP)=ZTG(:,INLAYERDEEP-1)
ENDIF
-
-
!
!--------------------------------------------------------------------------------
! 6. Set temperature profile and layer thicknesses
! ----------------------------------------------
-
+!
WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB: FILL PROFILE'
-
+!
CALL FILL_PFIELD(KLUOUT,'READ_GRIB_TG_NCEP',INLAYERDEEP,ZD,ZTG,PMASK,PTG,PD)
DEALLOCATE(ZD)
DEALLOCATE(ZTG)
!
-
-
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TG_NCEP',1,ZHOOK_HANDLE)
WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:END READ GRIB'
@@ -2455,6 +2587,7 @@ END SUBROUTINE READ_GRIB_TG_NCEP
! #######################
SUBROUTINE READ_GRIB_WG_NCEP(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
! #######################
+! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2
!
USE MODD_GRID_GRIB, ONLY : NNI
USE MODD_SURF_PAR, ONLY : XUNDEF
@@ -2504,7 +2637,8 @@ ILEV1 = -2
ILEV2 = -2
ZLEV1= 0.0
ZLEV2= 0.1
-CALL READ_GRIB(HGRIB,KLUOUT,260185,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
+CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,260185,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,&
+ HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
!
IF (IRET== 0) THEN
ZD(1)=(ZLEV2 - ZLEV1)
@@ -2524,7 +2658,8 @@ ZLEV1= 1.0
ZLEV2= 2.0
!
-CALL READ_GRIB(HGRIB,KLUOUT,260185,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
+CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,260185,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,&
+ HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
IF (IRET == 0) THEN
INLAYERDEEP = 4
@@ -2545,7 +2680,8 @@ ILEV2 = -2
ZLEV1= 0.4
ZLEV2= 1.0
- CALL READ_GRIB(HGRIB,KLUOUT,260185,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,260185,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,&
+ HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
!
IF (IRET == 0) THEN
ZD(3)=(ZLEV2 - ZLEV1)
@@ -2562,7 +2698,8 @@ ILEV1 = -2
ILEV2 = -2
ZLEV1= 0.1
ZLEV2= 0.4
- CALL READ_GRIB(HGRIB,KLUOUT,260185,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
+ CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,260185,IRET,ZFIELD,KLTYPE=ILTYPE,KLEV1=ILEV1,KLEV2=ILEV2,&
+ HTYPELEVEL=YTYPELEVEL,PLEV1=ZLEV1,PLEV2=ZLEV2)
!
IF (IRET== 0) THEN
ZD(2)=(ZLEV2 - ZLEV1)
diff --git a/src/SURFEX/prep_grib_grid.F90 b/src/SURFEX/prep_grib_grid.F90
index 7d148175e1c36065beab673558e926cf2e459e53..9461f2e4015f0d8afca12bdcadf43e548b8f9ab1 100644
--- a/src/SURFEX/prep_grib_grid.F90
+++ b/src/SURFEX/prep_grib_grid.F90
@@ -36,6 +36,8 @@
!! S. Faroux 01/2011 : to use library GRIB_API instead of GRIBEX (from
!! read_all_data_grib_case)
!! J. Escobar 1/06/2018 : remove IALLGRIB loop , useless & very memory consuming
+!! Gaelle Delautier (via Q.Rodier) 01/2019 : add GRIB 2 and adaptation to EPYGRAM-edited GRIB
+!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -56,7 +58,7 @@ USE MODD_GRID_ROTLATLON
USE MODD_GRID_GAUSS, ONLY : XILA1, XILO1, XILA2, XILO2, NINLA, NINLO, NILEN, LROTPOLE, XCOEF, XLAP, XLOP, &
XLAT, XLON
USE MODD_GRID_AROME, ONLY : XX, XY, NX, NY, XLAT0, XLON0, XLATOR, XLONOR, XRPK, XBETA, XZX, XZY, NIX
-USE MODD_GRID_GRIB, ONLY : NNI, CGRIB_FILE
+USE MODD_GRID_GRIB, ONLY : NNI, CGRIB_FILE,NGRIB_VERSION
USE MODD_SURF_PAR, ONLY : XUNDEF, NUNDEF
USE MODD_CSTS, ONLY : XPI
!
@@ -113,6 +115,9 @@ INTEGER :: JLOOP1 ! Dummy counter
!JUAN
!JUAN
INTEGER :: INFOMPI, J
+INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE :: INLO_GRIB ! Number of points along a parallel
+
+!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!---------------------------------------------------------------------------------------
@@ -148,6 +153,11 @@ END IF
IF (IRET /= 0) THEN
CALL ABOR1_SFX('PREP_GRIB_GRID: Error in reading center')
END IF
+!
+ CALL GRIB_GET(IGRIB,'editionNumber',NGRIB_VERSION,IRET)
+IF (IRET /= 0) THEN
+ CALL ABOR1_SFX('PREP_GRIB_GRID: Error in reading editionNumber')
+END IF
!
CALL GRIB_GET(IGRIB,'typeOfGrid',HGRID,IRET)
IF (IRET /= 0) THEN
@@ -464,14 +474,32 @@ SELECT CASE (HGRIDTYPE)
!
HGRIDTYPE = 'GAUSS '
IF (NRANK==NPIO) THEN
+! PROBLEME AVEC LES GRIB d'EPYGRAM
+! dans longitudeOfLastGridPointInDegrees la longitude du dernier point du
+! tableau (donc au pole sud)
+! dans les GRIB1 on a la valeur max du tableau des longitudes (donc Ã
+! l'equateur) REMARQUE : c'est ce qu'on a pour le CEP
+! pour éviter de changer le GRIB on calcule différemment XILO2 en le recalculant
+! le max de la longitude au niveau de l'équateur garce au champ 'pl'
+!
CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',XILA1)
CALL GRIB_GET(IGRIB,'longitudeOfFirstGridPointInDegrees',XILO1)
CALL GRIB_GET(IGRIB,'latitudeOfLastGridPointInDegrees',XILA2)
CALL GRIB_GET(IGRIB,'longitudeOfLastGridPointInDegrees',XILO2)
-
CALL GRIB_GET(IGRIB,'stretchingFactor',XCOEF)
CALL GRIB_GET(IGRIB,'latitudeOfStretchingPoleInDegrees',XLAP)
CALL GRIB_GET(IGRIB,'longitudeOfStretchingPoleInDegrees',XLOP)
+
+ IF (NRANK==NPIO) THEN
+ ALLOCATE (ININLO_GRIB(NINLA))
+ CALL GRIB_IS_MISSING(IGRIB,'pl',IMISSING,IRET)
+ IF (IRET == 0 .OR. IMISSING/=1) THEN ! quasi-regular
+ CALL GRIB_GET(IGRIB,'pl',ININLO_GRIB)
+ XILO2=360.-360./(MAXVAL(INLO_GRIB))
+ print*,"XILO2=",XILO2
+ ENDIF
+ DEALLOCATE(ININLO_GRIB)
+ ENDIF
ENDIF
LROTPOLE = .TRUE.
IF (NPROC>1) THEN