From 5cc0c0ad148784877adcdd919a7a42b097166dd1 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Wed, 16 Sep 2020 16:17:01 +0200
Subject: [PATCH] Quentin 16/09/2020: TSWATER bugfix: major correction on
 reading the land mask for interpolation of inland temperature TSWATER from
 GRIB file for AROME, ARPEGE, IFS and MOCAGE. The bug is present since MNH541
 July 2018 when GFS was introduced. For all models, the land mask is 1=land ;
 0=sea/inland water. Then, the surface temperature of the lakes is XUNDEF for
 PMASK/=0. With PMASK/=1 the surface temperature of the lakes was equal to
 land surface temperature (T2 in the grib)

---
 src/SURFEX/mode_read_grib.F90 | 25 ++++---------------------
 1 file changed, 4 insertions(+), 21 deletions(-)

diff --git a/src/SURFEX/mode_read_grib.F90 b/src/SURFEX/mode_read_grib.F90
index 786f33c2b..2751bff06 100644
--- a/src/SURFEX/mode_read_grib.F90
+++ b/src/SURFEX/mode_read_grib.F90
@@ -63,11 +63,9 @@ 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)        
+   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")
@@ -136,20 +134,16 @@ 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)
@@ -170,7 +164,6 @@ print*,IRET
 
     ENDIF
     !
-print*,KFOUND
     IF (IRET.EQ.0) THEN
       !
       KFOUND = KFOUND + 1
@@ -185,7 +178,6 @@ print*,KFOUND
       !
       IF (IRET.EQ.0) KFOUND = KFOUND + 1
       !
-print*,KFOUND
     ENDIF
     !
   ENDIF
@@ -348,11 +340,9 @@ 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)         
+   CALL GRIB_INDEX_SELECT(NIDX,'paramId',KPARAM,KRET)
  ENDIF
 END IF
  CALL GRIB_NEW_FROM_INDEX(NIDX,IGRIB,KRET)
@@ -363,7 +353,6 @@ 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
@@ -390,8 +379,6 @@ IF (PRESENT(KPARAM2)) THEN
   ENDIF
 ENDIF
 !
-print*,"IFOUND=",IFOUND
-
 IF (IFOUND==3) THEN
   !
   IF (PRESENT(KLTYPE)) KLTYPE = ILTYPE
@@ -443,8 +430,6 @@ 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
 !
@@ -454,14 +439,12 @@ SELECT CASE (HINMODEL)
   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)   
+      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    
@@ -758,7 +741,7 @@ SELECT CASE (HINMODEL)
      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 (SIZE(PMASK)==SIZE(PTS)) WHERE ((PMASK(:)/=0.) .OR. ((PMASK(:)==0.) .AND.(PTS(:)==9999.))) PTS = XUNDEF
 
 !
 IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TSWATER',1,ZHOOK_HANDLE)
-- 
GitLab