Newer
Older
!MNH_LIC Copyright 1998-2018 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.
!-----------------------------------------------------------------
! #################################
MODULE MODI_READ_ALL_DATA_GRIB_CASE
! #################################
INTERFACE
WAUTELET Philippe
committed
SUBROUTINE READ_ALL_DATA_GRIB_CASE(HFILE,TPPRE_REAL1,HGRIB,TPPGDFILE, &
WAUTELET Philippe
committed
USE MODD_IO_ll, ONLY: TFILEDATA
!
CHARACTER(LEN=4), INTENT(IN) :: HFILE ! which file ('ATM0','ATM1' or 'CHEM')
TYPE(TFILEDATA),POINTER,INTENT(INOUT) :: TPPRE_REAL1 ! PRE_REAL1 file
WAUTELET Philippe
committed
CHARACTER(LEN=28), INTENT(IN) :: HGRIB ! name of the GRIB file
TYPE(TFILEDATA), INTENT(IN) :: TPPGDFILE ! physiographic data file
INTEGER, INTENT(IN) :: KVERB ! verbosity level
LOGICAL, INTENT(IN) :: ODUMMY_REAL ! flag to interpolate dummy fields
REAL, INTENT(INOUT) :: PTIME_HORI ! time spent in hor. interpolations
!
END SUBROUTINE READ_ALL_DATA_GRIB_CASE
!
END INTERFACE
END MODULE MODI_READ_ALL_DATA_GRIB_CASE
! ##########################################################################
WAUTELET Philippe
committed
SUBROUTINE READ_ALL_DATA_GRIB_CASE(HFILE,TPPRE_REAL1,HGRIB,TPPGDFILE, &
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
PTIME_HORI,KVERB,ODUMMY_REAL )
! ##########################################################################
!
!!**** *READ_ALL_DATA_GRIB_CASE* - reads data for the initialization of real cases.
!!
!! PURPOSE
!! -------
! This routine reads the two input files :
! The PGD which is closed after reading
! The GRIB file
! Projection is read in READ_LFIFM_PGD (MODD_GRID).
! Grid and definition of large domain are read in PGD file and Grib files.
! The PGD files are also read in READ_LFIFM_PGD.
! The PGD file is closed.
! The MESO-NH domain is defined from PRE_REAL1.nam inputs in SET_SUBDOMAIN_CEP.
! Vertical grid is defined in READ_VER_GRID.
! PGD fields are stored on MESO-NH domain (in TRUNC_PGD).
!!
!!** METHOD
!! ------
!! 0. Declarations
!! 1. Declaration of arguments
!! 2. Declaration of local variables
!! 1. Read PGD file
!! 1. Domain restriction
!! 2. Coordinate conversion to lat,lon system
!! 2. Read Grib fields
!! 3. Vertical grid
!! 4. Free all temporary allocations
!!
!! EXTERNAL
!! --------
!! subroutine READ_LFIFM_PGD : to read PGD file
!! subroutine SET_SUBDOMAIN : to define the horizontal MESO-NH domain.
!! subroutine READ_VER_GRID : to read the vertical grid in namelist file.
!! subroutine HORIBL : horizontal bilinear interpolation
!! subroutine XYTOLATLON : projection from conformal to lat,lon
!!
!! Module MODI_SET_SUBDOMAIN : interface for subroutine SET_SUBDOMAIN
!! Module MODI_READ_VER_GRID : interface for subroutine READ_VER_GRID
!! Module MODI_HORIBL : interface for subroutine HORIBL
!! Module MODI_XYTOLATLON : interface for subroutine XYTOLATLON
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!!
!! Module MODD_CONF : contains configuration variables for all models.
!! NVERB : verbosity level for output-listing
!! Module MODD_LUNIT : contains logical unit names for all models
!! CLUOUT0 : name of output-listing
!! Module MODD_PGDDIM : contains dimension of PGD fields
!! NPGDIMAX: dimension along x (no external point)
!! NPGDJMAX: dimension along y (no external point)
!! Module MODD_PARAMETERS
!! JPHEXT
!!
!! REFERENCE
!! ---------
!!
!! Book 1 : Informations on ISBA model (soil moisture)
!! "Encoding and decoding Grib data", John D.Chambers, ECMWF, October 95
!! "A guide to Grib", John D.Stackpole, National weather service, March 94
!!
!! AUTHOR
!! ------
!!
!! J. Pettre and V. Bousquet
!!
!! MODIFICATIONS
!! -------------
!! Original 20/11/98
!! 15/03/99 (V. Masson) phasing with new PGD fields
!! 21/04/99 (V. Masson) bug in mask definitions for max Y index
!! 22/04/99 (V. Masson) optimizer bug in u,v loop
!! --> splitting of the loop
!! and splitting of the routine in more
!! contains
!! 28/05/99 (V. Bousquet) bug in wind interpolated variable for
!! Arpege
!! 31/05/99 (V. Masson) set pressure points (given on a regular grid at ECMWF)
!! on orography points (assuming the last are included in the former)
!! pressure computation from parameters A and B
!! (instead of interpolation from grib grid)
!! 20/07/00 (V. Masson) increase the threshold for land_sea index
!! 22/11/00 (P. Tulet) add INTERPOL_SV to initialize SV fields
!! (I. Mallet) from MOCAGE model (IMODE=3)
!! 01/02/01 (D. Gazen) add INI_NSV
!! 18/05/01 (P. Jabouille) problem with 129 grib code
!! 05/12/01 (I. Mallet) add Aladin reunion model
!! 02/10/02 (I. Mallet) 2 orography fields for CEP (SFC, ML=1)
!! 01/12/03 (D. Gazen) change Chemical scheme interface
!! 01/2004 (V. Masson) removes surface (externalization)
!! 01/06/02 (O.Nuissier) filtering of tropical cyclone
!! 01/05/04 (P. Tulet) add INTERPOL_SV to initialize SV dust
!! and aerosol fields
!! 08/06/2010 (G. Tanguy) replace GRIBEX by GRIB_API : change
!! of all the subroutine
!! 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
RODIER Quentin
committed
!! 01/2019 (G.Delautier via Q.Rodier) for GRIB2 ARPEGE and AROME from EPYGRAM
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
!------------
WAUTELET Philippe
committed
!
USE MODE_DATETIME
WAUTELET Philippe
committed
USE MODE_FM, ONLY: IO_FILE_CLOSE_ll
USE MODE_IO_ll, ONLY: UPCASE
WAUTELET Philippe
committed
USE MODE_THERMO
!
USE MODI_READ_HGRID_n
USE MODI_READ_VER_GRID
USE MODI_XYTOLATLON
USE MODI_HORIBL
USE MODI_INI_NSV
USE MODI_CH_INIT_SCHEME_n
USE MODI_REMOVAL_VORTEX
USE MODI_CH_INIT_CCS
USE MODI_CH_AER_INIT_SOA
USE MODI_INI_CTURB
WAUTELET Philippe
committed
USE MODI_CH_OPEN_INPUT
WAUTELET Philippe
committed
USE MODD_IO_ll, ONLY: TFILEDATA
USE MODD_CONF
USE MODD_CONF_n
USE MODD_CST
USE MODD_LUNIT
USE MODD_PARAMETERS
USE MODD_GRID
USE MODD_GRID_n
USE MODD_DIM_n
USE MODD_PARAM_n, ONLY : CTURB
USE MODD_TIME
USE MODD_TIME_n
USE MODD_CH_MNHC_n, ONLY : LUSECHEM,LUSECHAQ,LUSECHIC,LCH_PH
USE MODD_CH_M9_n, ONLY : NEQ , CNAMES
USE MODD_CH_AEROSOL, ONLY: CORGANIC, NCARB, NSOA, NSP, LORILAM,&
JPMODE, LVARSIGI, LVARSIGJ
USE MODD_NSV , ONLY : NSV
USE MODD_HURR_CONF, ONLY : LFILTERING,CFILTERING
USE MODD_PREP_REAL
USE MODE_MODELN_HANDLER
!JUAN REALZ
USE MODE_MPPDB
!JUAN REALZ
WAUTELET Philippe
committed
!
USE GRIB_API
!
IMPLICIT NONE
!
!* 0.1. Declaration of arguments
! ------------------------
!
WAUTELET Philippe
committed
CHARACTER(LEN=4), INTENT(IN) :: HFILE ! which file ('ATM0','ATM1' or 'CHEM')
TYPE(TFILEDATA),POINTER,INTENT(INOUT) :: TPPRE_REAL1! PRE_REAL1 file
WAUTELET Philippe
committed
CHARACTER(LEN=28), INTENT(IN) :: HGRIB ! name of the GRIB file
TYPE(TFILEDATA), INTENT(IN) :: TPPGDFILE ! physiographic data file
INTEGER, INTENT(IN) :: KVERB ! verbosity level
LOGICAL, INTENT(IN) :: ODUMMY_REAL ! flag to interpolate dummy fields
REAL, INTENT(INOUT) :: PTIME_HORI ! time spent in hor. interpolations
!
!* 0.2 Declaration of local variables
! ------------------------------
! General purpose variables
INTEGER :: ILUOUT0 ! Unit used for output msg.
INTEGER :: IRESP ! Return code of FM-routines
INTEGER :: IRET ! Return code from subroutines
WAUTELET Philippe
committed
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 ! |
INTEGER :: JI,JJ ! Dummy counters
INTEGER :: JLOOP1,JLOOP2 ! |
INTEGER :: JLOOP3,JLOOP4 ! |
INTEGER :: JLOOP ! |
! Variables used by the PGD reader
CHARACTER(LEN=28) :: YPGD_NAME ! not used - dummy argument
CHARACTER(LEN=28) :: YPGD_DAD_NAME ! not used - dummy argument
CHARACTER(LEN=2) :: YPGD_TYPE ! not used - dummy argument
! PGD Grib definition variables
INTEGER :: INO ! Number of points of the grid
INTEGER :: IIU ! Number of points along X
INTEGER :: IJU ! Number of points along Y
REAL, DIMENSION(:), ALLOCATABLE :: ZXOUT ! mapping PGD -> Grib (lon.)
REAL, DIMENSION(:), ALLOCATABLE :: ZYOUT ! mapping PGD -> Grib (lat.)
REAL, DIMENSION(:), ALLOCATABLE :: ZLONOUT ! mapping PGD -> Grib (lon.)
REAL, DIMENSION(:), ALLOCATABLE :: ZLATOUT ! mapping PGD -> Grib (lat.)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZXM ! X of PGD mass points
REAL, DIMENSION(:,:), ALLOCATABLE :: ZYM ! Y of PGD mass points
REAL, DIMENSION(:,:), ALLOCATABLE :: ZLATM ! Lat of PGD mass points
REAL, DIMENSION(:,:), ALLOCATABLE :: ZLONM ! Lon of PGD mass points
! Variable involved in the task of reading the grib file
INTEGER(KIND=kindOfInt) :: IUNIT ! unit of the grib file
CHARACTER(LEN=50) :: HGRID ! type of grid
INTEGER :: IPAR ! Parameter identifier
INTEGER :: ITYP ! type of level (Grib code table 3)
INTEGER :: ILEV1 ! level definition
INTEGER :: ILEV2 ! level definition
INTEGER :: IMODEL ! Type of Grib file :
! 0 -> ECMWF
! 1 -> METEO FRANCE - ALADIN/AROME
! 2 -> METEO FRANCE - ALADIN-REUNION
! 3 -> METEO FRANCE - ARPEGE
! 4 -> METEO FRANCE - ARPEGE
! 5 -> METEO FRANCE - MOCAGE
INTEGER :: ICENTER ! number of center
INTEGER :: ISIZE ! size of grib message
INTEGER(KIND=kindOfInt) :: ICOUNT ! number of messages in the file
INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE :: IGRIB ! number of the grib in memory
INTEGER :: INUM ,INUM_ZS ! number of a grib message
REAL,DIMENSION(:),ALLOCATABLE :: ZPARAM ! parameter of girb grid
INTEGER,DIMENSION(:),ALLOCATABLE :: IINLO ! longitude of grib grid
INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE :: IINLO_GRIB ! longitude of grib grid
REAL,DIMENSION(:),ALLOCATABLE :: ZPARAM_ZS ! parameter of girb grid for ZS
INTEGER,DIMENSION(:),ALLOCATABLE :: IINLO_ZS ! longitude of grib grid for ZS
REAL,DIMENSION(:),ALLOCATABLE :: ZVALUE ! Intermediate array
REAL,DIMENSION(:),ALLOCATABLE :: ZOUT ! Intermediate arrays
! Grib grid definition variables
INTEGER :: INI ! Number of points
INTEGER :: INLEVEL ! Number of levels
INTEGER :: ISTARTLEVEL ! First level (0 or 1)
TYPE(DATE_TIME) :: TPTCUR ! Date & time of the grib file data
INTEGER :: ITWOZS
! surface pressure
REAL, DIMENSION(:), ALLOCATABLE :: ZPS_G ! Grib data : Ps
REAL, DIMENSION(:), ALLOCATABLE :: ZLNPS_G ! Grib data : ln(Ps)
REAL, DIMENSION(:), ALLOCATABLE :: ZWORK_LNPS ! Grib data on zs grid: ln(Ps)
INTEGER :: INJ,INJ_ZS
! orography
CHARACTER(LEN=50) :: HGRID_ZS ! type of grid
!
! Reading and projection of the wind vectors u, v
REAL :: ZALPHA ! Angle of rotation
REAL, DIMENSION(:), ALLOCATABLE :: ZTU_LS ! Intermediate array for U
REAL, DIMENSION(:), ALLOCATABLE :: ZTV_LS ! | V
REAL :: ZLATPOLE ! Arpege stretching pole latitude
REAL :: ZLONPOLE ! Arpege stretching pole longitude
REAL :: ZLAT,ZLON ! Lat,lon of current point
REAL :: ZCOS,ZSIN ! cos,sin of rotation matrix
REAL, DIMENSION(:), ALLOCATABLE :: ZTU0_LS ! Arpege temp array for U
REAL, DIMENSION(:), ALLOCATABLE :: ZTV0_LS ! | V
!
! variables for hurricane filtering
REAL, DIMENSION(:,:), ALLOCATABLE :: ZTVF_LS,ZMSLP_LS
REAL :: ZGAMREF ! Standard atmosphere lapse rate (K/m)
! date
INTEGER :: ITIME
INTEGER :: IDATE
INTEGER :: ITIMESTEP
CHARACTER(LEN=10) :: CSTEPUNIT
!chemistery field
WAUTELET Philippe
committed
CHARACTER(LEN=16) :: YPRE_MOC="PRE_MOC1.nam"
INTEGER, DIMENSION(:), ALLOCATABLE :: INUMGRIB, INUMLEV ! grib
INTEGER, DIMENSION(:), ALLOCATABLE :: INUMLEV1, INUMLEV2 !numbers
INTEGER :: IMOC
INTEGER :: IVAR
INTEGER :: ICHANNEL
INTEGER :: INDX
INTEGER :: INACT
CHARACTER(LEN=40) :: YINPLINE ! input line
WAUTELET Philippe
committed
CHARACTER(LEN=16) :: YFIELD
CHARACTER, PARAMETER :: YPTAB = CHAR(9) ! TAB character is ASCII : 9
CHARACTER, PARAMETER :: YPCOM = CHAR(44)! COMma character is ASCII : 44
CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE :: YMNHNAME ! species names
INTEGER :: JN, JNREAL ! loop control variables
CHARACTER(LEN=40) :: YFORMAT
! temperature and humidity
INTEGER :: IT,IQ
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 :: ZT_G ! Temperature
REAL, DIMENSION(:,:), ALLOCATABLE :: ZQ_G ! Specific humidity
REAL, DIMENSION(:), ALLOCATABLE :: ZH_G ! Relative humidity
REAL, DIMENSION(:), ALLOCATABLE :: ZTHV_G ! Theta V
REAL, DIMENSION(:), ALLOCATABLE :: ZRV_G ! Vapor mixing ratio
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZPF_LS ! Pressure (flux point)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZPM_LS ! Pressure (mass point)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEXNF_LS ! Exner fct. (flux point)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEXNM_LS ! Exner fct. (mass point)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZH_LS ! Relative humidity
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRV_LS ! Vapor mixing ratio
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTHV_LS ! Theta V
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTEV_LS ! T V
REAL, DIMENSION(:), ALLOCATABLE :: ZPV ! vertical level in grib file
INTEGER :: IPVPRESENT ,IPV
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZR_DUM
INTEGER :: IMI
WAUTELET Philippe
committed
TYPE(TFILEDATA),POINTER :: TZFILE
INTEGER, DIMENSION(JP_GFS) :: IP_GFS ! list of pressure levels for GFS model
RODIER Quentin
committed
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/)!
WAUTELET Philippe
committed
TZFILE => NULL()
!
IMI = GET_CURRENT_MODEL_INDEX()
!
!* 1. READ PGD FILE
! -------------
!
WAUTELET Philippe
committed
CALL READ_HGRID_n(TPPGDFILE,YPGD_NAME,YPGD_DAD_NAME,YPGD_TYPE)
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
!
! 1.1 Domain restriction
!
!JUAN REALZ
CALL GET_DIM_EXT_ll('B',IIU,IJU)
!!$IIU=NIMAX + 2*JPHEXT
!!$IJU=NJMAX + 2*JPHEXT
!JUAN REALZ
INO = IIU * IJU
!
!
! 1.2 Coordinate conversion to lat,lon system
!
ALLOCATE (ZXM(IIU,IJU))
ALLOCATE (ZYM(IIU,IJU))
ALLOCATE (ZLONM(IIU,IJU))
ALLOCATE (ZLATM(IIU,IJU))
ZXM(1:IIU-1,1) = (XXHAT(1:IIU-1) + XXHAT(2:IIU) ) / 2.
ZXM(IIU,1) = XXHAT(IIU) - XXHAT(IIU-1) + ZXM(IIU-1,1)
ZXM(:,2:IJU) = SPREAD(ZXM(:,1),2,IJU-1)
ZYM(1,1:IJU-1) = (XYHAT(1:IJU-1) + XYHAT(2:IJU)) / 2.
ZYM(1,IJU) = XYHAT(IJU) - XYHAT(IJU-1) + ZYM(1,IJU-1)
ZYM(2:IIU,:) = SPREAD(ZYM(1,:),1,IIU-1)
CALL SM_XYTOLATLON_A (XLAT0,XLON0,XRPK,XLATORI,XLONORI,ZXM,ZYM,ZLATM,ZLONM, &
IIU,IJU)
ALLOCATE (ZLONOUT(INO))
ALLOCATE (ZLATOUT(INO))
JLOOP1 = 0
DO JJ = 1, IJU
ZLONOUT(JLOOP1+1:JLOOP1+IIU) = ZLONM(1:IIU,JJ)
ZLATOUT(JLOOP1+1:JLOOP1+IIU) = ZLATM(1:IIU,JJ)
JLOOP1 = JLOOP1 + IIU
ENDDO
DEALLOCATE (ZLATM)
DEALLOCATE (ZLONM)
DEALLOCATE (ZYM)
DEALLOCATE (ZXM)
!
ALLOCATE (ZXOUT(INO))
ALLOCATE (ZYOUT(INO))
!
!---------------------------------------------------------------------------------------
!
!* 2. READ GRIB FIELDS
! ----------------
!
IF (HFILE(1:3)=='ATM' .OR. HFILE=='CHEM') THEN
WRITE (ILUOUT0,'(A,A4)') ' -- Grib reader started for ',HFILE
ELSE
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','bad input argument')
END IF
!
!* 2.1 Charge in memory the grib messages
!
! open grib file
CALL GRIB_OPEN_FILE(IUNIT,HGRIB,'R',IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error opening the grib file ',TRIM(HGRIB),', error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
! count the messages in the file
CALL GRIB_COUNT_IN_FILE(IUNIT,ICOUNT,IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error in reading the grib file ',TRIM(HGRIB),', error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
ALLOCATE(IGRIB(ICOUNT))
! initialize the tabular with a negativ number
! ( all the IGRIB will be different )
IGRIB(:)=-12
!charge all the message in memory
DO JLOOP=1,ICOUNT
CALL GRIB_NEW_FROM_FILE(IUNIT,IGRIB(JLOOP),IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error in reading the grib file - ILOOP=',JLOOP,' - error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
END DO
! close the grib file
CALL GRIB_CLOSE_FILE(IUNIT)
!
!
!---------------------------------------------------------------------------------------
!* 2.2 Research center with the first message
!---------------------------------------------------------------------------------------
!
CALL GRIB_GET(IGRIB(1),'centre',ICENTER,IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error in reading center - error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
CALL GRIB_GET(IGRIB(1),'typeOfGrid',HGRID,IRET_GRIB)
IF (IRET_GRIB /= 0) THEN
!callabortstop
WRITE(YMSG,*) 'Error in reading type of grid - error code ', IRET_GRIB
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
END IF
!
IMODEL = -1
SELECT CASE (ICENTER)
CASE (98)
WRITE (ILUOUT0,'(A)') &
' | Grib file from European Center for Medium-range Weather Forecast'
IMODEL = 0
ALLOCATE(ZPARAM(6))
CASE (85)
SELECT CASE (HGRID)
CASE('lambert')
RODIER Quentin
committed
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'
IMODEL = 2
ALLOCATE(ZPARAM(9))
CASE('unknown_PLPresent','reduced_stretched_rotated_gg')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arpege model'
ALLOCATE(ZPARAM(10))
RODIER Quentin
committed
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'
WRITE (ILUOUT0,'(A)') 'but same grid as ECMWF model (unstretched)'
IMODEL = 4
ALLOCATE(ZPARAM(10))
CASE('regular_ll')
WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Mocage model'
IMODEL = 5
ALLOCATE(ZPARAM(6))
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
CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE','unsupported Grib file format')
END IF
!
!---------------------------------------------------------------------------------------
!* 2.3 Read and interpol orography
!---------------------------------------------------------------------------------------
!
WRITE (ILUOUT0,'(A)') ' | Searching orography'
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM_ZS,KPARAM=129)
RODIER Quentin
committed
WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
END IF
RODIER Quentin
committed
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
RODIER Quentin
committed
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
RODIER Quentin
committed
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)
WRITE (ILUOUT0,*) 'orog ',IGRIB(INUM_ZS)
END SELECT
ZPARAM(:)=-999.
CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ,IRET_GRIB)
ALLOCATE(IINLO(INJ))
CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM_ZS),IIU,IJU,ZLONOUT,ZLATOUT,&
ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
ALLOCATE(ZPARAM_ZS(SIZE(ZPARAM)))
ZPARAM_ZS=ZPARAM
ALLOCATE(IINLO_ZS(SIZE(IINLO)))
IINLO_ZS=IINLO
CALL GRIB_GET_SIZE(IGRIB(INUM_ZS),'values',ISIZE)
ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM_ZS),'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.)
DEALLOCATE(IINLO)
DEALLOCATE(ZVALUE)
!
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
ALLOCATE (XZS_LS(IIU,IJU))
ALLOCATE (XZSMT_LS(IIU,IJU))
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XZS_LS)
XZSMT_LS = XZS_LS ! no smooth orography in this case
ELSE IF (HFILE=='CHEM') THEN
ALLOCATE (XZS_SV_LS(IIU,IJU))
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XZS_SV_LS)
END IF
DEALLOCATE (ZOUT)
!
!---------------------------------------------------------------------------------------
!* 2.4 Interpolation surface pressure
!---------------------------------------------------------------------------------------
!
!* 2.4.1 Read pressure
!
WRITE (ILUOUT0,'(A)') ' | Searching pressure'
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=152)
CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=1)
CASE(6,7) ! NEW AROME,ARPEGE
CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=3,KNUMBER=25)
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=134)
END SELECT
IF(INUM < 0) THEN
WRITE (ILUOUT0,'(A)')'Surface pressure is missing - abort'
CALL ABORT
STOP
ENDIF
! recuperation du tableau de valeurs
CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
ALLOCATE(ZVALUE(ISIZE))
CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
! determination des tableaux ZPS_G et ZLNPS_G
SELECT CASE (IMODEL)
RODIER Quentin
committed
CASE(0,6,7) ! ECMWF
ALLOCATE (ZPS_G (ISIZE))
ALLOCATE (ZLNPS_G(ISIZE))
ZLNPS_G(:) = ZVALUE(1:ISIZE)
ZPS_G (:) = EXP(ZVALUE(1:ISIZE))
CASE(1,2,3,4,5,10) ! arpege mocage aladin aladin-reunion NCEP
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
ALLOCATE (ZPS_G (ISIZE))
ALLOCATE (ZLNPS_G(ISIZE))
ZPS_G (:) = ZVALUE(1:ISIZE)
ZLNPS_G(:) = LOG(ZVALUE(1:ISIZE))
END SELECT
DEALLOCATE (ZVALUE)
!
!* 2.4.2 Removes pressure points not on orography points
! (if pressure is on a regular grid)
CALL GRIB_GET(IGRIB(INUM),'typeOfGrid',HGRID)
CALL GRIB_GET(IGRIB(INUM_ZS),'typeOfGrid',HGRID_ZS)
CALL GRIB_GET(IGRIB(INUM),'Nj',INJ)
CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ_ZS)
!
IF ( HGRID(1:7)=='regular' .AND. HGRID_ZS(1:7)=='reduced' .AND.&
INJ == INJ_ZS) THEN
WRITE (ILUOUT0,'(A)')'HGRID(1:7)==regular .AND. HGRID_ZS(1:7)==reduced .AND. INJ == INJ_ZS - abort'
CALL ABORT
STOP
ELSE
ALLOCATE(ZWORK_LNPS(SIZE(ZLNPS_G)))
ZWORK_LNPS(:) = ZLNPS_G(:)
ENDIF
!
IF (HFILE(1:3)=='ATM') THEN
ALLOCATE (XPS_LS(IIU,IJU))
ELSE IF (HFILE=='CHEM') THEN
ALLOCATE (XPS_SV_LS(IIU,IJU))
END IF
!
ALLOCATE(IINLO(INJ))
CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
ALLOCATE(ZOUT(INO))
CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI,&
ZWORK_LNPS,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE. )
DEALLOCATE(ZWORK_LNPS)
DEALLOCATE(IINLO)
!
!* 2.4.3 conversion to surface pressure
!
ZOUT=EXP(ZOUT)
IF (HFILE(1:3)=='ATM') THEN
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XPS_LS(:,:))
ELSE IF (HFILE=='CHEM') THEN
CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XPS_SV_LS(:,:))
END IF
!JUAN REALZ
CALL MPPDB_CHECK2D(XZS_LS,"XZS_LS",PRECISION)
!JUAN REALZ
DEALLOCATE (ZOUT)
DEALLOCATE (ZLNPS_G)
!
!---------------------------------------------------------------------------------------
!* 2.5 Interpolation temperature and humidity
!---------------------------------------------------------------------------------------
!
!* 2.5.1 Read T and Q on each level
!
WRITE (ILUOUT0,'(A)') ' | Reading T and Q fields'
!
RODIER Quentin
committed
IF (IMODEL/=10) THEN
SELECT CASE (IMODEL)
CASE(0) ! ECMWF
ISTARTLEVEL=1
IT=130
IQ=133
RODIER Quentin
committed
CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
RODIER Quentin
committed
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
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)
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))
!
IF (IMODEL/=10) THEN ! others than NCEP
DO JLOOP1=1, INLEVEL
ILEV1 = JLOOP1-1+ISTARTLEVEL
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ILEV1)
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))
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ILEV1)
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)
RODIER Quentin
committed
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(:,JLOOP1),IRET_GRIB)
WRITE (ILUOUT0,*) 'Q ',ILEV1,IRET_GRIB
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ILEV1)
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 (IMODEL/=10) THEN ! others than NCEP
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)
RODIER Quentin
committed
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)
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
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 (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)
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)
!
ALLOCATE (ZEXNF_G(INI,INLEVEL))
ZEXNF_G(:,:) = (ZPF_G(:,:)/XP00)**(XRD/XCPD)
!
ALLOCATE (ZEXNM_G(INI,INLEVEL))
ZEXNM_G(:,1:INLEVEL-1) = (ZEXNF_G(:,1:INLEVEL-1)-ZEXNF_G(:,2:INLEVEL)) / &
(LOG(ZEXNF_G(:,1:INLEVEL-1))-LOG(ZEXNF_G(:,2:INLEVEL)))
ZEXNM_G(:,INLEVEL) = (ZPF_G(:,INLEVEL)/2./XP00)**(XRD/XCPD)
DEALLOCATE (ZEXNF_G)
DEALLOCATE (ZPF_G)
!
ALLOCATE (ZPM_G(INI,INLEVEL))
ZPM_G(:,:) = XP00*(ZEXNM_G(:,:))**(XCPD/XRD)
!
!* 2.5.4 Compute the relative humidity and the interpolate Thetav, P, Q and H
!
IF (IMODEL==1) THEN
! search cloud_water in Arome case (forecast)
ISTARTLEVEL = 1
IPAR=246
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
IF (INUM < 0) THEN
ISTARTLEVEL = 0
RODIER Quentin
committed
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)'
LCPL_AROME=.TRUE.
NRR=6
END IF
! search also turbulent kinetic energy
ISTARTLEVEL = 1
IPAR=251
RODIER Quentin
committed
CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IPAR,KLEV1=ISTARTLEVEL)
IF (INUM < 0) THEN
ISTARTLEVEL = 0
RODIER Quentin
committed
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)
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
END IF
IF (INUM > 0) CTURB='TKEL'
END IF
!
!
WRITE (ILUOUT0,'(A)') ' | Computing relative humidity on each level'
ALLOCATE (ZH_G(INI))
ALLOCATE (ZH_LS(IIU,IJU,INLEVEL))
IF (HFILE(1:3)=='ATM') THEN
ALLOCATE (XT_LS(IIU,IJU,INLEVEL))
ALLOCATE (XQ_LS(IIU,IJU,INLEVEL,NRR)) ; XQ_LS=0.
ELSE IF (HFILE=='CHEM') THEN
ALLOCATE (XT_SV_LS(IIU,IJU,INLEVEL))
ALLOCATE (XQ_SV_LS(IIU,IJU,INLEVEL,1))
END IF
IF (CTURB=='TKEL') THEN
IF (ALLOCATED(XTKE_LS)) DEALLOCATE(XTKE_LS)
ALLOCATE(XTKE_LS(IIU,IJU,INLEVEL)) ; XTKE_LS=0.
CALL INI_CTURB
ELSE
IF (ALLOCATED(XTKE_LS)) DEALLOCATE(XTKE_LS)
ALLOCATE(XTKE_LS(0,0,0))
END IF
ALLOCATE (ZTHV_LS(IIU,IJU,INLEVEL))
ALLOCATE (ZTHV_G(INI))
ALLOCATE (ZRV_G(INI))
ALLOCATE (ZOUT(INO))
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,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
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 (IMODEL/=10) THEN ! others than NCEP
IF (HFILE(1:3)=='ATM') THEN
ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_LS,1,IIU),2,IJU)*XP00_LS + &