diff --git a/src/MNH/change_gribex_var.f90 b/src/MNH/change_gribex_var.f90
index af6f23f3d298a5608e474b5bc1aef457e91d97cd..48a9741de41c7fdac18ca15f731ed4bc1698715b 100644
--- a/src/MNH/change_gribex_var.f90
+++ b/src/MNH/change_gribex_var.f90
@@ -157,6 +157,7 @@ END MODULE MODI_CHANGE_GRIBEX_VAR
!! Masson 12/06/97 add relative humidity
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Pergaud : 2018 add GFS
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -241,27 +242,39 @@ CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
!* 1. COMPUTATION OF PRESSURE AND EXNER FUNCTION
! ------------------------------------------
!
-ALLOCATE(ZPFLUX_LS(IIU,IJU,ILU),ZEXNFLUX_LS(IIU,IJU,ILU))
+IF (SIZE(PB_LS)/=0) THEN ! hybrid level
+
+ ALLOCATE(ZPFLUX_LS(IIU,IJU,ILU),ZEXNFLUX_LS(IIU,IJU,ILU))
!
-ZPFLUX_LS(:,:,:)=SPREAD(SPREAD(PA_LS,1,IIU),2,IJU)*PP00_LS &
+ ZPFLUX_LS(:,:,:)=SPREAD(SPREAD(PA_LS,1,IIU),2,IJU)*PP00_LS &
+SPREAD(SPREAD(PB_LS,1,IIU),2,IJU)*SPREAD(PPS_LS,3,ILU)
!
-ZEXNFLUX_LS(:,:,:)=(ZPFLUX_LS(:,:,:)/XP00)**(XRD/XCPD)
+ ZEXNFLUX_LS(:,:,:)=(ZPFLUX_LS(:,:,:)/XP00)**(XRD/XCPD)
!
!-------------------------------------------------------------------------------
!
!* 2. COMPUTATION OF EXNER FUNCTION AT MASS POINT
! -------------------------------------------
!
-ALLOCATE(ZEXNMASS_LS(IIU,IJU,ILU))
+ ALLOCATE(ZEXNMASS_LS(IIU,IJU,ILU))
!
-ZEXNMASS_LS(:,:,1:ILU-1)=(ZEXNFLUX_LS(:,:,1:ILU-1)-ZEXNFLUX_LS(:,:,2:ILU)) &
+ ZEXNMASS_LS(:,:,1:ILU-1)=(ZEXNFLUX_LS(:,:,1:ILU-1)-ZEXNFLUX_LS(:,:,2:ILU)) &
/(LOG(ZEXNFLUX_LS(:,:,1:ILU-1))-LOG(ZEXNFLUX_LS(:,:,2:ILU)))
!
-ZEXNMASS_LS(:,:,ILU) =(ZPFLUX_LS(:,:,ILU)/2./XP00)**(XRD/XCPD)
+ ZEXNMASS_LS(:,:,ILU) =(ZPFLUX_LS(:,:,ILU)/2./XP00)**(XRD/XCPD)
!
-PPMASS_LS(:,:,:)=XP00*(ZEXNMASS_LS(:,:,:))**(XCPD/XRD)
-!-------------------------------------------------------------------------------
+ PPMASS_LS(:,:,:)=XP00*(ZEXNMASS_LS(:,:,:))**(XCPD/XRD)
+
+ELSE
+ PPMASS_LS(:,:,:)=SPREAD(SPREAD(PA_LS,1,IIU),2,IJU)
+ ALLOCATE(ZEXNMASS_LS(IIU,IJU,ILU))
+ ZEXNMASS_LS(:,:,:)=(PPMASS_LS(:,:,:)/XP00)**(XRD/XCPD)
+
+ ALLOCATE(ZEXNFLUX_LS(IIU,IJU,ILU))
+ ZEXNFLUX_LS(:,:,1:ILU-1)=(ZEXNMASS_LS(:,:,1:ILU-1)-ZEXNMASS_LS (:,:,2:ILU)) &
+ /(LOG(ZEXNMASS_LS(:,:,1:ILU-1))-LOG(ZEXNMASS_LS (:,:,2:ILU)))
+ ZEXNFLUX_LS(:,:,ILU) =(PPMASS_LS(:,:,ILU)/2./XP00)**(XRD/XCPD)
+END IF!-------------------------------------------------------------------------------
!
!* 3. COMPUTATION OF RELATIVE HUMIDITY
! --------------------------------
@@ -326,7 +339,10 @@ PZMASS_LS(:,:,1) =PZFLUX_LS(:,:,2) - XCPD/XG &
!* 8. VERTICAL WIND
! -------------
!
-IF (PRESENT(PW_LS)) THEN
+IF (PRESENT(PW_LS) .AND. SIZE(PB_LS)==0) THEN
+ PW_LS(:,:,:)=0. ! NCEP case not treated
+
+ELSEIF (PRESENT(PW_LS)) THEN
!
!* 8.0 allocations
! -----------
diff --git a/src/MNH/read_all_data_grib_case.f90 b/src/MNH/read_all_data_grib_case.f90
index 67f66348eab5e1e9cc717150fcf74145c6e6f004..77235ddbfb7937a934d8013c186badecee6b40ab 100644
--- a/src/MNH/read_all_data_grib_case.f90
+++ b/src/MNH/read_all_data_grib_case.f90
@@ -126,6 +126,7 @@ END MODULE MODI_READ_ALL_DATA_GRIB_CASE
!! 05/12/2016 (G.Delautier) length of HGRID for grib_api > 1.14
!! 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
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -196,6 +197,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
REAL :: ZA,ZB,ZC ! Dummy variables
REAL :: ZD,ZE,ZF ! |
REAL :: ZTEMP ! |
@@ -233,6 +235,7 @@ INTEGER :: IMODEL ! Type of Grib file :
! 3 -> METEO FRANCE - ARPEGE
! 4 -> METEO FRANCE - ARPEGE
! 5 -> METEO FRANCE - MOCAGE
+ ! 10 -> NCEP - GFS
INTEGER :: ICENTER ! number of center
INTEGER :: ISIZE ! size of grib message
INTEGER(KIND=kindOfInt) :: ICOUNT ! number of messages in the file
@@ -319,8 +322,10 @@ INTEGER :: IPVPRESENT ,IPV
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZR_DUM
INTEGER :: IMI
TYPE(TFILEDATA),POINTER :: TZFILE
-!
+INTEGER, DIMENSION(JP_GFS) :: IP_GFS ! list of pressure levels for GFS model
!---------------------------------------------------------------------------------------
+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/)!
!
TZFILE => NULL()
!
@@ -465,7 +470,11 @@ SELECT CASE (ICENTER)
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Mocage model'
IMODEL = 5
ALLOCATE(ZPARAM(6))
- END SELECT
+ END SELECT
+ CASE (7)
+ WRITE (ILUOUT0,'(A)') ' | Grib file from National Center for Environmental Prediction'
+ IMODEL = 10
+ ALLOCATE(ZPARAM(6))
END SELECT
IF (IMODEL==-1) THEN
!callabortstop
@@ -496,6 +505,16 @@ SELECT CASE (IMODEL)
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
+ END DO
+ INUM_ZS=218
+ WRITE (ILUOUT0,*) 'lsm ',IGRIB(350)
+ WRITE (ILUOUT0,*) 'orog ',IGRIB(INUM_ZS)
END SELECT
ZPARAM(:)=-999.
CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ,IRET_GRIB)
@@ -515,8 +534,10 @@ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
DEALLOCATE(IINLO)
DEALLOCATE(ZVALUE)
!
-! Datas given in archives are multiplied by the gravity acceleration
-ZOUT = ZOUT / XG
+IF (IMODEL/=10) THEN ! others than NCEP
+ ! Data given in archives are multiplied by the gravity acceleration
+ ZOUT = ZOUT / XG
+END IF
!
! Stores the field in a 2 dimension array
IF (HFILE(1:3)=='ATM') THEN
@@ -543,6 +564,8 @@ SELECT CASE (IMODEL)
CALL SEARCH_FIELD(152,-1,-1,-1,IGRIB,INUM)
CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
CALL SEARCH_FIELD(1,-1,-1,-1,IGRIB,INUM)
+ CASE(10) ! NCEP
+ CALL SEARCH_FIELD(134,1,0,0,IGRIB,INUM)
END SELECT
IF(INUM < 0) THEN
WRITE (ILUOUT0,'(A)')'Surface pressure is missing - abort'
@@ -560,7 +583,7 @@ SELECT CASE (IMODEL)
ALLOCATE (ZLNPS_G(ISIZE))
ZLNPS_G(:) = ZVALUE(1:ISIZE)
ZPS_G (:) = EXP(ZVALUE(1:ISIZE))
- CASE(1,2,3,4,5) ! arpege mocage aladin aladin-reunion
+ CASE(1,2,3,4,5,10) ! arpege mocage aladin aladin-reunion NCEP
ALLOCATE (ZPS_G (ISIZE))
ALLOCATE (ZLNPS_G(ISIZE))
ZPS_G (:) = ZVALUE(1:ISIZE)
@@ -663,103 +686,158 @@ SELECT CASE (IMODEL)
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
!
-CALL GRIB_GET(IGRIB(INUM),'NV',INLEVEL)
-INLEVEL = NINT(INLEVEL / 2.) - 1
-CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
+IF (IMODEL/=10) THEN ! others than NCEP
+ CALL GRIB_GET(IGRIB(INUM),'NV',INLEVEL)
+ INLEVEL = NINT(INLEVEL / 2.) - 1
+ CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
+ELSE
+ INLEVEL=JP_GFS
+END IF
!
ALLOCATE (ZT_G(ISIZE,INLEVEL))
ALLOCATE (ZQ_G(ISIZE,INLEVEL))
!
-DO JLOOP1=1, INLEVEL
- ILEV1 = JLOOP1-1+ISTARTLEVEL
- CALL SEARCH_FIELD(IQ,109,ILEV1,-1,IGRIB,INUM)
- 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)
- IF (INUM< 0) THEN
+IF (IMODEL/=10) THEN ! others than NCEP
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = JLOOP1-1+ISTARTLEVEL
+ CALL SEARCH_FIELD(IQ,109,ILEV1,-1,IGRIB,INUM)
+ IF (INUM< 0) THEN
!callabortstop
- WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
- CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
- END IF
- CALL GRIB_GET(IGRIB(INUM),'values',ZT_G(:,INLEVEL-JLOOP1+1))
- CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
-END DO
+ 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)
+ IF (INUM< 0) THEN
+ !callabortstop
+ WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET(IGRIB(INUM),'values',ZT_G(:,INLEVEL-JLOOP1+1))
+ CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
+ END DO
+ELSE ! NCEP
+ DO JLOOP1=1, INLEVEL
+ ILEV1 = IP_GFS(JLOOP1)
+ CALL SEARCH_FIELD(IQ,100,ILEV1,-1,IGRIB,INUM)
+ 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(:,JLOOP1),IRET_GRIB)
+ WRITE (ILUOUT0,*) 'Q ',ILEV1,IRET_GRIB
+ CALL SEARCH_FIELD(IT,100,ILEV1,-1,IGRIB,INUM)
+ IF (INUM< 0) THEN
+ !callabortstop
+ WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+ END IF
+ CALL GRIB_GET(IGRIB(INUM),'values',ZT_G(:,JLOOP1),IRET_GRIB)
+ WRITE (ILUOUT0,*) 'T ',ILEV1,IRET_GRIB
+ CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
+ END DO
+END IF
+
ALLOCATE(IINLO(INJ))
CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
!
!* 2.5.2 Load level definition parameters A and B
!
-IF (HFILE(1:3)=='ATM') THEN
- XP00_LS = 101325.
-ELSE IF (HFILE=='CHEM') THEN
- XP00_SV_LS = 101325.
-END IF
-!
-IF (INLEVEL > 0) THEN
+IF (IMODEL/=10) THEN ! others than NCEP
+
IF (HFILE(1:3)=='ATM') THEN
- ALLOCATE (XA_LS(INLEVEL))
- ALLOCATE (XB_LS(INLEVEL))
+ XP00_LS = 101325.
ELSE IF (HFILE=='CHEM') THEN
- ALLOCATE (XA_SV_LS(INLEVEL))
- ALLOCATE (XB_SV_LS(INLEVEL))
+ XP00_SV_LS = 101325.
END IF
!
- CALL GRIB_GET(IGRIB(INUM),'PVPresent',IPVPRESENT)
- IF (IPVPRESENT==1) THEN
- CALL GRIB_GET_SIZE(IGRIB(INUM),'pv',IPV)
- ALLOCATE(ZPV(IPV))
- CALL GRIB_GET(IGRIB(INUM),'pv',ZPV)
- ELSE
- !callabortstop
- 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)
- DO JLOOP1 = 1, INLEVEL
- XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+JLOOP1) / XP00_LS
- XB_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
- END DO
- CASE (1,2)
- JLOOP2 = 2
- DO JLOOP1 = 1, INLEVEL
- JLOOP2 = JLOOP2 + 1
- XA_LS(1 + INLEVEL - JLOOP1) = ZPV(JLOOP2)
- JLOOP2 = JLOOP2 + 1
- XB_LS(1 + INLEVEL - JLOOP1) = ZPV(JLOOP2)
- END DO
- CASE (5)
- DO JLOOP1 = 1, INLEVEL
- IF (HFILE(1:3)=='ATM') THEN
- XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+ JLOOP1) / XP00_LS**2
+ IF (INLEVEL > 0) THEN
+ IF (HFILE(1:3)=='ATM') THEN
+ ALLOCATE (XA_LS(INLEVEL))
+ ALLOCATE (XB_LS(INLEVEL))
+ ELSE IF (HFILE=='CHEM') THEN
+ ALLOCATE (XA_SV_LS(INLEVEL))
+ ALLOCATE (XB_SV_LS(INLEVEL))
+ END IF
+!
+ CALL GRIB_GET(IGRIB(INUM),'PVPresent',IPVPRESENT)
+ IF (IPVPRESENT==1) THEN
+ CALL GRIB_GET_SIZE(IGRIB(INUM),'pv',IPV)
+ ALLOCATE(ZPV(IPV))
+ CALL GRIB_GET(IGRIB(INUM),'pv',ZPV)
+ ELSE
+ !callabortstop
+ 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)
+ DO JLOOP1 = 1, INLEVEL
+ XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+JLOOP1) / XP00_LS
XB_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
- ELSE IF (HFILE=='CHEM') THEN
- XA_SV_LS(1 + INLEVEL - JLOOP1) = ZPV(1+ JLOOP1) / XP00_LS**2
- XB_SV_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
- END IF
- END DO
- END SELECT
+ END DO
+ CASE (1,2)
+ JLOOP2 = 2
+ DO JLOOP1 = 1, INLEVEL
+ JLOOP2 = JLOOP2 + 1
+ XA_LS(1 + INLEVEL - JLOOP1) = ZPV(JLOOP2)
+ JLOOP2 = JLOOP2 + 1
+ XB_LS(1 + INLEVEL - JLOOP1) = ZPV(JLOOP2)
+ END DO
+ CASE (5)
+ DO JLOOP1 = 1, INLEVEL
+ IF (HFILE(1:3)=='ATM') THEN
+ XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+ JLOOP1) / XP00_LS**2
+ XB_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
+ ELSE IF (HFILE=='CHEM') THEN
+ XA_SV_LS(1 + INLEVEL - JLOOP1) = ZPV(1+ JLOOP1) / XP00_LS**2
+ XB_SV_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
+ END IF
+ END DO
+ END SELECT
+ ELSE
+ !callabortstop
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','level definition section is missing')
+ END IF
ELSE
- !callabortstop
- CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','level definition section is missing')
+ ALLOCATE (XA_LS(INLEVEL))
+ ALLOCATE (XB_LS(0))
+ XA_LS = 100.*IP_GFS
END IF
!
!* 2.5.3 Compute atmospheric pressure on grib grid
!
WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on Grib grid is being computed'
ALLOCATE (ZPF_G(INI,INLEVEL))
-IF (HFILE(1:3)=='ATM') THEN
+IF (IMODEL/=10) THEN ! others than NCEP
+ IF (HFILE(1:3)=='ATM') THEN
ZPF_G(:,:) = SPREAD(XA_LS,1,INI)*XP00_LS + &
SPREAD(XB_LS,1,INI)*SPREAD(ZPS_G(1:INI),2,INLEVEL)
-ELSE IF (HFILE=='CHEM') THEN
+ ELSE IF (HFILE=='CHEM') THEN
ZPF_G(:,:) = SPREAD(XA_SV_LS,1,INI)*XP00_SV_LS + &
SPREAD(XB_SV_LS,1,INI)*SPREAD(ZPS_G(1:INI),2,INLEVEL)
+ END IF
+ELSE
+ ZPF_G(:,:) = 100.*SPREAD(IP_GFS,1,INI)
END IF
DEALLOCATE (ZPS_G)
!
@@ -826,47 +904,78 @@ ALLOCATE (ZTHV_LS(IIU,IJU,INLEVEL))
ALLOCATE (ZTHV_G(INI))
ALLOCATE (ZRV_G(INI))
ALLOCATE (ZOUT(INO))
-DO JLOOP1=1, INLEVEL
- !
- ! Compute Theta V and relative humidity on grib grid
- !
- ! (1/rv) = (1/q) - 1
- ! Thetav = T . (P0/P)^(Rd/Cpd) . ( (1 + (Rv/Rd).rv) / (1 + rv) )
- ! Hu = P / ( ( (Rd/Rv) . ((1/rv) - 1) + 1) . Es(T) )
- !
- ZRV_G(:) = 1. / (1./MAX(ZQ_G(:,JLOOP1),1.E-12) - 1.)
- !
- ZTHV_G(:)=ZT_G(:,JLOOP1) * ((XP00/ZPM_G(:,JLOOP1))**(XRD/XCPD)) * &
+IF (IMODEL/=10) THEN ! others than NCEP
+ DO JLOOP1=1, INLEVEL
+ !
+ ! Compute Theta V and relative humidity on grib grid
+ !
+ ! (1/rv) = (1/q) - 1
+ ! Thetav = T . (P0/P)^(Rd/Cpd) . ( (1 + (Rv/Rd).rv) / (1 + rv) )
+ ! Hu = P / ( ( (Rd/Rv) . ((1/rv) - 1) + 1) . Es(T) )
+ !
+ ZRV_G(:) = 1. / (1./MAX(ZQ_G(:,JLOOP1),1.E-12) - 1.)
+ !
+ ZTHV_G(:)=ZT_G(:,JLOOP1) * ((XP00/ZPM_G(:,JLOOP1))**(XRD/XCPD)) * &
((1. + XRV*ZRV_G(:)/XRD) / (1. + ZRV_G(:)))
- !
- ZH_G(1:INI) = 100. * ZPM_G(:,JLOOP1) / ( (XRD/XRV)*(1./ZRV_G(:)+1.)*SM_FOES(ZT_G(:,JLOOP1)) )
- ZH_G(:) = MAX(MIN(ZH_G(:),100.),0.)
- !
- ! Interpolation : H
- CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ !
+ ZH_G(1:INI) = 100. * ZPM_G(:,JLOOP1) / ( (XRD/XRV)*(1./ZRV_G(:)+1.)*SM_FOES(ZT_G(:,JLOOP1)) )
+ ZH_G(:) = MAX(MIN(ZH_G(:),100.),0.)
+ !
+ ! Interpolation : H
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
ZH_G,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
- CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZH_LS(:,:,JLOOP1))
- ZH_LS(:,:,JLOOP1) = MAX(MIN(ZH_LS(:,:,JLOOP1),100.),0.)
- !
- ! interpolation : Theta V
- CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZH_LS(:,:,JLOOP1))
+ ZH_LS(:,:,JLOOP1) = MAX(MIN(ZH_LS(:,:,JLOOP1),100.),0.)
+ !
+ ! interpolation : Theta V
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
ZTHV_G,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
- CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZTHV_LS(:,:,JLOOP1))
- !
-END DO
+ CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZTHV_LS(:,:,JLOOP1))
+ !
+ END DO
+ELSE !NCEP
+ DO JLOOP1=1, INLEVEL
+ WRITE (ILUOUT0,*), 'JLOOP1=',JLOOP1,MINVAL(ZPM_G(:,JLOOP1)),MINVAL(ZT_G(:,JLOOP1)),MINVAL(ZQ_G(:,JLOOP1))
+ WRITE (ILUOUT0,*), ' ',MAXVAL(ZPM_G(:,JLOOP1)),MAXVAL(ZT_G(:,JLOOP1)),MAXVAL(ZQ_G(:,JLOOP1))
+ ZH_G(:) =ZQ_G(:,JLOOP1)
+ ZRV_G(:) = (XRD/XRV)*SM_FOES(ZT_G(:,JLOOP1))*0.01*ZH_G(:) &
+ /(ZPM_G(:,JLOOP1) -SM_FOES(ZT_G(:,JLOOP1))*0.01*ZH_G(:))
+ WRITE (ILUOUT0,*), ' ',MINVAL(ZRV_G(:)),MAXVAL(ZRV_G(:))
+ ZTHV_G(:)=ZT_G(:,JLOOP1) * ((XP00/ZPM_G(:,JLOOP1))**(XRD/XCPD)) * &
+ ((1. + XRV*ZRV_G(:)/XRD) / (1. + ZRV_G(:)))
+ WRITE (ILUOUT0,*), ' ',MINVAL(ZTHV_G(:)),MAXVAL(ZTHV_G(:))
+ !
+ ! Interpolation : H
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ ZH_G,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+ CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZH_LS(:,:,JLOOP1))
+ ZH_LS(:,:,JLOOP1) = MAX(MIN(ZH_LS(:,:,JLOOP1),100.),0.)
+ !
+ ! interpolation : Theta V
+ CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+ ZTHV_G,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+ CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZTHV_LS(:,:,JLOOP1))
+ !
+ END DO
+END IF
+
DEALLOCATE (ZOUT)
!
!* 2.5.5 Compute atmospheric pressure on MESO-NH grid
!
WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on MesoNH grid is being computed'
ALLOCATE (ZPF_LS(IIU,IJU,INLEVEL))
-IF (HFILE(1:3)=='ATM') THEN
- ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_LS,1,IIU),2,IJU)*XP00_LS + &
- SPREAD(SPREAD(XB_LS,1,IIU),2,IJU)*SPREAD(XPS_LS,3,INLEVEL)
-ELSE IF (HFILE=='CHEM') THEN
- ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_SV_LS,1,IIU),2,IJU)*XP00_LS + &
- SPREAD(SPREAD(XB_SV_LS,1,IIU),2,IJU)*SPREAD(XPS_SV_LS,3,INLEVEL)
-END IF
+IF (IMODEL/=10) THEN ! others than NCEP
+ IF (HFILE(1:3)=='ATM') THEN
+ ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_LS,1,IIU),2,IJU)*XP00_LS + &
+ SPREAD(SPREAD(XB_LS,1,IIU),2,IJU)*SPREAD(XPS_LS,3,INLEVEL)
+ ELSE IF (HFILE=='CHEM') THEN
+ ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_SV_LS,1,IIU),2,IJU)*XP00_LS + &
+ SPREAD(SPREAD(XB_SV_LS,1,IIU),2,IJU)*SPREAD(XPS_SV_LS,3,INLEVEL)
+ END IF
+ELSE
+ ZPF_LS(:,:,:) = 100.*SPREAD(SPREAD(IP_GFS,1,IIU),2,IJU)
+END IF
!
ALLOCATE (ZEXNF_LS(IIU,IJU,INLEVEL))
ZEXNF_LS(:,:,:) = (ZPF_LS(:,:,:)/XP00)**(XRD/XCPD)
@@ -1101,7 +1210,7 @@ END IF
!
! The way winds are processed depends upon the type of archive :
!
-! -> ECMWF
+! -> ECMWF, NCEP
! Winds are projected from a standard lat,lon grid to MesoNH grid. This correcponds to
! a rotation of an angle :
! Alpha = k.(L-L0) - Beta
@@ -1123,8 +1232,13 @@ END IF
! After this projection, the file is simil
!
IF (HFILE(1:3)=='ATM') THEN
-ITYP = 109
-ISTARTLEVEL = 1
+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))
@@ -1146,10 +1260,17 @@ SELECT CASE (IMODEL)
ISTARTLEVEL = 0
CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
END IF
+ CASE (10)
+ IPAR = 131
+ ISTARTLEVEL = 1
END SELECT
DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
- ILEV1 = JLOOP1
+ IF (IMODEL/=10) THEN ! others than NCEP
+ ILEV1 = JLOOP1
+ ELSE
+ ILEV1 = IP_GFS(JLOOP1)
+ END IF
! read component u
CALL SEARCH_FIELD(IPAR,ITYP,ILEV1,ILEV2,IGRIB,INUM)
IF (INUM < 0) THEN
@@ -1178,7 +1299,11 @@ DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
END IF
DEALLOCATE (ZVALUE)
! read component v and perform interpolation if not Arpege grid
- ILEV1 = JLOOP1
+ IF (IMODEL/=10) THEN ! others than NCEP
+ ILEV1 = JLOOP1
+ ELSE
+ ILEV1 = IP_GFS(JLOOP1)
+ END IF
CALL SEARCH_FIELD(IPAR+1,ITYP,ILEV1,ILEV2,IGRIB,INUM)
IF (INUM < 0) THEN
!callabortstop
@@ -1619,7 +1744,7 @@ INTEGER :: ILEV1 ! Level parameter 1
INTEGER :: ILEV2 ! Level parameter 2
INTEGER :: JLOOP ! Dummy counter
INTEGER :: IVERSION
-CHARACTER(LEN=20) :: YLTYPELU
+CHARACTER(LEN=24) :: YLTYPELU
CHARACTER(LEN=20) :: YLTYPE
!
! Variables used to display messages
@@ -1628,6 +1753,8 @@ INTEGER :: ILUOUT0 ! Logical unit number of the listing
ILUOUT0 = TLUOUT0%NLU
!
SELECT CASE (KLTYPE)
+CASE(100)
+ YLTYPE='isobaricInhPa'
CASE(109)
YLTYPE='hybrid'
CASE(1)
@@ -2004,6 +2131,30 @@ CASE(3,4) ! ARPEGE
PPARAM(8)=ZILOSP
PPARAM(9)=ZSTRECH
END IF
+!
+CASE(10) ! NCEP
+!
+ CALL GRIB_GET(KGRIB,'latitudeOfFirstGridPointInDegrees',ZILA1)
+ CALL GRIB_GET(KGRIB,'longitudeOfFirstGridPointInDegrees',ZILO1)
+ CALL GRIB_GET(KGRIB,'latitudeOfLastGridPointInDegrees',ZILA2)
+ CALL GRIB_GET(KGRIB,'longitudeOfLastGridPointInDegrees',ZILO2)
+ CALL GRIB_GET(KGRIB,'Nj',IINLA,IRET_GRIB)
+ CALL GRIB_GET(KGRIB,'Ni',INLO_GRIB(1),IRET_GRIB)
+ INLO_GRIB(2:)=INLO_GRIB(1)
+ KNI=IINLA*INLO_GRIB(1)
+ GREADY= (PPARAM(1)==INLO_GRIB(1) .AND. PPARAM(2)==IINLA .AND.&
+ PPARAM(3)==ZILA1 .AND. PPARAM(4)==ZILO1 .AND.&
+ PPARAM(5)==ZILA2 .AND. PPARAM(6)==ZILO2)
+ PPARAM(1)=INLO_GRIB(1)
+ PPARAM(2)=IINLA
+ PPARAM(3)=ZILA1
+ PPARAM(4)=ZILO1
+ PPARAM(5)=ZILA2
+ PPARAM(6)=ZILO2
+ IF (.NOT. GREADY) THEN
+ PLXOUT=PLONOUT
+ PLYOUT=PLATOUT
+ ENDIF
END SELECT
!JUAN
KINLO=INLO_GRIB
diff --git a/src/MNH/ver_prep_gribex_case.f90 b/src/MNH/ver_prep_gribex_case.f90
index 3abba8e8646cc39ffce4f72f7015e9162403240c..4368fc4c10aabfd7cfdaa393f936683276450f38 100644
--- a/src/MNH/ver_prep_gribex_case.f90
+++ b/src/MNH/ver_prep_gribex_case.f90
@@ -1,8 +1,13 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision
+! MASDEV4_7 prep_real 2006/05/23 14:49:51
+!-----------------------------------------------------------------
! ################################
MODULE MODI_VER_PREP_GRIBEX_CASE
! ################################
@@ -83,6 +88,7 @@ END MODULE MODI_VER_PREP_GRIBEX_CASE
!! Nov, 22 2000 (I. Mallet) add scalar variables
!! Nov, 22 2000 (P. Jabouille) change routine name
!! 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
!-------------------------------------------------------------------------------
!
@@ -259,6 +265,8 @@ IF (HFILE(1:3)=='ATM') THEN
END IF
!
IF (HFILE(1:3)=='ATM') THEN
+
+ IF (SIZE(XB_LS)/=0) THEN ! hybrid level (w at flux points)
CALL VER_INTERP_TO_MIXED_GRID('ATM ',.TRUE.,XZS_LS,XZSMT_LS, &
ZZMASS_LS,ZSV_LS, &
ZZFLUX_LS,XPS_LS,ZPMHP_LS, &
@@ -267,6 +275,16 @@ IF (HFILE(1:3)=='ATM') THEN
ZTKE_LS, &
ZU_LS,ZV_LS, &
ZW_LS,'FLUX' )
+ ELSE ! isobaric surfaces (w at mass points)
+ CALL VER_INTERP_TO_MIXED_GRID('ATM ',.TRUE.,XZS_LS,XZSMT_LS, &
+ ZZMASS_LS,ZSV_LS, &
+ ZZFLUX_LS,XPS_LS,ZPMHP_LS, &
+ ZTHV_LS,ZR_LS, &
+ ZHU_LS, &
+ ZTKE_LS, &
+ ZU_LS,ZV_LS, &
+ ZW_LS,'MASS' )
+ END IF
ELSE IF (HFILE=='CHEM') THEN
CALL VER_INTERP_TO_MIXED_GRID(HFILE,.TRUE.,XZS_SV_LS,XZS_SV_LS,&
ZZMASS_LS,ZSV_LS )
diff --git a/src/SURFEX/modd_grid_grib.F90 b/src/SURFEX/modd_grid_grib.F90
index 9ae078577cfa550f3be06be9e9ec9ac8ad93e753..0a653fb366c1a9a594bb7652c2f1a20a6dde9042 100644
--- a/src/SURFEX/modd_grid_grib.F90
+++ b/src/SURFEX/modd_grid_grib.F90
@@ -36,8 +36,10 @@ USE GRIB_API, ONLY : kindOfInt
IMPLICIT NONE
!
INTEGER :: NNI ! total number of physical points
+INTEGER :: NGRIB_VERSION ! GRIB-API version (1 or 2)
+
!
- CHARACTER(LEN=6) :: CINMODEL!
+ CHARACTER(LEN=6) :: CINMODEL !
CHARACTER(LEN=28) :: CGRIB_FILE
INTEGER(KIND=kindOfInt) :: NIDX
!
diff --git a/src/SURFEX/mode_read_grib.F90 b/src/SURFEX/mode_read_grib.F90
index 7594cbfc1e5d516243e7e3fcac82a8896e8c977d..2eb32506329b395fc5f067472fdda84b334668f5 100644
--- a/src/SURFEX/mode_read_grib.F90
+++ b/src/SURFEX/mode_read_grib.F90
@@ -19,7 +19,7 @@ CONTAINS
SUBROUTINE MAKE_GRIB_INDEX(HGRIB)
! ####################
!
-USE MODD_GRID_GRIB, ONLY : CGRIB_FILE, NIDX
+USE MODD_GRID_GRIB, ONLY : CGRIB_FILE, NIDX, NGRIB_VERSION
!
IMPLICIT NONE
!
@@ -27,6 +27,8 @@ IMPLICIT NONE
!
INTEGER(KIND=kindOfInt) :: IRET
REAL(KIND=JPRB) :: ZHOOK_HANDLE
+INTEGER, DIMENSION(:), ALLOCATABLE :: IVERSION
+INTEGER :: ILEN
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:MAKE_GRIB_INDEX',0,ZHOOK_HANDLE)
!
@@ -35,7 +37,25 @@ IF (CGRIB_FILE==HGRIB) RETURN
!
CGRIB_FILE=HGRIB
!
+ CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'editionNumber',IRET)
+IF (IRET/=0) THEN
+ CALL ABOR1_SFX("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while searching edition number")
+ELSE
+ call GRIB_INDEX_GET_SIZE(NIDX,'editionNumber',ILEN)
+ ALLOCATE(IVERSION(ILEN))
+ CALL GRIB_INDEX_GET(NIDX, 'editionNumber', IVERSION , IRET)
+ NGRIB_VERSION=IVERSION(1)
+
+ CALL GRIB_INDEX_RELEASE(NIDX,IRET)
+ IF (IRET/=0) CALL ABOR1_SFX("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while deleting the grib index")
+ENDIF
+
+!test version
+IF (NGRIB_VERSION==1) THEN
CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'indicatorOfParameter',IRET)
+ELSEIF (NGRIB_VERSION==2) THEN
+ CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'paramId',IRET)
+ENDIF
IF (IRET/=0) CALL ABOR1_SFX("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while creating the grib index")
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:MAKE_GRIB_INDEX',1,ZHOOK_HANDLE)
@@ -67,7 +87,7 @@ IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:CLEAR_GRIB_INDEX',1,ZHOOK_HANDLE)
END SUBROUTINE CLEAR_GRIB_INDEX
!-------------------------------------------------------------------
! ####################
- SUBROUTINE GET_GRIB_MESSAGE(KLUOUT,KLTYPE,KLEV1,KLEV2,KGRIB,KFOUND)
+ SUBROUTINE GET_GRIB_MESSAGE(KLUOUT,KLTYPE,KLEV1,KLEV2,KGRIB,KFOUND,HTYPELEVEL,PLEV1,PLEV2)
! ####################
!
USE MODD_GRID_GRIB, ONLY : NIDX
@@ -80,10 +100,16 @@ INTEGER, INTENT(INOUT) :: KLEV1
INTEGER, INTENT(INOUT) :: KLEV2
INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KGRIB
INTEGER, INTENT(OUT) :: KFOUND
+CHARACTER(LEN=*), INTENT(INOUT), OPTIONAL :: HTYPELEVEL ! TypeOfLevel JPMODIF
+REAL, INTENT(INOUT), OPTIONAL :: PLEV1 ! top level of soil
+REAL, INTENT(INOUT), OPTIONAL :: PLEV2 ! Bottom level of soil
!
+
INTEGER :: ILTYPE
INTEGER :: ILEV1
INTEGER :: ILEV2
+CHARACTER(LEN=50) :: YTYPELEVEL ! TypeOfLevel JPMODIF
+REAL :: ZLEV1,ZLEV2
INTEGER(KIND=kindOfInt) :: IRET
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -101,6 +127,11 @@ DO WHILE (IRET /= GRIB_END_OF_INDEX .AND. KFOUND/=3)
IF (KLTYPE/=-2) THEN
CALL GRIB_GET(KGRIB,'indicatorOfTypeOfLevel',ILTYPE,IRET)
CALL TEST_IRET(KLUOUT,ILTYPE,KLTYPE,IRET)
+ ELSE
+ IF (PRESENT(HTYPELEVEL)) THEN
+ CALL GRIB_GET(KGRIB,'typeOfLevel',YTYPELEVEL,IRET)
+ CALL TEST_IRET_STR(KLUOUT,TRIM(YTYPELEVEL),HTYPELEVEL,IRET)
+ ENDIF
ENDIF
!
IF (IRET.EQ.0) THEN
@@ -110,6 +141,10 @@ DO WHILE (IRET /= GRIB_END_OF_INDEX .AND. KFOUND/=3)
IF (KLEV1/=-2) THEN
CALL GRIB_GET(KGRIB,'topLevel',ILEV1,IRET)
CALL TEST_IRET(KLUOUT,ILEV1,KLEV1,IRET)
+ ELSE IF (PRESENT(PLEV1)) THEN !JP
+ CALL GRIB_GET(KGRIB,'topLevel',ZLEV1,IRET)
+ CALL TEST_IRET_FLOAT(KLUOUT,ZLEV1,PLEV1,IRET)
+
ENDIF
!
IF (IRET.EQ.0) THEN
@@ -119,6 +154,9 @@ DO WHILE (IRET /= GRIB_END_OF_INDEX .AND. KFOUND/=3)
IF (KLEV2/=-2) THEN
CALL GRIB_GET(KGRIB,'bottomLevel',ILEV2,IRET)
CALL TEST_IRET(KLUOUT,ILEV2,KLEV2,IRET)
+ ELSE IF (PRESENT(PLEV2)) THEN !JP
+ CALL GRIB_GET(KGRIB,'bottomLevel',ZLEV2,IRET)
+ CALL TEST_IRET_FLOAT(KLUOUT,ZLEV2,PLEV2,IRET)
ENDIF
!
IF (IRET.EQ.0) KFOUND = KFOUND + 1
@@ -167,14 +205,76 @@ ENDIF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',1,ZHOOK_HANDLE)
END SUBROUTINE TEST_IRET
+
+! ##############
+ SUBROUTINE TEST_IRET_STR(KLUOUT,VAL1,VAL0,KRET)
+! ##############
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
+CHARACTER(LEN=*), INTENT(IN) :: VAL1
+CHARACTER(LEN=*), INTENT(INOUT) :: VAL0
+INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KRET ! number of the message researched
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',0,ZHOOK_HANDLE)
+!
+IF (KRET > 0) THEN
+ WRITE (KLUOUT,'(A)')' | Error encountered in the Grib file, skipping field'
+ELSE IF (KRET == -6) THEN
+ WRITE (KLUOUT,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
+ELSEIF (VAL1 /= VAL0) THEN
+ IF (VAL0 == '') THEN
+ VAL0 = VAL1
+ ELSE
+ KRET=1
+ ENDIF
+ENDIF
+!
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',1,ZHOOK_HANDLE)
+END SUBROUTINE TEST_IRET_STR
+
+! ##############
+ SUBROUTINE TEST_IRET_FLOAT(KLUOUT,VAL1,VAL0,KRET)
+! ##############
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
+REAL, INTENT(IN) :: VAL1
+REAL, INTENT(INOUT) :: VAL0
+INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KRET ! number of the message researched
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',0,ZHOOK_HANDLE)
+!
+IF (KRET > 0) THEN
+ WRITE (KLUOUT,'(A)')' | Error encountered in the Grib file, skipping field'
+ELSE IF (KRET == -6) THEN
+ WRITE (KLUOUT,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
+ELSEIF (VAL1 /= VAL0) THEN
+ IF (VAL0 == -1.0) THEN
+ VAL0 = VAL1
+ ELSE
+ KRET=1
+ ENDIF
+ENDIF
+!
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:TEST_IRET',1,ZHOOK_HANDLE)
+END SUBROUTINE TEST_IRET_FLOAT
+
!
END SUBROUTINE GET_GRIB_MESSAGE
!-------------------------------------------------------------------
! ####################
- SUBROUTINE READ_GRIB(HGRIB,KLUOUT,KPARAM,KRET,PFIELD,KLTYPE,KLEV1,KLEV2,KPARAM2)
+ SUBROUTINE READ_GRIB(HGRIB,KLUOUT,KPARAM,KRET,PFIELD,KLTYPE,KLEV1,KLEV2,KPARAM2, &
+ HTYPELEVEL,PLEV1,PLEV2)
! ####################
!
-USE MODD_GRID_GRIB, ONLY : NIDX
+USE MODD_GRID_GRIB, ONLY : NIDX, NGRIB_VERSION
!
IMPLICIT NONE
!
@@ -187,11 +287,15 @@ 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
+CHARACTER(LEN=*), INTENT(INOUT), OPTIONAL :: HTYPELEVEL
+!
+REAL, INTENT(INOUT), OPTIONAL :: PLEV1,PLEV2
!
INTEGER :: ILTYPE, ILEV1, ILEV2
INTEGER(KIND=kindOfInt) :: IGRIB
INTEGER :: ISIZE, IFOUND
REAL(KIND=JPRB) :: ZHOOK_HANDLE
+REAL :: ZLEV1,ZLEV2
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB',0,ZHOOK_HANDLE)
!
@@ -201,16 +305,36 @@ ILEV1=-2
IF (PRESENT(KLEV1)) ILEV1=KLEV1
ILEV2=-2
IF (PRESENT(KLEV2)) ILEV2=KLEV2
+ZLEV1=-2.0
+IF (PRESENT(PLEV1)) ZLEV1=PLEV1
+ZLEV2=-2.0
+IF (PRESENT(PLEV2)) ZLEV2=PLEV2
!
CALL MAKE_GRIB_INDEX(HGRIB)
!
IFOUND=0
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)
+END IF
CALL GRIB_NEW_FROM_INDEX(NIDX,IGRIB,KRET)
-IF (KRET.EQ.0) CALL GET_GRIB_MESSAGE(KLUOUT,ILTYPE,ILEV1,ILEV2,IGRIB,IFOUND)
+
+
+!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
+ CALL GET_GRIB_MESSAGE(KLUOUT,ILTYPE,ILEV1,ILEV2,IGRIB,IFOUND)
+ ENDIF
+ENDIF
!
+
+!WRITE (KLUOUT,*) 'READ_GRIB GRIB_NEW_FROM_INDEX ',KPARAM,IGRIB,KRET,IFOUND
+
IF (PRESENT(KPARAM2)) THEN
IF (IFOUND/=3) THEN
CALL GRIB_INDEX_SELECT(NIDX,'indicatorOfParameter',KPARAM2,KRET)
@@ -225,12 +349,15 @@ IF (PRESENT(KPARAM2)) THEN
ENDIF
ENDIF
!
+
+
IF (IFOUND==3) THEN
!
IF (PRESENT(KLTYPE)) KLTYPE = ILTYPE
IF (PRESENT(KLEV1)) KLEV1 = ILEV1
IF (PRESENT(KLEV2)) KLEV2 = ILEV2
!
+
IF (.NOT.ASSOCIATED(PFIELD)) THEN
CALL GRIB_GET_SIZE(IGRIB,'values',ISIZE,KRET)
IF (KRET.NE.0) CALL ABOR1_SFX("MODE_READ_GRIB:READ_GRIB: Problem getting size of values")
@@ -267,6 +394,8 @@ REAL, DIMENSION(:), POINTER :: PMASK ! Land mask
INTEGER(KIND=kindOfInt) :: IRET ! return code
INTEGER :: ILTYPE ! leveltype
INTEGER :: ILEV ! level
+INTEGER :: INUM_ZS,ISIZE,ICOUNT,JLOOP,IPARAM,IGRIB
+CHARACTER(LEN=24) :: YLTYPELU
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_LAND_MASK',0,ZHOOK_HANDLE)
@@ -275,6 +404,10 @@ WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_LAND_MASK: | Reading land mask fr
SELECT CASE (HINMODEL)
CASE ('ECMWF ')
CALL READ_GRIB(HGRIB,KLUOUT,172,IRET,PMASK)
+ CASE ('NCEP ')
+ CALL READ_GRIB(HGRIB,KLUOUT,172,IRET,PMASK)
+
+
CASE ('ARPEGE','ALADIN','MOCAGE')
CALL READ_GRIB(HGRIB,KLUOUT,81,IRET,PMASK)
CASE ('HIRLAM')
@@ -312,6 +445,10 @@ INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
REAL, DIMENSION(:), POINTER :: PZS !
!
INTEGER(KIND=kindOfInt) :: IRET ! return code
+
+INTEGER :: ISIZE,JLOOP,ICOUNT
+CHARACTER(LEN=24) :: YLTYPELU
+
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------
!* Read orography
@@ -321,6 +458,8 @@ 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)
+ CASE ('NCEP ')
+ CALL READ_GRIB(HGRIB,KLUOUT,228002,IRET,PZS)
CASE ('ARPEGE','MOCAGE')
CALL READ_GRIB(HGRIB,KLUOUT,8,IRET,PZS)
CASE ('HIRLAM','ALADIN')
@@ -402,6 +541,8 @@ INTEGER(KIND=kindOfInt) :: IRET ! return code
INTEGER :: ILTYPE ! type of level (Grib code table 3)
INTEGER :: ILEV ! level definition
REAL(KIND=JPRB) :: ZHOOK_HANDLE
+CHARACTER(LEN=7) :: YTYPELEVEL ! Type of searched level
+
!-------------------------------------------------------------------
!* Read surface temperature
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_T',0,ZHOOK_HANDLE)
@@ -430,6 +571,12 @@ SELECT CASE (HINMODEL)
ILEV=904
CALL READ_GRIB(HGRIB,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)
+
CASE DEFAULT
CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_T:OPTION NOT SUPPORTED '//HINMODEL)
END SELECT
@@ -489,7 +636,7 @@ SELECT CASE (HINMODEL)
CASE ('ECMWF ')
CALL READ_GRIB(HGRIB,KLUOUT,34,IRET,PSST)
IF (IRET /= 0) CALL READ_GRIB_T(HGRIB,KLUOUT,HINMODEL,PSST)
- CASE ('ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_T(HGRIB,KLUOUT,HINMODEL,PSST)
CASE DEFAULT
CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_SST:OPTION NOT SUPPORTED '//HINMODEL)
@@ -523,13 +670,13 @@ SELECT CASE (HINMODEL)
CASE ('ECMWF ')
CALL READ_GRIB(HGRIB,KLUOUT,3080,IRET,PTS)
IF (IRET /= 0) CALL READ_GRIB_T2(HGRIB,KLUOUT,HINMODEL,PMASK,PTS)
- CASE ('ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_T2(HGRIB,KLUOUT,HINMODEL,PMASK,PTS)
CASE DEFAULT
CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_TSWATER:OPTION NOT SUPPORTED '//HINMODEL)
END SELECT
!
-IF (SIZE(PMASK)==SIZE(PTS)) WHERE (PMASK(:)/=0.) PTS = XUNDEF
+IF (SIZE(PMASK)==SIZE(PTS)) WHERE (PMASK(:)/=1.) PTS = XUNDEF
!
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TSWATER',1,ZHOOK_HANDLE)
END SUBROUTINE READ_GRIB_TSWATER
@@ -541,7 +688,7 @@ END SUBROUTINE READ_GRIB_TSWATER
USE MODD_SURF_PAR, ONLY : XUNDEF
!
IMPLICIT NONE
-!
+ !
CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
@@ -551,6 +698,8 @@ REAL, DIMENSION(:), POINTER :: PT2 !
INTEGER(KIND=kindOfInt) :: IRET
INTEGER :: ILTYPE, ILEV ! type of level (Grib code table 3)
REAL(KIND=JPRB) :: ZHOOK_HANDLE
+CHARACTER(LEN=19) :: YTYPELEVEL ! Type of searched level
+REAL :: ZLEV1,ZLEV2
!-------------------------------------------------------------------
!* Read deep soil temperature
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_T2',0,ZHOOK_HANDLE)
@@ -565,7 +714,13 @@ SELECT CASE (HINMODEL)
IF (IRET /= 0) THEN
ILTYPE=105
CALL READ_GRIB(HGRIB,KLUOUT,11,IRET,PT2,KLTYPE=ILTYPE)
- ENDIF
+ 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)
CASE ('HIRLAM ')
ILTYPE=105
ILEV=954
@@ -632,7 +787,7 @@ IF (KLTYPE==112) THEN
PD = (KLEV2 - KLEV1) / 100.
ELSE
IF (KNLAYERDEEP == 4) THEN
- PD = PV4
+ PD = PV4
ELSE
PD = PV
END IF
@@ -713,6 +868,8 @@ INTEGER :: INLAYERDEEP! number of deep moisture layers
REAL, DIMENSION(:), POINTER :: ZFIELD => NULL() ! first layer temperature
REAL, DIMENSION(:,:), ALLOCATABLE:: ZTG ! first layer temperature
REAL, DIMENSION(:) , ALLOCATABLE:: ZD
+
+
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!--------------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TG_ECMWF',0,ZHOOK_HANDLE)
@@ -1882,7 +2039,9 @@ SELECT CASE(HINMODEL)
CASE('ECMWF ')
CALL READ_GRIB(HGRIB,KLUOUT,141,IRET,ZFIELD)
CASE('ARPEGE','ALADIN','MOCAGE','HIRLAM')
- CALL READ_GRIB(HGRIB,KLUOUT,66,IRET,ZFIELD)
+ CALL READ_GRIB(HGRIB,KLUOUT,66,IRET,ZFIELD)
+ CASE('NCEP ')
+ CALL READ_GRIB(HGRIB,KLUOUT,3066,IRET,ZFIELD)
CASE DEFAULT
CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_SNOW_VEG_AND_DEPTH: OPTION NOT SUPPORTED '//HINMODEL)
END SELECT
@@ -2128,4 +2287,344 @@ IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TF_TEB',1,ZHOOK_HANDLE)
!------------------------------------------------------------------------------
END SUBROUTINE READ_GRIB_TF_TEB
!-------------------------------------------------------------------
+! #######################
+ SUBROUTINE READ_GRIB_TG_NCEP(HGRIB,KLUOUT,HINMODEL,PMASK,PTG,PD)
+! #######################
+!
+IMPLICIT NONE
+!
+!* dummy arguments
+! ---------------
+ CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
+INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
+ CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
+REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
+REAL, DIMENSION(:,:), POINTER :: PTG ! field to initialize
+REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
+!
+!* local variables
+! ---------------
+INTEGER(KIND=kindOfInt) :: IRET ! return code
+INTEGER :: ILTYPE ! type of level (Grib code table 3)
+INTEGER :: ILEV1 ! level definition
+INTEGER :: ILEV2 ! level definition
+INTEGER :: JL ! layer loop counter
+INTEGER :: INLAYERDEEP! number of deep moisture layers
+REAL, DIMENSION(:), POINTER :: ZFIELD => NULL() ! first layer temperature
+REAL, DIMENSION(:,:), ALLOCATABLE:: ZTG ! first layer temperature
+REAL, DIMENSION(:) , ALLOCATABLE:: ZD
+REAL :: ZLEV1,ZLEV2 ! level definition in float
+CHARACTER(LEN=19) :: YTYPELEVEL ! Type of searched level
+
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!--------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TG_NCEP',0,ZHOOK_HANDLE)
+WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_TG_NCEP: | Reading soil temperature'
+!
+ALLOCATE(ZD(5))
+!
+! 1. Search and read level 1 (and its depth)
+! --------------------------------------
+
+YTYPELEVEL='depthBelowLandLayer'
+
+ILTYPE= -2
+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)
+!
+IF (IRET== 0) THEN
+ ZD(1)=(ZLEV2 - ZLEV1)
+
+ ALLOCATE(ZTG(SIZE(ZFIELD),5))
+ ZTG(:,1)=ZFIELD
+
+ELSE
+ CALL ABOR1_SFX('MODE_READ_GRIB: SOIL TEMPERATURE LEVEL 1 MISSING (READ_GRIB_TG_NCEP)')
+ENDIF
+!
+! 2. Search and read level 4 (and its depth) This level is optionnal
+! ---------------------------------------------------------------
+ILTYPE= -2
+ILEV1 = -2
+ILEV2 = -2
+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)
+
+IF (IRET == 0) THEN
+ INLAYERDEEP = 4
+ ZD(4)=(ZLEV2 - ZLEV1)
+
+ ZTG(:,4)=ZFIELD
+ELSE
+ INLAYERDEEP = 3
+ ZD(4) = 0.
+ENDIF
+!
+! 3. Search and read level 3 (and its depth) This level is optionnal
+! ---------------------------------------------------------------
+! ---------------------------------------------------------------
+ILTYPE= -2
+ILEV1 = -2
+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)
+!
+IF (IRET == 0) THEN
+ ZD(3)=(ZLEV2 - ZLEV1)
+ ZTG(:,3)=ZFIELD
+ELSE
+ INLAYERDEEP = 2
+ ZD(3) = 0.
+ENDIF
+!
+! 4. Search and read level 2 (and its depth)
+! ---------------------------------------
+ILTYPE= -2
+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)
+!
+IF (IRET== 0) THEN
+ ZD(2)=(ZLEV2 - ZLEV1)
+ ZTG(:,2)=ZFIELD
+ DEALLOCATE(ZFIELD)
+ELSE
+ CALL ABOR1_SFX('MODE_READ_GRIB: SOIL TEMPERATURE LEVEL 2 MISSING (READ_GRIB_TG_NCEP)')
+ENDIF
+!--------------------------------------------------------------------------------
+! 5. Assumes uniform temperature profile up to 3m depth
+! -------------------------------------------------
+!
+WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:UP TO 3m'
+ write(KLUOUT,'(a,i3,i3)') 'editionNumber MAKE_GRIB_INDEX',INLAYERDEEP
+
+IF(SUM(ZD(1:INLAYERDEEP)) < 3.) THEN
+ !We add a temperature layer
+ INLAYERDEEP=INLAYERDEEP+1
+ write(KLUOUT,'(a,i3,i3)') 'editionNumber MAKE_GRIB_INDEX',INLAYERDEEP
+
+ ZD(INLAYERDEEP)=3.-SUM(ZD(1:INLAYERDEEP-1))
+ write(KLUOUT,*) 'ZD',ZD
+
+ 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)
+!
+!write(KLUOUT,*) 'OUT FILL ZD',PD,PTG
+
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TG_NCEP',1,ZHOOK_HANDLE)
+
+WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:END READ GRIB'
+
+END SUBROUTINE READ_GRIB_TG_NCEP
+!--------------------------------------------------------------------------------
+! #######################
+ SUBROUTINE READ_GRIB_WG_NCEP(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
+! #######################
+!
+USE MODD_GRID_GRIB, ONLY : NNI
+USE MODD_SURF_PAR, ONLY : XUNDEF
+!
+IMPLICIT NONE
+!
+!* dummy arguments
+! ---------------
+ CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
+INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
+ CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
+REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
+REAL, DIMENSION(:,:), POINTER :: PFIELD ! field to initialize
+REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
+!* local variables
+! ---------------
+INTEGER(KIND=kindOfInt) :: IRET ! return code
+INTEGER :: ILTYPE ! type of level (Grib code table 3)
+INTEGER :: ILEV1 ! level definition
+INTEGER :: ILEV2 ! level definition
+REAL, DIMENSION(:), POINTER :: ZFIELD => NULL()
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZWG ! first water reservoir
+REAL, DIMENSION(:), ALLOCATABLE :: ZD ! Height of each layer
+INTEGER :: INLAYERDEEP! number of deep moisture layers
+REAL :: ZWWILT ! ECMWF wilting point
+REAL :: ZWFC ! ECMWF field capacity
+REAL :: ZWSAT ! ECMWF saturation
+INTEGER :: JL ! loop counter on layers
+REAL :: ZLEV1,ZLEV2 ! level definition in float
+CHARACTER(LEN=19) :: YTYPELEVEL ! Type of searched level
+
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!--------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_WG_NCEP',0,ZHOOK_HANDLE)
+WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_WG_NCEP: | Reading soil moisture'
+WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_WG_NCEP: | WARNING READING LOW VEGETATION TILE (NR 4) ONLY'
+!
+ALLOCATE(ZD(5))
+
+! 1. Search and read level 1 (and its depth)
+! --------------------------------------
+
+YTYPELEVEL='depthBelowLandLayer'
+
+ILTYPE= -2
+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)
+!
+IF (IRET== 0) THEN
+ ZD(1)=(ZLEV2 - ZLEV1)
+
+ ALLOCATE(ZWG(SIZE(ZFIELD),5))
+ ZWG(:,1)=ZFIELD
+ELSE
+ CALL ABOR1_SFX('MODE_READ_GRIB: SOIL MOISTURE LEVEL 1 MISSING (READ_GRIB_TG_NCEP)')
+ENDIF
+!
+! 2. Search and read level 4 (and its depth) This level is optionnal
+! ---------------------------------------------------------------
+ILTYPE= -2
+ILEV1 = -2
+ILEV2 = -2
+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)
+
+IF (IRET == 0) THEN
+ INLAYERDEEP = 4
+ ZD(4)=(ZLEV2 - ZLEV1)
+
+ ZWG(:,4)=ZFIELD
+ELSE
+ INLAYERDEEP = 3
+ ZD(4) = 0.
+ENDIF
+!
+! 3. Search and read level 3 (and its depth) This level is optionnal
+! ---------------------------------------------------------------
+! ---------------------------------------------------------------
+ILTYPE= -2
+ILEV1 = -2
+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)
+!
+IF (IRET == 0) THEN
+ ZD(3)=(ZLEV2 - ZLEV1)
+ ZWG(:,3)=ZFIELD
+ELSE
+ INLAYERDEEP = 2
+ ZD(3) = 0.
+ENDIF
+!
+! 4. Search and read level 2 (and its depth)
+! ---------------------------------------
+ILTYPE= -2
+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)
+!
+IF (IRET== 0) THEN
+ ZD(2)=(ZLEV2 - ZLEV1)
+ ZWG(:,2)=ZFIELD
+ DEALLOCATE(ZFIELD)
+ELSE
+ CALL ABOR1_SFX('MODE_READ_GRIB: SOIL MOISTURE LEVEL 2 MISSING (READ_GRIB_TG_NCEP)')
+ENDIF
+!--------------------------------------------------------------------------------
+! 5. Assumes uniform temperature profile up to 3m depth
+! -------------------------------------------------
+!
+IF(SUM(ZD(1:INLAYERDEEP)) < 3.) THEN
+ !We add a humidity layer
+ INLAYERDEEP=INLAYERDEEP+1
+ ZD(INLAYERDEEP)=3.-SUM(ZD(1:INLAYERDEEP-1))
+ ZWG(:,INLAYERDEEP)=ZWG(:,INLAYERDEEP-1)
+ENDIF
+!
+!--------------------------------------------------------------------------------
+! 6. Set temperature profile and layer thicknesses
+! ----------------------------------------------
+ CALL FILL_PFIELD(KLUOUT,'READ_GRIB_WG_NCEP',INLAYERDEEP,ZD,ZWG,PMASK,PFIELD,PD)
+DEALLOCATE(ZD)
+DEALLOCATE(ZWG)
+
+!--------------------------------------------------------------------------------
+! 7. Convert from specific humidity to relative humidity
+! ---------------------------------------------------
+! Compute model's constants
+ZWFC = 0.171
+ZWWILT = 0.086
+!
+! Then perform conversion
+DO JL=1,INLAYERDEEP
+ WHERE (PFIELD(:,JL).NE.XUNDEF) PFIELD(:,JL) = (PFIELD(:,JL) - ZWWILT) / (ZWFC - ZWWILT)
+ENDDO
+
+!
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_WG_NCEP',1,ZHOOK_HANDLE)
+END SUBROUTINE READ_GRIB_WG_NCEP
+!----------------------------------------------------------------------------
+! #######################
+ SUBROUTINE READ_GRIB_WGI_NCEP(HGRIB,KLUOUT,PFIELD,PD)
+! #######################
+!
+USE MODD_GRID_GRIB, ONLY : NNI
+!
+IMPLICIT NONE
+!
+!* dummy arguments
+! ---------------
+ CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
+INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
+REAL, DIMENSION(:,:), POINTER :: PFIELD ! field to initialize
+REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
+!
+!* local variables
+! ---------------
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!--------------------------------------------------------------------------------
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_WGI_NCEP',0,ZHOOK_HANDLE)
+!
+ALLOCATE (PFIELD(NNI,4))
+ALLOCATE (PD (NNI,4))
+PFIELD(:,:) = 0.
+!
+PD (:,1) = 0.1
+PD (:,2) = 0.3
+PD (:,3) = 0.6
+PD (:,4) = 1.0
+
+
+!
+IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_WGI_NCEP',1,ZHOOK_HANDLE)
+END SUBROUTINE READ_GRIB_WGI_NCEP
+!-------------------------------------------------------------------
END MODULE MODE_READ_GRIB
diff --git a/src/SURFEX/prep_flake_grib.F90 b/src/SURFEX/prep_flake_grib.F90
index 30989a93c8c9b586d145f64dc8ebe6c404d7c6ff..1e09421b3879b452d2204baa216a317ac1db8597 100644
--- a/src/SURFEX/prep_flake_grib.F90
+++ b/src/SURFEX/prep_flake_grib.F90
@@ -81,7 +81,7 @@ SELECT CASE(HSURF)
!
CASE('ZS ')
SELECT CASE (YINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_ZS_LAND(HFILE,KLUOUT,YINMODEL,ZMASK,ZFIELD)
ALLOCATE(PFIELD(SIZE(ZFIELD),1))
PFIELD(:,1) = ZFIELD(:)
@@ -94,7 +94,7 @@ SELECT CASE(HSURF)
!
CASE('TS ')
SELECT CASE (YINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_TSWATER(HFILE,KLUOUT,YINMODEL,ZMASK,ZFIELD)
ALLOCATE(PFIELD(SIZE(ZFIELD),1))
PFIELD(:,1) = ZFIELD(:)
diff --git a/src/SURFEX/prep_grib_grid.F90 b/src/SURFEX/prep_grib_grid.F90
index e33be8e3a6078b50d0b609fb8b7260ac0fef1501..28868257930e1e6011949f1ecca4cbb282fe3bd7 100644
--- a/src/SURFEX/prep_grib_grid.F90
+++ b/src/SURFEX/prep_grib_grid.F90
@@ -114,6 +114,11 @@ INTEGER :: JLOOP1 ! Dummy counter
INTEGER :: INFOMPI, J
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
+INTEGER :: JLOOP
+INTEGER(KIND=kindOfInt) :: ICOUNT ! number of messages in the file
+INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE :: IALLGRIB ! number of the grib in memory
+!
+!
!---------------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('PREP_GRIB_GRID_1',0,ZHOOK_HANDLE)
@@ -135,6 +140,23 @@ END IF
IF (IRET /= 0) THEN
CALL ABOR1_SFX('PREP_GRIB_GRID: Error in reading the grib file')
END IF
+!*JPC*GFS*AUG2015
+ CALL GRIB_COUNT_IN_FILE(IUNIT,ICOUNT)
+IF (IRET /= 0) THEN
+ CALL ABOR1_SFX('PREP_GRIB_GRID: Error in reading the grib file')
+END IF
+ALLOCATE(IALLGRIB(ICOUNT))
+! initialize the tabular with a negativ number
+! ( all the IALLGRIB will be different )
+IALLGRIB(:)=-12
+!charge all the message in memory
+DO JLOOP=1,ICOUNT
+CALL GRIB_NEW_FROM_FILE(IUNIT,IALLGRIB(JLOOP),IRET)
+IF (IRET /= 0) THEN
+ CALL ABOR1_SFX('PREP_GRIB_GRID: Error in reading the grib file')
+END IF
+END DO
+!*JPC*GFS*AUG2015
!
! close the grib file
CALL GRIB_CLOSE_FILE(IUNIT)
@@ -156,6 +178,10 @@ END IF
!
SELECT CASE (ICENTER)
+ CASE (7)
+ WRITE (KLUOUT,'(A)') ' | Grib file from National Centres for Environmental Prediction'
+ HINMODEL = 'NCEP '
+ HGRIDTYPE= 'LATLON '
CASE (96)
SELECT CASE (HGRID)
diff --git a/src/SURFEX/prep_isba_grib.F90 b/src/SURFEX/prep_isba_grib.F90
index 5b2430a6e402b66d451c4e951fff63574274672d..f21f91208b3ac31fd80b7418aec91cbf4ba012e6 100644
--- a/src/SURFEX/prep_isba_grib.F90
+++ b/src/SURFEX/prep_isba_grib.F90
@@ -97,6 +97,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_TG_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_TG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CASE('NCEP ')
+ CALL READ_GRIB_TG_NCEP(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
@@ -110,6 +112,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_WG_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_WG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CASE('NCEP ')
+ CALL READ_GRIB_WG_NCEP(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
@@ -125,6 +129,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_WGI_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_WGI_HIRLAM(HFILE,KLUOUT,ZFIELD,ZD)
+ CASE('NCEP ')
+ CALL READ_GRIB_WGI_NCEP(HFILE,KLUOUT,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
!
@@ -175,6 +181,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_TG_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_TG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CASE('NCEP ')
+ CALL READ_GRIB_TG_NCEP(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
END SELECT
ALLOCATE(PFIELD(NNI,1,1))
PFIELD(:,1,1) =ZFIELD(:,1)
diff --git a/src/SURFEX/prep_seaflux_grib.F90 b/src/SURFEX/prep_seaflux_grib.F90
index 2267fcabe72db7622ff0cea0bb7761198b55e03b..f4ee510ad67f87f25be2a57db850b76a59782361 100644
--- a/src/SURFEX/prep_seaflux_grib.F90
+++ b/src/SURFEX/prep_seaflux_grib.F90
@@ -75,7 +75,7 @@ SELECT CASE(HSURF)
!
CASE('ZS ')
SELECT CASE (CINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_ZS_SEA(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD)
ALLOCATE(PFIELD(SIZE(ZFIELD),1))
PFIELD(:,1) = ZFIELD(:)
@@ -88,7 +88,7 @@ SELECT CASE(HSURF)
!
CASE('SST ')
SELECT CASE (CINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_SST(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD)
ALLOCATE(PFIELD(SIZE(ZFIELD),1))
PFIELD(:,1) = ZFIELD(:)
diff --git a/src/SURFEX/prep_teb_garden_grib.F90 b/src/SURFEX/prep_teb_garden_grib.F90
index 49cf5861c407a0c7ef4dc50ae6de57ad58b35649..40dff4d0659c57a570589df5e96c10ba5bc2d45a 100644
--- a/src/SURFEX/prep_teb_garden_grib.F90
+++ b/src/SURFEX/prep_teb_garden_grib.F90
@@ -95,6 +95,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_TG_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_TG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CASE('NCEP')
+ CALL READ_GRIB_TG_NCEP(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
@@ -107,6 +109,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_WG_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_WG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CASE('NCEP')
+ CALL READ_GRIB_WG_NCEP(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
@@ -122,6 +126,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_WGI_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_WGI_HIRLAM(HFILE,KLUOUT,ZFIELD,ZD)
+ CASE('NCEP')
+ CALL READ_GRIB_WGI_NCEP(HFILE,KLUOUT,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
!
diff --git a/src/SURFEX/prep_teb_greenroof_grib.F90 b/src/SURFEX/prep_teb_greenroof_grib.F90
index 4bbd19b523c875e2bf8cd9c5fd8e5405e7326cd0..abf4f19298774cecc3ed6e6a09cc075994674ca5 100644
--- a/src/SURFEX/prep_teb_greenroof_grib.F90
+++ b/src/SURFEX/prep_teb_greenroof_grib.F90
@@ -96,6 +96,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_TG_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_TG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CASE('NCEP ')
+ CALL READ_GRIB_TG_NCEP(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
@@ -108,6 +110,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_WG_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_WG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CASE('NCEP ')
+ CALL READ_GRIB_WG_NCEP(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
@@ -123,6 +127,8 @@ SELECT CASE(HSURF)
CALL READ_GRIB_WGI_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
CALL READ_GRIB_WGI_HIRLAM(HFILE,KLUOUT,ZFIELD,ZD)
+ CASE('NCEP ')
+ CALL READ_GRIB_WGI_NCEP(HFILE,KLUOUT,ZFIELD,ZD)
END SELECT
CALL SOIL_PROFILE_GRIB
!
diff --git a/src/SURFEX/prep_teb_grib.F90 b/src/SURFEX/prep_teb_grib.F90
index d2fb790e664191ad2a710ad8ba5d3b84d3bec8ad..2892d0aa3ef73403028d237a87380557b339c974 100644
--- a/src/SURFEX/prep_teb_grib.F90
+++ b/src/SURFEX/prep_teb_grib.F90
@@ -88,7 +88,7 @@ SELECT CASE(HSURF)
!
CASE('ZS ')
SELECT CASE (CINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_ZS_LAND(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD1D)
ALLOCATE(PFIELD(SIZE(ZFIELD1D),1))
PFIELD(:,1) = ZFIELD1D(:)
@@ -106,7 +106,9 @@ SELECT CASE(HSURF)
CASE('ARPEGE','ALADIN','MOCAGE')
CALL READ_GRIB_TG_METEO_FRANCE(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
CASE('HIRLAM')
- CALL READ_GRIB_TG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CALL READ_GRIB_TG_HIRLAM(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
+ CASE('NCEP ')
+ CALL READ_GRIB_TG_NCEP(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD,ZD)
END SELECT
!* if deep road temperature is prescribed
IF (XTI_ROAD/=XUNDEF) THEN
@@ -120,7 +122,7 @@ SELECT CASE(HSURF)
CASE('T_FLOOR')
!* reading of the profile and its depth definition
SELECT CASE(CINMODEL)
- CASE('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_TF_TEB(HFILE,KLUOUT,CINMODEL,ZTI_BLD,ZMASK,ZFIELD,ZD)
END SELECT
!* if deep road temperature is prescribed
@@ -138,7 +140,7 @@ SELECT CASE(HSURF)
CASE('T_WIN1')
SELECT CASE (CINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_TS(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD1D)
ALLOCATE(PFIELD(NNI,1))
PFIELD(:,1) = ZFIELD1D(:)
@@ -164,7 +166,7 @@ SELECT CASE(HSURF)
!
CASE('T_CAN ')
SELECT CASE (CINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_T2_LAND(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD1D)
ALLOCATE(PFIELD(SIZE(ZFIELD1D),1))
PFIELD(:,1) = ZFIELD1D(:)
@@ -176,7 +178,7 @@ SELECT CASE(HSURF)
!
CASE('Q_CAN ')
SELECT CASE (CINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
ALLOCATE(PFIELD(NNI,1))
PFIELD(:,1) = 0.01
END SELECT
diff --git a/src/SURFEX/prep_watflux_grib.F90 b/src/SURFEX/prep_watflux_grib.F90
index 4971e2ed6cc9c215376457ca37ffcfee06d15506..40899d1d51a00e40490fcfc5799fa75760d0f65f 100644
--- a/src/SURFEX/prep_watflux_grib.F90
+++ b/src/SURFEX/prep_watflux_grib.F90
@@ -78,7 +78,7 @@ SELECT CASE(HSURF)
!
CASE('ZS ')
SELECT CASE (CINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_ZS_LAND(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD)
ALLOCATE(PFIELD(SIZE(ZFIELD),1))
PFIELD(:,1) = ZFIELD(:)
@@ -90,7 +90,7 @@ SELECT CASE(HSURF)
!
CASE('TSWATER')
SELECT CASE (CINMODEL)
- CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM')
+ CASE ('ECMWF ','ARPEGE','ALADIN','MOCAGE','HIRLAM','NCEP ')
CALL READ_GRIB_TSWATER(HFILE,KLUOUT,CINMODEL,ZMASK,ZFIELD)
ALLOCATE(PFIELD(SIZE(ZFIELD),1))
PFIELD(:,1) = ZFIELD(:)