From 231fe82f3ee58a82efd8818dd6c95b47c2a76325 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Thu, 3 Sep 2020 17:29:55 +0200
Subject: [PATCH] Quentin 03/09/2020: GFS bugfix: vertical interpolation
(ver_inter_to_mixed_grid.f90) from the GFS isobaric pressure grid was wrong.
Previously, the computation of the altitude of mass points XZMASS_LS
(ver_prep_gribex.f90) was done from XZFLUX_LS computed from the orography
PZS_LS (in change_gribex_var.f90) as for hybrid vertical coordinate. However,
the vertical grid of GFS is isobaric pressure (no orography). Thus,
previously the first level of GFS data (at 1000hPa) was interpolated at the
first grid level of MNH (orography + ZHAT) leading to strong biais of
temperature, wind and so on near high topography. The corrections include : -
reading of all levels (31), - add tyeOfFirstFixedSurface argument for
SEARCH_FIELD routine to read correctly the orography and geopotential height
from the GRIB, - reading the geopotential height in
read_all_data_grib_case.f90 for GFS and inject it in the case of vertical
interpolation of an isobaric pressure grid in ver_prep_gribex_case.
---
src/MNH/modd_prep_real.f90 | 2 +
src/MNH/read_all_data_grib_case.f90 | 81 ++++++++++++++++++++++++-----
src/MNH/ver_prep_gribex_case.f90 | 5 +-
src/SURFEX/mode_read_grib.F90 | 6 ++-
4 files changed, 77 insertions(+), 17 deletions(-)
diff --git a/src/MNH/modd_prep_real.f90 b/src/MNH/modd_prep_real.f90
index d2d063217..6933ae26a 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 882dc12bb..bc1711560 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 4368fc4c1..593ba064d 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 344ea5967..786f33c2b 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
--
GitLab