diff --git a/src/MNH/modd_prep_real.f90 b/src/MNH/modd_prep_real.f90
index d2d06321744568f7b2f1caf138f96388ff092c75..6933ae26a8f64486b93121fa44923c6a8a67774a 100644
--- a/src/MNH/modd_prep_real.f90
+++ b/src/MNH/modd_prep_real.f90
@@ -28,6 +28,7 @@
!! Original 05/05
!! 05/06 (I.Mallet) add *_SV_* variables to allow chemical
!! initialization from HCHEMFILE
+!! 09/20 (Q. Rodier) add geopotential height for GFS GRIB read
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -70,6 +71,7 @@ REAL :: XLEN2_LS ! Decay scale for small-scale
REAL,DIMENSION(:,:), ALLOCATABLE :: XPS_LS ! surface pressure
REAL,DIMENSION(:,:), ALLOCATABLE :: XZS_LS ! orography
REAL,DIMENSION(:,:), ALLOCATABLE :: XZSMT_LS ! smooth orography
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: XGH_LS ! geopotential height
REAL,DIMENSION(:,:,:), ALLOCATABLE :: XZFLUX_LS! altitude of pressure points
REAL,DIMENSION(:,:,:), ALLOCATABLE :: XZMASS_LS! altitude of mass points
REAL,DIMENSION(:,:,:), ALLOCATABLE :: XPMHP_LS ! pressure minus hyd. pressure
diff --git a/src/MNH/read_all_data_grib_case.f90 b/src/MNH/read_all_data_grib_case.f90
index 882dc12bb5ac3f637c8a6d139022b98d860bf931..bc171156015b1a1ad9cad8f8a3d48f5487afa48f 100644
--- a/src/MNH/read_all_data_grib_case.f90
+++ b/src/MNH/read_all_data_grib_case.f90
@@ -132,6 +132,7 @@ END MODULE MODI_READ_ALL_DATA_GRIB_CASE
! P. Wautelet 14/03/2019: correct ZWS when variable not present in file
! Q. Rodier 27/01/2020: switch of GRIB number ID for Orograpgy and hydrometeors in ARPEGE/AROME in EPyGrAM v1.3.7
! Q. Rodier 21/04/2020: correction GFS u and v wind component written in the right vertical order
+! Q. Rodier 02/09/2020 : Read and interpol geopotential height for interpolation on isobaric surface Grid of NCEP
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -203,7 +204,7 @@ INTEGER :: ILUOUT0 ! Unit used for output msg.
INTEGER :: IRESP ! Return code of FM-routines
INTEGER :: IRET ! Return code from subroutines
INTEGER(KIND=kindOfInt) :: IRET_GRIB ! Return code from subroutines
-INTEGER, PARAMETER :: JP_GFS=26 ! number of pressure levels for GFS model
+INTEGER, PARAMETER :: JP_GFS=31 ! number of pressure levels for GFS model
REAL :: ZA,ZB,ZC ! Dummy variables
REAL :: ZD,ZE,ZF ! |
REAL :: ZTEMP ! |
@@ -311,6 +312,7 @@ REAL, DIMENSION(:,:), ALLOCATABLE :: ZPF_G ! Pressure (flux point)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZPM_G ! Pressure (mass point)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZEXNF_G ! Exner fct. (flux point)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZEXNM_G ! Exner fct. (mass point)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZGH_G ! Geopotential Height
REAL, DIMENSION(:,:), ALLOCATABLE :: ZT_G ! Temperature
REAL, DIMENSION(:,:), ALLOCATABLE :: ZQ_G ! Specific humidity
REAL, DIMENSION(:), ALLOCATABLE :: ZH_G ! Relative humidity
@@ -334,7 +336,7 @@ INTEGER :: IVERSION,ILEVTYPE
LOGICAL :: GFIND ! to test if sea wave height is found
!---------------------------------------------------------------------------------------
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/)!
+ 250,200,150,100,70,50,30,20,10,7,5,3,2,1/)!
!
TZFILE => NULL()
!
@@ -531,15 +533,10 @@ SELECT CASE (IMODEL)
END IF
ENDIF
CASE(10) ! NCEP
- DO IVAR=0,222
- CALL SEARCH_FIELD(IGRIB,INUM_ZS,KPARAM=IVAR)
+ CALL SEARCH_FIELD(IGRIB,INUM_ZS,KDIS=0,KCAT=3,KNUMBER=5,KTFFS=1)
IF(INUM_ZS < 0) THEN
- WRITE (ILUOUT0,'(A)')'Orography is missing'
- ENDIF
- END DO
- INUM_ZS=218
- WRITE (ILUOUT0,*) 'lsm ',IGRIB(350)
- WRITE (ILUOUT0,*) 'orog ',IGRIB(INUM_ZS)
+ WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
+ ENDIF
END SELECT
ZPARAM(:)=-999.
CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ,IRET_GRIB)
@@ -755,7 +752,7 @@ IF (IMODEL/=10) THEN
STOP
ENDIF
ELSE ! NCEP
- ISTARTLEVEL=10
+ ISTARTLEVEL=1000
IT=130
IQ=157
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ISTARTLEVEL)
@@ -813,7 +810,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(IGRIB,INUM,KPARAM=IT,KLEV1=ILEV1)
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=0,KNUMBER=0,KLEV1=ILEV1,KTFFS=100)
IF (INUM< 0) THEN
!callabortstop
WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
@@ -1055,6 +1052,43 @@ ELSE !NCEP
END IF
DEALLOCATE (ZOUT)
+
+
+!---------------------------------------------------------------------------------------
+!* 2.5.4.2 Read and interpol geopotential height for interpolation on isobaric surface Grid of NCEP
+!---------------------------------------------------------------------------------------
+!
+ALLOCATE (ZGH_G(ISIZE,INLEVEL))
+!
+IF(IMODEL==10) THEN !NCEP with pressure grid only
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = IP_GFS(JLOOP1)
+ WRITE (ILUOUT0,'(A)') ' | Searching geopotential height'
+ CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=3,KNUMBER=5,KLEV1=ILEV1)
+ IF (INUM< 0) THEN
+ !callabortstop
+ WRITE(YMSG,*) 'Geopoential height level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ !
+ CALL GRIB_GET(IGRIB(INUM),'values',ZGH_G(:,JLOOP1),IRET_GRIB)
+ CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
+ !
+ END DO
+ !
+ CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM_ZS),IIU,IJU,ZLONOUT,ZLATOUT,&
+ ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
+ !
+ ALLOCATE(ZOUT(INO))
+ ALLOCATE(XGH_LS(IIU,IJU,INLEVEL))
+ !
+ DO JLOOP1=1, INLEVEL
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ ZGH_G(:,JLOOP1),INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+ CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XGH_LS(:,:,JLOOP1))
+ END DO
+ DEALLOCATE(ZOUT)
+END IF
!
!* 2.5.5 Compute atmospheric pressure on MESO-NH grid
!
@@ -1911,7 +1945,7 @@ END SUBROUTINE ARRAY_1D_TO_2D
!---------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------
!#################################################################################
-SUBROUTINE SEARCH_FIELD(KGRIB,KNUM,KPARAM,KDIS,KCAT,KNUMBER,KLEV1)
+SUBROUTINE SEARCH_FIELD(KGRIB,KNUM,KPARAM,KDIS,KCAT,KNUMBER,KLEV1,KTFFS)
!#################################################################################
! search the grib message corresponding to KPARAM,KLTYPE,KLEV1,KLEV2 in all
! the KGIRB messages
@@ -1931,13 +1965,14 @@ 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
+INTEGER,INTENT(IN),OPTIONAL :: KTFFS ! TypeOfFirstFixedSurface
!
! Declaration of local variables
!
INTEGER :: IFOUND ! Number of correct parameters
INTEGER :: ISEARCH ! Number of correct parameters to find
INTEGER :: IRET ! error code
-INTEGER :: IPARAM,IDIS,ICAT,INUMBER
+INTEGER :: IPARAM,IDIS,ICAT,INUMBER,ITFFS
INTEGER :: ILEV1 ! Level parameter 1
INTEGER :: JLOOP ! Dummy counter
INTEGER :: IVERSION
@@ -1955,6 +1990,7 @@ IF (PRESENT(KDIS)) ISEARCH=ISEARCH+1
IF (PRESENT(KCAT)) ISEARCH=ISEARCH+1
IF (PRESENT(KNUMBER)) ISEARCH=ISEARCH+1
IF (PRESENT(KLEV1)) ISEARCH=ISEARCH+1
+IF(PRESENT(KTFFS)) ISEARCH=ISEARCH+1
!
DO JLOOP=1,SIZE(KGRIB)
IFOUND = 0
@@ -1967,6 +2003,23 @@ DO JLOOP=1,SIZE(KGRIB)
WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
CYCLE
ENDIF
+ !
+ IF (PRESENT(KTFFS)) THEN
+ CALL GRIB_GET(KGRIB(JLOOP),'typeOfFirstFixedSurface',ITFFS,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 (ITFFS==KTFFS) THEN
+ IFOUND = IFOUND + 1
+ ELSE
+ CYCLE
+ ENDIF
+ ENDIF
+
IF (PRESENT(KPARAM)) THEN
IF (IVERSION == 2) THEN
CALL GRIB_GET(KGRIB(JLOOP),'paramId',IPARAM,IRET_GRIB)
diff --git a/src/MNH/ver_prep_gribex_case.f90 b/src/MNH/ver_prep_gribex_case.f90
index 4368fc4c10aabfd7cfdaa393f936683276450f38..593ba064da6e33f3c6f81e51ee300baf6f00eeae 100644
--- a/src/MNH/ver_prep_gribex_case.f90
+++ b/src/MNH/ver_prep_gribex_case.f90
@@ -90,6 +90,8 @@ END MODULE MODI_VER_PREP_GRIBEX_CASE
!! May 2006 Remove EPS
!! Apr, 09 2018 (J.-P. Chaboureau) add isobaric surface
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Sep, 02, 2020 (Q. Rodier) use of geopotential height instead of
+!! height above orography for isobaric surface interpolation
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -276,8 +278,9 @@ IF (HFILE(1:3)=='ATM') THEN
ZU_LS,ZV_LS, &
ZW_LS,'FLUX' )
ELSE ! isobaric surfaces (w at mass points)
+ !Warning, in that case (NCEP only for now) ZZFLUX_LS is not correct (but not used)
CALL VER_INTERP_TO_MIXED_GRID('ATM ',.TRUE.,XZS_LS,XZSMT_LS, &
- ZZMASS_LS,ZSV_LS, &
+ XGH_LS,ZSV_LS, &
ZZFLUX_LS,XPS_LS,ZPMHP_LS, &
ZTHV_LS,ZR_LS, &
ZHU_LS, &
diff --git a/src/SURFEX/mode_read_grib.F90 b/src/SURFEX/mode_read_grib.F90
index 344ea5967b0658f1a362b4eef77366bc302691fa..786f33c2b7da63b340276d881c88c5127e1f8d74 100644
--- a/src/SURFEX/mode_read_grib.F90
+++ b/src/SURFEX/mode_read_grib.F90
@@ -514,7 +514,7 @@ 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)
+ 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,4,IRET,PZS)
@@ -536,7 +536,9 @@ IF (IRET /= 0) THEN
END IF
!
! Datas given in archives are multiplied by the gravity acceleration
-PZS = PZS / XG
+IF(HINMODEL /= 'NCEP') THEN
+ PZS = PZS / XG
+END IF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_ZS',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_ZS