b7c9ac81447cee60a2365fba62a7a1f81b9ce336
[MNH-git_open_source-lfs.git] / src / MNH / read_all_data_grib_case.f90
1 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
4 !MNH_LIC for details. version 1.
5 !-----------------------------------------------------------------
6 !--------------- special set of characters for RCS information
7 !-----------------------------------------------------------------
8 ! $Source$ $Revision$
9 ! masdev4_7 BUG3 2007/11/20 12:25:23
10 !-----------------------------------------------------------------
11 !     #################################
12       MODULE MODI_READ_ALL_DATA_GRIB_CASE
13 !     #################################
14 INTERFACE
15 SUBROUTINE READ_ALL_DATA_GRIB_CASE(HFILE,HPRE_REAL1,HGRIB,HPGDFILE,      &
16                     PTIME_HORI,KVERB,ODUMMY_REAL                         ) 
17 !
18 CHARACTER(LEN=4),   INTENT(IN) :: HFILE    !which file ('ATM0','ATM1' or 'CHEM')
19 CHARACTER(LEN=28),  INTENT(IN) :: HPRE_REAL1 ! name of the PRE_REAL1 file
20 CHARACTER(LEN=28),  INTENT(IN) :: HGRIB    ! name of the GRIB file
21 CHARACTER(LEN=28),  INTENT(IN) :: HPGDFILE ! name of the physiographic data file
22 INTEGER,           INTENT(IN)  :: KVERB    ! verbosity level
23 LOGICAL,           INTENT(IN)  :: ODUMMY_REAL! flag to interpolate dummy fields
24 REAL,           INTENT(INOUT)  :: PTIME_HORI ! time spent in hor. interpolations
25 END SUBROUTINE READ_ALL_DATA_GRIB_CASE
26 !
27 END INTERFACE
28 END MODULE MODI_READ_ALL_DATA_GRIB_CASE
29 !     ##########################################################################
30       SUBROUTINE READ_ALL_DATA_GRIB_CASE(HFILE,HPRE_REAL1,HGRIB,HPGDFILE,      &
31                        PTIME_HORI,KVERB,ODUMMY_REAL                            )
32 !     ##########################################################################
33 !
34 !!****  *READ_ALL_DATA_GRIB_CASE* - reads data for the initialization of real cases.
35 !!
36 !!    PURPOSE
37 !!    -------
38 !     This routine reads the two input files :
39 !       The PGD which is closed after reading
40 !       The GRIB file
41 !     Projection is read in READ_LFIFM_PGD (MODD_GRID).
42 !     Grid and definition of large domain are read in PGD file and Grib files.
43 !     The PGD files are also read in READ_LFIFM_PGD.
44 !     The PGD file is closed.
45 !     The MESO-NH domain is defined from PRE_REAL1.nam inputs in SET_SUBDOMAIN_CEP.
46 !     Vertical grid is defined in READ_VER_GRID.
47 !     PGD fields are stored on MESO-NH domain (in TRUNC_PGD).
48 !!
49 !!**  METHOD
50 !!    ------
51 !!  0. Declarations
52 !!    1. Declaration of arguments
53 !!    2. Declaration of local variables
54 !!  1. Read PGD file
55 !!    1. Domain restriction
56 !!    2. Coordinate conversion to lat,lon system
57 !!  2. Read Grib fields
58 !!  3. Vertical grid
59 !!  4. Free all temporary allocations
60 !!
61 !!    EXTERNAL
62 !!    --------
63 !!    subroutine READ_LFIFM_PGD    : to read PGD file
64 !!    subroutine SET_SUBDOMAIN     : to define the horizontal MESO-NH domain.
65 !!    subroutine READ_VER_GRID     : to read the vertical grid in namelist file.
66 !!    subroutine HORIBL            : horizontal bilinear interpolation
67 !!    subroutine XYTOLATLON        : projection from conformal to lat,lon
68 !!
69 !!    function   FMLOOK_ll         : to retrieve the logical unit associated with a file
70 !!
71 !!    Module     MODI_SET_SUBDOMAIN     : interface for subroutine SET_SUBDOMAIN
72 !!    Module     MODI_READ_VER_GRID     : interface for subroutine READ_VER_GRID
73 !!    Module     MODI_HORIBL            : interface for subroutine HORIBL
74 !!    Module     MODI_XYTOLATLON        : interface for subroutine XYTOLATLON
75 !!
76 !!    IMPLICIT ARGUMENTS
77 !!    ------------------
78 !!
79 !!      Module MODD_CONF      : contains configuration variables for all models.
80 !!         NVERB : verbosity level for output-listing
81 !!      Module MODD_LUNIT     : contains logical unit names for all models
82 !!         CLUOUT0 : name of output-listing
83 !!      Module MODD_PGDDIM    : contains dimension of PGD fields
84 !!         NPGDIMAX: dimension along x (no external point)
85 !!         NPGDJMAX: dimension along y (no external point)
86 !!      Module MODD_PARAMETERS
87 !!         JPHEXT
88 !!
89 !!    REFERENCE
90 !!    ---------
91 !!
92 !!      Book 1 : Informations on ISBA model (soil moisture)
93 !!      "Encoding and decoding Grib data", John D.Chambers, ECMWF, October 95
94 !!      "A guide to Grib", John D.Stackpole, National weather service, March 94
95 !!
96 !!    AUTHOR
97 !!    ------
98 !!
99 !!      J. Pettre and V. Bousquet
100 !!
101 !!    MODIFICATIONS
102 !!    -------------
103 !!      Original    20/11/98
104 !!                  15/03/99 (V. Masson) phasing with new PGD fields
105 !!                  21/04/99 (V. Masson) bug in mask definitions for max Y index
106 !!                  22/04/99 (V. Masson) optimizer bug in u,v loop
107 !!                                       --> splitting of the loop
108 !!                                       and splitting of the routine in more
109 !!                                       contains
110 !!                  28/05/99 (V. Bousquet) bug in wind interpolated variable for
111 !!                                         Arpege
112 !!                  31/05/99 (V. Masson) set pressure points (given on a regular grid at ECMWF)
113 !!                             on orography points (assuming the last are included in the former)
114 !!                                       pressure computation from parameters A and B
115 !!                                        (instead of interpolation from grib grid)
116 !!                  20/07/00 (V. Masson) increase the threshold for land_sea index
117 !!                  22/11/00 (P. Tulet) add INTERPOL_SV to initialize SV fields
118 !!                           (I. Mallet) from MOCAGE model (IMODE=3)
119 !!                  01/02/01 (D. Gazen) add INI_NSV
120 !!                  18/05/01 (P. Jabouille) problem with 129 grib code
121 !!                  05/12/01 (I. Mallet) add Aladin reunion model
122 !!                  02/10/02 (I. Mallet) 2 orography fields for CEP (SFC, ML=1)
123 !!                  01/12/03 (D. Gazen)   change Chemical scheme interface
124 !!                  01/2004  (V. Masson) removes surface (externalization)
125 !!                  01/06/02 (O.Nuissier) filtering of tropical cyclone
126 !!                  01/05/04 (P. Tulet) add INTERPOL_SV to initialize SV dust
127 !!                                                         and aerosol fields
128 !!                  08/06/2010 (G. Tanguy) replace GRIBEX by GRIB_API : change
129 !!                                         of all the subroutine
130 !-------------------------------------------------------------------------------
131 !
132 !*      0. DECLARATIONS
133 !------------
134 USE MODE_FM
135 USE MODE_IO_ll
136 USE MODE_TIME
137 !
138 USE MODI_ADD_FORECAST_TO_DATE
139 USE MODI_READ_HGRID_n
140 USE MODI_READ_VER_GRID
141 USE MODI_XYTOLATLON
142 USE MODI_HORIBL
143 USE MODI_INI_NSV
144 USE MODI_CH_INIT_SCHEME_n
145 USE MODI_REMOVAL_VORTEX
146 USE MODI_CH_INIT_CCS
147 USE MODI_CH_AER_INIT_SOA
148 USE MODI_INI_CTURB
149 !
150 USE MODD_CONF
151 USE MODD_CONF_n
152 USE MODD_CST
153 USE MODD_LUNIT
154 USE MODD_PARAMETERS
155 USE MODD_GRID
156 USE MODD_GRID_n
157 USE MODD_DIM_n
158 USE MODD_PARAM_n, ONLY : CTURB
159 USE MODD_TIME
160 USE MODD_TIME_n
161 USE MODD_CH_MNHC_n, ONLY : LUSECHEM,LUSECHAQ,LUSECHIC,LCH_PH
162 USE MODD_CH_M9_n,   ONLY : NEQ ,  CNAMES
163 USE MODD_CH_AEROSOL, ONLY: CORGANIC, NCARB, NSOA, NSP, LORILAM,&
164                            JPMODE, LVARSIGI, LVARSIGJ
165 USE MODD_NSV      , ONLY : NSV
166 USE MODD_HURR_CONF, ONLY : LFILTERING,CFILTERING
167 USE MODD_PREP_REAL
168 USE MODE_MODELN_HANDLER
169 !JUAN REALZ
170 USE MODE_MPPDB
171 !JUAN REALZ
172 USE GRIB_API
173 USE MODI_CH_OPEN_INPUT  
174 USE MODE_THERMO
175
176 !
177 IMPLICIT NONE
178 !
179 !* 0.1. Declaration of arguments
180 !       ------------------------
181 !
182 CHARACTER(LEN=4),       INTENT(IN) :: HFILE !which file ('ATM0','ATM1' or 'CHEM')
183 CHARACTER(LEN=28),      INTENT(IN) :: HPRE_REAL1 ! name of the PRE_REAL1 file
184 CHARACTER(LEN=28),      INTENT(IN) :: HGRIB      ! name of the GRIB file
185 CHARACTER(LEN=28),      INTENT(IN) :: HPGDFILE   ! name of the physiographic data file
186 INTEGER,               INTENT(IN)  :: KVERB    ! verbosity level
187 LOGICAL,               INTENT(IN)  :: ODUMMY_REAL! flag to interpolate dummy fields
188 REAL,                INTENT(INOUT) :: PTIME_HORI ! time spent in hor. interpolations
189 !
190 !* 0.2 Declaration of local variables
191 !      ------------------------------
192 ! General purpose variables
193 INTEGER                            :: ILUOUT0       ! Unit used for output msg.
194 INTEGER                            :: IRESP   ! Return code of FM-routines
195 INTEGER                            :: IRET          ! Return code from subroutines
196 INTEGER(KIND=kindOfInt)                  :: IRET_GRIB          ! Return code from subroutines
197 REAL                               :: ZA,ZB,ZC      ! Dummy variables
198 REAL                               :: ZD,ZE,ZF      !  |
199 REAL                               :: ZTEMP         !  |
200 INTEGER                            :: JI,JJ         ! Dummy counters
201 INTEGER                            :: JLOOP1,JLOOP2 !  |
202 INTEGER                            :: JLOOP3,JLOOP4 !  |
203 INTEGER                            :: JLOOP         !  |
204 ! Variables used by the PGD reader
205 CHARACTER(LEN=28)                  :: YPGD_NAME     ! not used - dummy argument
206 CHARACTER(LEN=28)                  :: YPGD_DAD_NAME ! not used - dummy argument
207 CHARACTER(LEN=2)                   :: YPGD_TYPE     ! not used - dummy argument
208 ! PGD Grib definition variables
209 INTEGER                            :: INO           ! Number of points of the grid
210 INTEGER                            :: IIU           ! Number of points along X
211 INTEGER                            :: IJU           ! Number of points along Y
212 REAL, DIMENSION(:), ALLOCATABLE    :: ZXOUT         ! mapping PGD -> Grib (lon.)
213 REAL, DIMENSION(:), ALLOCATABLE    :: ZYOUT         ! mapping PGD -> Grib (lat.)
214 REAL, DIMENSION(:), ALLOCATABLE    :: ZLONOUT       ! mapping PGD -> Grib (lon.)
215 REAL, DIMENSION(:), ALLOCATABLE    :: ZLATOUT       ! mapping PGD -> Grib (lat.)
216 REAL, DIMENSION(:,:), ALLOCATABLE  :: ZXM           ! X of PGD mass points
217 REAL, DIMENSION(:,:), ALLOCATABLE  :: ZYM           ! Y of PGD mass points
218 REAL, DIMENSION(:,:), ALLOCATABLE  :: ZLATM         ! Lat of PGD mass points
219 REAL, DIMENSION(:,:), ALLOCATABLE  :: ZLONM         ! Lon of PGD mass points
220 ! Variable involved in the task of reading the grib file
221 INTEGER(KIND=kindOfInt)                            :: IUNIT         ! unit of the grib file
222 CHARACTER(LEN=20)                  :: HGRID         ! type of grid
223 INTEGER                            :: IPAR          ! Parameter identifier
224 INTEGER                            :: ITYP          ! type of level (Grib code table 3)
225 INTEGER                            :: ILEV1         ! level definition
226 INTEGER                            :: ILEV2         ! level definition
227 INTEGER                            :: IMODEL        ! Type of Grib file :
228                                                     ! 0 -> ECMWF
229                                                     ! 1 -> METEO FRANCE - ALADIN/AROME
230                                                     ! 2 -> METEO FRANCE - ALADIN-REUNION
231                                                     ! 3 -> METEO FRANCE - ARPEGE
232                                                     ! 4 -> METEO FRANCE - ARPEGE
233                                                     ! 5 -> METEO FRANCE - MOCAGE
234 INTEGER                            :: ICENTER       ! number of center
235 INTEGER                            :: ISIZE         ! size of grib message
236 INTEGER(KIND=kindOfInt)                            :: ICOUNT        ! number of messages in the file
237 INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE   :: IGRIB         ! number of the grib in memory
238 INTEGER                            :: INUM ,INUM_ZS ! number of a grib message 
239 REAL,DIMENSION(:),ALLOCATABLE      :: ZPARAM        ! parameter of girb grid
240 INTEGER,DIMENSION(:),ALLOCATABLE   :: IINLO         ! longitude of grib grid
241 INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE   :: IINLO_GRIB         ! longitude of grib grid
242 REAL,DIMENSION(:),ALLOCATABLE      :: ZPARAM_ZS     ! parameter of girb grid for ZS
243 INTEGER,DIMENSION(:),ALLOCATABLE   :: IINLO_ZS      ! longitude of grib grid for ZS
244 REAL,DIMENSION(:),ALLOCATABLE      :: ZVALUE        ! Intermediate array
245 REAL,DIMENSION(:),ALLOCATABLE      :: ZOUT          ! Intermediate arrays
246 ! Grib grid definition variables
247 INTEGER                            :: INI           ! Number of points
248 INTEGER                            :: INLEVEL       ! Number of levels
249 INTEGER                            :: ISTARTLEVEL   ! First level (0 or 1)
250 TYPE(DATE_TIME)                    :: TPTCUR        ! Date & time of the grib file data
251 INTEGER                            :: ITWOZS
252 ! surface pressure
253 REAL, DIMENSION(:), ALLOCATABLE    :: ZPS_G         ! Grib data : Ps
254 REAL, DIMENSION(:), ALLOCATABLE    :: ZLNPS_G  ! Grib data : ln(Ps)
255 REAL, DIMENSION(:), ALLOCATABLE    :: ZWORK_LNPS ! Grib data on zs grid: ln(Ps)
256 INTEGER                            :: INJ,INJ_ZS
257 ! orography
258 CHARACTER(LEN=20)                  :: HGRID_ZS         ! type of grid
259 !
260 ! Reading and projection of the wind vectors u, v
261 REAL                               :: ZALPHA        ! Angle of rotation
262 REAL, DIMENSION(:), ALLOCATABLE    :: ZTU_LS        ! Intermediate array for U
263 REAL, DIMENSION(:), ALLOCATABLE    :: ZTV_LS        !  |                     V
264 REAL                               :: ZLATPOLE      ! Arpege stretching pole latitude
265 REAL                               :: ZLONPOLE      ! Arpege stretching pole longitude
266 REAL                               :: ZLAT,ZLON     ! Lat,lon of current point
267 REAL                               :: ZCOS,ZSIN     ! cos,sin of rotation matrix
268 REAL, DIMENSION(:), ALLOCATABLE    :: ZTU0_LS       ! Arpege temp array for U
269 REAL, DIMENSION(:), ALLOCATABLE    :: ZTV0_LS       !  |                    V
270 !
271 ! variables for hurricane filtering
272 REAL, DIMENSION(:,:), ALLOCATABLE  :: ZTVF_LS,ZMSLP_LS
273 REAL                               :: ZGAMREF       ! Standard atmosphere lapse rate (K/m)
274 ! date
275 INTEGER  :: ITIME
276 INTEGER  :: IDATE
277 INTEGER  :: ITIMESTEP
278 CHARACTER(LEN=10) :: CSTEPUNIT
279 !chemistery field
280 CHARACTER(LEN=12)                  :: YPRE_MOC="PRE_MOC1.nam"
281 INTEGER, DIMENSION(:), ALLOCATABLE :: INUMGRIB, INUMLEV  ! grib
282 INTEGER, DIMENSION(:), ALLOCATABLE :: INUMLEV1, INUMLEV2 !numbers
283 INTEGER                            :: IMOC
284 INTEGER                            :: IVAR
285 INTEGER                            :: ICHANNEL
286 INTEGER                            :: INDX
287 INTEGER                            :: INACT
288 CHARACTER(LEN=40)                  :: YINPLINE        ! input line
289 CHARACTER(LEN=16)                  :: YFIELD
290 CHARACTER, PARAMETER               :: YPTAB = CHAR(9) ! TAB character is ASCII : 9
291 CHARACTER, PARAMETER               :: YPCOM = CHAR(44)! COMma character is ASCII : 44
292 CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE :: YMNHNAME ! species names
293 INTEGER                            :: JN, JNREAL ! loop control variables
294 CHARACTER(LEN=40)                  :: YFORMAT
295 ! temperature and humidity
296 INTEGER                             :: IT,IQ
297 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZPF_G    ! Pressure (flux point)
298 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZPM_G    ! Pressure (mass point)
299 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZEXNF_G  ! Exner fct. (flux point)
300 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZEXNM_G  ! Exner fct. (mass point)
301 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZT_G     ! Temperature
302 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZQ_G     ! Specific humidity
303 REAL, DIMENSION(:), ALLOCATABLE     :: ZH_G     ! Relative humidity
304 REAL, DIMENSION(:), ALLOCATABLE     :: ZTHV_G   ! Theta V
305 REAL, DIMENSION(:), ALLOCATABLE     :: ZRV_G    ! Vapor mixing ratio
306 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZPF_LS   ! Pressure (flux point)
307 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZPM_LS   ! Pressure (mass point)
308 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEXNF_LS ! Exner fct. (flux point)
309 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEXNM_LS ! Exner fct. (mass point)
310 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZH_LS    ! Relative humidity
311 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRV_LS   ! Vapor mixing ratio
312 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTHV_LS  ! Theta V
313 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTEV_LS  ! T V
314 REAL, DIMENSION(:), ALLOCATABLE     :: ZPV      ! vertical level in grib file
315 INTEGER                             :: IPVPRESENT ,IPV
316 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZR_DUM
317 INTEGER                            :: IMI
318 !
319 !---------------------------------------------------------------------------------------
320 !
321 IMI = GET_CURRENT_MODEL_INDEX()
322 !
323 !* 1. READ PGD FILE
324 !     -------------
325 !
326 CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRET)
327 CALL READ_HGRID_n(HPGDFILE,YPGD_NAME,YPGD_DAD_NAME,YPGD_TYPE)
328 !
329 ! 1.1 Domain restriction
330 !
331 !JUAN REALZ
332 CALL GET_DIM_EXT_ll('B',IIU,IJU)
333 !!$IIU=NIMAX + 2*JPHEXT
334 !!$IJU=NJMAX + 2*JPHEXT
335 !JUAN REALZ
336 INO = IIU * IJU
337 !
338 !
339 ! 1.2 Coordinate conversion to lat,lon system
340 !
341 ALLOCATE (ZXM(IIU,IJU))
342 ALLOCATE (ZYM(IIU,IJU))
343 ALLOCATE (ZLONM(IIU,IJU))
344 ALLOCATE (ZLATM(IIU,IJU))
345 ZXM(1:IIU-1,1) = (XXHAT(1:IIU-1) + XXHAT(2:IIU) ) / 2.
346 ZXM(IIU,1)     = XXHAT(IIU) - XXHAT(IIU-1) + ZXM(IIU-1,1)
347 ZXM(:,2:IJU)   = SPREAD(ZXM(:,1),2,IJU-1)
348 ZYM(1,1:IJU-1) = (XYHAT(1:IJU-1) + XYHAT(2:IJU)) / 2.
349 ZYM(1,IJU)     = XYHAT(IJU) - XYHAT(IJU-1) + ZYM(1,IJU-1)
350 ZYM(2:IIU,:)   = SPREAD(ZYM(1,:),1,IIU-1)
351 CALL SM_XYTOLATLON_A (XLAT0,XLON0,XRPK,XLATORI,XLONORI,ZXM,ZYM,ZLATM,ZLONM, &
352                       IIU,IJU)
353 ALLOCATE (ZLONOUT(INO))
354 ALLOCATE (ZLATOUT(INO))
355 JLOOP1 = 0
356 DO JJ = 1, IJU
357   ZLONOUT(JLOOP1+1:JLOOP1+IIU) = ZLONM(1:IIU,JJ)
358   ZLATOUT(JLOOP1+1:JLOOP1+IIU) = ZLATM(1:IIU,JJ)
359   JLOOP1 = JLOOP1 + IIU
360 ENDDO
361 DEALLOCATE (ZLATM)
362 DEALLOCATE (ZLONM)
363 DEALLOCATE (ZYM)
364 DEALLOCATE (ZXM)
365 !
366 ALLOCATE (ZXOUT(INO))
367 ALLOCATE (ZYOUT(INO))
368 !
369 !---------------------------------------------------------------------------------------
370 !
371 !* 2. READ GRIB FIELDS
372 !     ----------------
373 !
374 IF (HFILE(1:3)=='ATM' .OR. HFILE=='CHEM') THEN
375   WRITE (ILUOUT0,'(A,A4)') ' -- Grib reader started for ',HFILE
376 ELSE
377   WRITE (ILUOUT0,'(A)') ' -> Bad input argument in read_all_data_grib_case - abort'
378  !callabortstop
379   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
380   CALL ABORT
381   STOP
382 END IF
383 !
384 !* 2.1 Charge in memory the grib messages
385 !
386 ! open grib file
387 CALL GRIB_OPEN_FILE(IUNIT,HGRIB,'R',IRET_GRIB)
388 IF (IRET_GRIB /= 0) THEN
389   WRITE (ILUOUT0,'(A,A,A,I2)') ' -> Error opening the grib file ',HGRIB,', error code ', IRET_GRIB
390   !callabortstop
391   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
392   CALL ABORT
393   STOP
394 END IF
395 ! count the messages in the file
396 CALL GRIB_COUNT_IN_FILE(IUNIT,ICOUNT,IRET_GRIB)
397 IF (IRET_GRIB /= 0) THEN
398   WRITE (ILUOUT0,'(A,I2)')' -> Error in reading the grib file - error code ', IRET_GRIB
399   !callabortstop
400   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
401   CALL ABORT
402   STOP
403 END IF
404 ALLOCATE(IGRIB(ICOUNT))
405 ! initialize the tabular with a negativ number 
406 ! ( all the IGRIB will be  different )
407 IGRIB(:)=-12
408 !charge all the message in memory
409 DO JLOOP=1,ICOUNT
410 CALL GRIB_NEW_FROM_FILE(IUNIT,IGRIB(JLOOP),IRET_GRIB)
411 IF (IRET_GRIB /= 0) THEN
412   WRITE (ILUOUT0,'(A,I3,A,I2)') ' -> Error in reading the grib file - ILOOP= ',JLOOP,' - error code ', IRET_GRIB
413   !callabortstop
414   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
415   CALL ABORT
416   STOP
417 END IF
418 END DO
419 ! close the grib file
420 CALL GRIB_CLOSE_FILE(IUNIT)
421 !
422 !
423 !---------------------------------------------------------------------------------------
424 !* 2.2 Research center with the first message
425 !---------------------------------------------------------------------------------------
426 !
427 CALL GRIB_GET(IGRIB(1),'centre',ICENTER,IRET_GRIB)
428 IF (IRET_GRIB /= 0) THEN
429   WRITE (ILUOUT0,'(A,I2)')' -> Error in reading center - error code ', IRET_GRIB
430   !callabortstop
431   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
432   CALL ABORT
433   STOP
434 END IF
435 CALL GRIB_GET(IGRIB(1),'typeOfGrid',HGRID,IRET_GRIB)
436 IF (IRET_GRIB /= 0) THEN
437   WRITE (ILUOUT0,'(A,I2)')' -> Error in reading type of grid - error code ', IRET_GRIB
438   !callabortstop
439   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
440   CALL ABORT
441   STOP
442 END IF
443 !
444 IMODEL = -1
445 SELECT CASE (ICENTER)
446   CASE (98)
447     WRITE (ILUOUT0,'(A)') &
448         ' | Grib file from European Center for Medium-range Weather Forecast'
449     IMODEL = 0
450     ALLOCATE(ZPARAM(6))
451   CASE (85)
452     SELECT CASE (HGRID)      
453       CASE('lambert')
454         WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Aladin france model'
455         IMODEL = 1
456         ALLOCATE(ZPARAM(10))
457       CASE('mercator')
458         WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Aladin reunion model'
459       IMODEL = 2
460       ALLOCATE(ZPARAM(9))
461
462     CASE('unknown_PLPresent')
463       WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arpege model'
464       IMODEL = 3
465       ALLOCATE(ZPARAM(10))
466
467     CASE('regular_gg')
468       WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arpege model'
469       WRITE (ILUOUT0,'(A)') 'but same grid as ECMWF model (unstretched)'
470       IMODEL = 4
471       ALLOCATE(ZPARAM(10))
472     CASE('regular_ll')  
473       WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Mocage model'
474       IMODEL = 5
475        ALLOCATE(ZPARAM(6))
476     END SELECT
477 END SELECT
478 IF (IMODEL==-1) THEN
479   WRITE (ILUOUT0,'(A)') ' -> Unsupported Grib file format - abort'
480 !callabortstop
481   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
482   CALL ABORT
483   STOP
484 END IF
485 !
486 !---------------------------------------------------------------------------------------
487 !* 2.3 Read and interpol orography 
488 !---------------------------------------------------------------------------------------
489 !
490 WRITE (ILUOUT0,'(A)') ' | Searching orography'
491 SELECT CASE (IMODEL)
492   CASE(0) ! ECMWF
493      CALL SEARCH_FIELD(129,109,1,0,IGRIB,INUM_ZS)
494      IF(INUM_ZS < 0) THEN
495      CALL SEARCH_FIELD(129,1,0,0,IGRIB,INUM_ZS)
496        IF(INUM_ZS < 0) THEN
497          WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
498        END IF 
499       END IF
500   CASE(3,4,5) ! arpege et mocage
501       CALL SEARCH_FIELD(8,-1,-1,-1,IGRIB,INUM_ZS)
502        IF(INUM_ZS < 0) THEN
503          WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
504        ENDIF 
505   CASE(1,2) ! aladin et aladin reunion
506       CALL SEARCH_FIELD(6,-1,-1,-1,IGRIB,INUM_ZS)
507        IF(INUM_ZS < 0) THEN
508          WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
509        ENDIF 
510 END SELECT
511 ZPARAM(:)=-999.
512 CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ,IRET_GRIB)
513 ALLOCATE(IINLO(INJ))
514 CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM_ZS),IIU,IJU,ZLONOUT,ZLATOUT,&
515         ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
516 ALLOCATE(ZPARAM_ZS(SIZE(ZPARAM)))
517 ZPARAM_ZS=ZPARAM
518 ALLOCATE(IINLO_ZS(SIZE(IINLO)))
519 IINLO_ZS=IINLO
520 CALL GRIB_GET_SIZE(IGRIB(INUM_ZS),'values',ISIZE)
521 ALLOCATE(ZVALUE(ISIZE))
522 CALL GRIB_GET(IGRIB(INUM_ZS),'values',ZVALUE)
523 ALLOCATE(ZOUT(INO))
524 CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
525             ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
526 DEALLOCATE(IINLO)
527 DEALLOCATE(ZVALUE)
528 !
529 ! Datas given in archives are multiplied by the gravity acceleration
530 ZOUT = ZOUT / XG
531 !
532 ! Stores the field in a 2 dimension array
533 IF (HFILE(1:3)=='ATM') THEN
534   ALLOCATE (XZS_LS(IIU,IJU))
535   ALLOCATE (XZSMT_LS(IIU,IJU))
536   CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XZS_LS)
537   XZSMT_LS = XZS_LS   ! no smooth orography in this case
538 ELSE IF (HFILE=='CHEM') THEN
539   ALLOCATE (XZS_SV_LS(IIU,IJU))
540   CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XZS_SV_LS)
541 END IF
542 DEALLOCATE (ZOUT)
543 !
544 !---------------------------------------------------------------------------------------
545 !* 2.4 Interpolation surface pressure
546 !---------------------------------------------------------------------------------------
547 !
548 !*  2.4.1  Read pressure
549 !   
550 WRITE (ILUOUT0,'(A)') ' | Searching pressure'
551
552 SELECT CASE (IMODEL)
553   CASE(0) ! ECMWF
554      CALL SEARCH_FIELD(152,-1,-1,-1,IGRIB,INUM)
555   CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
556       CALL SEARCH_FIELD(1,-1,-1,-1,IGRIB,INUM)
557 END SELECT
558 IF(INUM < 0) THEN
559    WRITE (ILUOUT0,'(A)')'Surface pressure is missing - abort'
560    CALL ABORT
561    STOP
562 ENDIF 
563 ! recuperation du tableau de valeurs
564 CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
565 ALLOCATE(ZVALUE(ISIZE))
566 CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
567 ! determination des tableaux ZPS_G et ZLNPS_G
568 SELECT CASE (IMODEL)
569   CASE(0) ! ECMWF
570     ALLOCATE (ZPS_G  (ISIZE))
571     ALLOCATE (ZLNPS_G(ISIZE))
572     ZLNPS_G(:) =     ZVALUE(1:ISIZE)
573     ZPS_G  (:) = EXP(ZVALUE(1:ISIZE))
574   CASE(1,2,3,4,5) ! arpege mocage aladin aladin-reunion
575     ALLOCATE (ZPS_G  (ISIZE))
576     ALLOCATE (ZLNPS_G(ISIZE))
577     ZPS_G  (:) =     ZVALUE(1:ISIZE)
578     ZLNPS_G(:) = LOG(ZVALUE(1:ISIZE)) 
579 END SELECT
580 DEALLOCATE (ZVALUE)
581 !
582 !*  2.4.2  Removes pressure points not on orography points
583 !       (if pressure is on a regular grid)
584 CALL GRIB_GET(IGRIB(INUM),'typeOfGrid',HGRID)
585 CALL GRIB_GET(IGRIB(INUM_ZS),'typeOfGrid',HGRID_ZS)
586 CALL GRIB_GET(IGRIB(INUM),'Nj',INJ)
587 CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ_ZS)
588 !
589 IF ( HGRID(1:7)=='regular' .AND. HGRID_ZS(1:7)=='reduced' .AND.&
590      INJ == INJ_ZS) THEN
591    WRITE (ILUOUT0,'(A)')'HGRID(1:7)==regular .AND. HGRID_ZS(1:7)==reduced .AND. INJ == INJ_ZS - abort'
592    CALL ABORT
593    STOP
594 ELSE
595    ALLOCATE(ZWORK_LNPS(SIZE(ZLNPS_G)))
596    ZWORK_LNPS(:) = ZLNPS_G(:)
597 ENDIF   
598 !
599 IF (HFILE(1:3)=='ATM') THEN
600   ALLOCATE (XPS_LS(IIU,IJU))
601 ELSE IF (HFILE=='CHEM') THEN
602   ALLOCATE (XPS_SV_LS(IIU,IJU))
603 END IF
604 !
605 ALLOCATE(IINLO(INJ))
606 CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
607         ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
608 ALLOCATE(ZOUT(INO))
609 CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI,&
610          ZWORK_LNPS,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE. )
611 DEALLOCATE(ZWORK_LNPS)
612 DEALLOCATE(IINLO)
613 !
614 !* 2.4.3 conversion to surface pressure
615 !
616 ZOUT=EXP(ZOUT)
617 IF (HFILE(1:3)=='ATM') THEN
618   CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XPS_LS(:,:))
619 ELSE IF (HFILE=='CHEM') THEN
620   CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XPS_SV_LS(:,:))
621 END IF
622 !JUAN REALZ
623 CALL MPPDB_CHECK2D(XZS_LS,"XZS_LS",PRECISION)
624 !JUAN REALZ
625 DEALLOCATE (ZOUT)
626 DEALLOCATE (ZLNPS_G)
627 !
628 !---------------------------------------------------------------------------------------
629 !* 2.5 Interpolation temperature and humidity
630 !---------------------------------------------------------------------------------------
631 !
632 !* 2.5.1 Read T and Q on each level
633 !
634 WRITE (ILUOUT0,'(A)') ' | Reading T and Q fields'
635 !
636 SELECT CASE (IMODEL)
637   CASE(0) ! ECMWF
638           ISTARTLEVEL=1
639           IT=130
640           IQ=133
641      CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
642      IF(INUM < 0) THEN 
643              ISTARTLEVEL=0
644        CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
645      ENDIF
646      IF(INUM < 0) THEN
647         WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
648         CALL ABORT
649         STOP
650      ENDIF      
651      CALL SEARCH_FIELD(IQ,109,ISTARTLEVEL,-1,IGRIB,INUM)
652      IF(INUM < 0) THEN
653         WRITE (ILUOUT0,'(A)')'Atmospheric specific humidity is missing - abort'
654         CALL ABORT
655         STOP
656      ENDIF 
657   CASE(1,2,3,4,5) ! arpege mocage aladin et aladin reunion
658           IT=11
659           IQ=51
660           ISTARTLEVEL=1
661       CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
662       IF(INUM < 0) THEN 
663           ISTARTLEVEL=0    
664        CALL SEARCH_FIELD(IT,109,ISTARTLEVEL,-1,IGRIB,INUM)
665      ENDIF
666      IF(INUM < 0) THEN
667         WRITE (ILUOUT0,'(A)')'Air temperature is missing - abort'
668         CALL ABORT
669         STOP
670      ENDIF 
671      CALL SEARCH_FIELD(IQ,109,ISTARTLEVEL,-1,IGRIB,INUM)
672      IF(INUM < 0) THEN
673         WRITE (ILUOUT0,'(A)')'Atmospheric specific humidity is missing - abort'
674         CALL ABORT
675         STOP
676      ENDIF 
677 END SELECT
678 !
679 CALL GRIB_GET(IGRIB(INUM),'NV',INLEVEL)
680 INLEVEL = NINT(INLEVEL / 2.) - 1
681 CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
682 !
683 ALLOCATE (ZT_G(ISIZE,INLEVEL))
684 ALLOCATE (ZQ_G(ISIZE,INLEVEL))
685 !
686 DO JLOOP1=1, INLEVEL
687   ILEV1 = JLOOP1-1+ISTARTLEVEL
688   CALL SEARCH_FIELD(IQ,109,ILEV1,-1,IGRIB,INUM)
689   IF (INUM< 0) THEN
690     WRITE (ILUOUT0,'(A,I2,A)') ' -> Atmospheric humidity level ',JLOOP1,' is missing - abort'
691 !callabortstop
692     CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
693     CALL ABORT
694     STOP
695   END IF
696   CALL GRIB_GET(IGRIB(INUM),'values',ZQ_G(:,INLEVEL-JLOOP1+1))
697   CALL SEARCH_FIELD(IT,109,ILEV1,-1,IGRIB,INUM)
698   IF (INUM< 0) THEN
699     WRITE (ILUOUT0,'(A,I2,A)') ' -> Atmospheric temperature level ',JLOOP1,' is missing - abort'
700     !callabortstop
701     CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
702     CALL ABORT
703     STOP
704   END IF
705   CALL GRIB_GET(IGRIB(INUM),'values',ZT_G(:,INLEVEL-JLOOP1+1))
706   CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
707 END DO
708 ALLOCATE(IINLO(INJ))
709 CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
710         ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
711 !
712 !*  2.5.2  Load level definition parameters A and B
713 !
714 IF (HFILE(1:3)=='ATM') THEN
715   XP00_LS = 101325.
716 ELSE IF (HFILE=='CHEM') THEN
717   XP00_SV_LS = 101325.
718 END IF
719 !
720 IF (INLEVEL > 0) THEN
721   IF (HFILE(1:3)=='ATM') THEN
722     ALLOCATE (XA_LS(INLEVEL))
723     ALLOCATE (XB_LS(INLEVEL))
724   ELSE IF (HFILE=='CHEM') THEN
725     ALLOCATE (XA_SV_LS(INLEVEL))
726     ALLOCATE (XB_SV_LS(INLEVEL))
727   END IF
728 !
729   CALL GRIB_GET(IGRIB(INUM),'PVPresent',IPVPRESENT)
730   IF (IPVPRESENT==1) THEN
731      CALL GRIB_GET_SIZE(IGRIB(INUM),'pv',IPV)
732      ALLOCATE(ZPV(IPV))
733      CALL GRIB_GET(IGRIB(INUM),'pv',ZPV)
734   ELSE
735      WRITE (ILUOUT0,'(A)') "THERE IS NO PV VALUE IN THIS MESSAGE : abort"
736      !callabortstop
737      CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
738      CALL ABORT
739      STOP
740   ENDIF
741   SELECT CASE (IMODEL)
742     CASE (0,3,4)
743       DO JLOOP1 = 1, INLEVEL
744         XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+JLOOP1) / XP00_LS
745         XB_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
746       END DO
747     CASE (1,2)
748       JLOOP2 = 2
749       DO JLOOP1 = 1, INLEVEL
750         JLOOP2 = JLOOP2 + 1
751         XA_LS(1 + INLEVEL - JLOOP1) = ZPV(JLOOP2)
752         JLOOP2 = JLOOP2 + 1
753         XB_LS(1 + INLEVEL - JLOOP1) = ZPV(JLOOP2)
754       END DO
755     CASE (5)
756       DO JLOOP1 = 1, INLEVEL
757         IF (HFILE(1:3)=='ATM') THEN
758           XA_LS(1 + INLEVEL - JLOOP1) = ZPV(1+        JLOOP1) / XP00_LS**2 
759           XB_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
760         ELSE IF (HFILE=='CHEM') THEN
761           XA_SV_LS(1 + INLEVEL - JLOOP1) = ZPV(1+        JLOOP1) / XP00_LS**2 
762           XB_SV_LS(1 + INLEVEL - JLOOP1) = ZPV(2+INLEVEL+JLOOP1)
763         END IF
764       END DO
765   END SELECT
766 ELSE
767   WRITE (ILUOUT0,'(A)') ' -> Level definition section is missing - abort'
768  !callabortstop
769   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
770   CALL ABORT
771   STOP
772 END IF
773 !
774 !*  2.5.3  Compute atmospheric pressure on grib grid
775 !
776 WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on Grib grid is being computed'
777 ALLOCATE (ZPF_G(INI,INLEVEL))
778 IF (HFILE(1:3)=='ATM') THEN
779     ZPF_G(:,:) = SPREAD(XA_LS,1,INI)*XP00_LS + &
780     SPREAD(XB_LS,1,INI)*SPREAD(ZPS_G(1:INI),2,INLEVEL)
781 ELSE IF (HFILE=='CHEM') THEN
782     ZPF_G(:,:) = SPREAD(XA_SV_LS,1,INI)*XP00_SV_LS + &
783     SPREAD(XB_SV_LS,1,INI)*SPREAD(ZPS_G(1:INI),2,INLEVEL)
784 END IF
785 DEALLOCATE (ZPS_G)
786 !
787 ALLOCATE (ZEXNF_G(INI,INLEVEL))
788 ZEXNF_G(:,:) = (ZPF_G(:,:)/XP00)**(XRD/XCPD)
789 !
790 ALLOCATE (ZEXNM_G(INI,INLEVEL))
791 ZEXNM_G(:,1:INLEVEL-1) = (ZEXNF_G(:,1:INLEVEL-1)-ZEXNF_G(:,2:INLEVEL)) / &
792                      (LOG(ZEXNF_G(:,1:INLEVEL-1))-LOG(ZEXNF_G(:,2:INLEVEL)))
793 ZEXNM_G(:,INLEVEL) = (ZPF_G(:,INLEVEL)/2./XP00)**(XRD/XCPD)
794 DEALLOCATE (ZEXNF_G)
795 DEALLOCATE (ZPF_G)
796 !
797 ALLOCATE (ZPM_G(INI,INLEVEL))
798 ZPM_G(:,:) = XP00*(ZEXNM_G(:,:))**(XCPD/XRD)
799 !
800 !*  2.5.4  Compute the relative humidity and the interpolate Thetav, P, Q and H
801 !
802 IF (IMODEL==1) THEN
803   ! search cloud_water in Arome case (forecast)
804   ISTARTLEVEL = 1
805   IPAR=246
806   CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
807   IF (INUM < 0) THEN
808     ISTARTLEVEL = 0
809     CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
810   END IF
811   IF (INUM > 0) THEN
812     WRITE (ILUOUT0,'(A)') ' | Grib file from French Weather Service - Arome model (forecast)'
813     LCPL_AROME=.TRUE.
814     NRR=6
815   END IF
816   ! search also turbulent kinetic energy 
817   ISTARTLEVEL = 1
818   IPAR=251
819   CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
820   IF (INUM < 0) THEN
821     ISTARTLEVEL = 0
822     CALL SEARCH_FIELD(IPAR,109,ISTARTLEVEL,-1,IGRIB,INUM)
823   END IF
824   IF (INUM > 0) CTURB='TKEL'
825 END IF
826 !
827 !
828 WRITE (ILUOUT0,'(A)') ' | Computing relative humidity on each level'
829 ALLOCATE (ZH_G(INI))
830 ALLOCATE (ZH_LS(IIU,IJU,INLEVEL))
831 IF (HFILE(1:3)=='ATM') THEN
832   ALLOCATE (XT_LS(IIU,IJU,INLEVEL))
833   ALLOCATE (XQ_LS(IIU,IJU,INLEVEL,NRR)) ; XQ_LS=0.
834 ELSE IF (HFILE=='CHEM') THEN
835   ALLOCATE (XT_SV_LS(IIU,IJU,INLEVEL))
836   ALLOCATE (XQ_SV_LS(IIU,IJU,INLEVEL,1))
837 END IF
838 IF (CTURB=='TKEL') THEN
839   IF (ALLOCATED(XTKE_LS)) DEALLOCATE(XTKE_LS)
840   ALLOCATE(XTKE_LS(IIU,IJU,INLEVEL)) ; XTKE_LS=0.
841   CALL INI_CTURB
842 ELSE
843   IF (ALLOCATED(XTKE_LS)) DEALLOCATE(XTKE_LS)
844   ALLOCATE(XTKE_LS(0,0,0))
845 END IF
846 ALLOCATE (ZTHV_LS(IIU,IJU,INLEVEL))
847 ALLOCATE (ZTHV_G(INI))
848 ALLOCATE (ZRV_G(INI))
849 ALLOCATE (ZOUT(INO))
850 DO JLOOP1=1, INLEVEL
851   !
852   ! Compute Theta V and relative humidity on grib grid
853   !
854   !   (1/rv) = (1/q)  -  1
855   !   Thetav = T . (P0/P)^(Rd/Cpd) . ( (1 + (Rv/Rd).rv) / (1 + rv) )
856   !   Hu = P / ( ( (Rd/Rv) . ((1/rv) - 1) + 1) . Es(T) )
857   !
858   ZRV_G(:) = 1. / (1./MAX(ZQ_G(:,JLOOP1),1.E-12) - 1.)
859   !
860   ZTHV_G(:)=ZT_G(:,JLOOP1) * ((XP00/ZPM_G(:,JLOOP1))**(XRD/XCPD)) * &
861                              ((1. + XRV*ZRV_G(:)/XRD) / (1. + ZRV_G(:)))
862   !
863   ZH_G(1:INI) = 100. * ZPM_G(:,JLOOP1) / ( (XRD/XRV)*(1./ZRV_G(:)+1.)*SM_FOES(ZT_G(:,JLOOP1)) )
864   ZH_G(:) = MAX(MIN(ZH_G(:),100.),0.)
865   !
866   ! Interpolation : H           
867   CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
868         ZH_G,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
869   CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZH_LS(:,:,JLOOP1))
870   ZH_LS(:,:,JLOOP1) = MAX(MIN(ZH_LS(:,:,JLOOP1),100.),0.)
871   !
872   ! interpolation : Theta V
873   CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
874         ZTHV_G,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
875   CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,ZTHV_LS(:,:,JLOOP1))
876   !
877 END DO
878 DEALLOCATE (ZOUT)
879 !
880 !*  2.5.5  Compute atmospheric pressure on MESO-NH grid
881 !
882 WRITE (ILUOUT0,'(A)') ' | Atmospheric pressure on MesoNH grid is being computed'
883 ALLOCATE (ZPF_LS(IIU,IJU,INLEVEL))
884 IF (HFILE(1:3)=='ATM') THEN
885   ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_LS,1,IIU),2,IJU)*XP00_LS + &
886   SPREAD(SPREAD(XB_LS,1,IIU),2,IJU)*SPREAD(XPS_LS,3,INLEVEL)
887 ELSE IF (HFILE=='CHEM') THEN
888   ZPF_LS(:,:,:) = SPREAD(SPREAD(XA_SV_LS,1,IIU),2,IJU)*XP00_LS + &
889   SPREAD(SPREAD(XB_SV_LS,1,IIU),2,IJU)*SPREAD(XPS_SV_LS,3,INLEVEL)
890 END IF
891 !
892 ALLOCATE (ZEXNF_LS(IIU,IJU,INLEVEL))
893 ZEXNF_LS(:,:,:) = (ZPF_LS(:,:,:)/XP00)**(XRD/XCPD)
894 !
895 ALLOCATE (ZEXNM_LS(IIU,IJU,INLEVEL))
896 ZEXNM_LS(:,:,1:INLEVEL-1) = (ZEXNF_LS(:,:,1:INLEVEL-1)-ZEXNF_LS(:,:,2:INLEVEL)) / &
897                      (LOG(ZEXNF_LS(:,:,1:INLEVEL-1))-LOG(ZEXNF_LS(:,:,2:INLEVEL)))
898 ZEXNM_LS(:,:,INLEVEL) = (ZPF_LS(:,:,INLEVEL)/2./XP00)**(XRD/XCPD)
899 DEALLOCATE (ZEXNF_LS)
900 DEALLOCATE (ZPF_LS)
901 !
902 ALLOCATE (ZPM_LS(IIU,IJU,INLEVEL))
903 ZPM_LS(:,:,:) = XP00*(ZEXNM_LS(:,:,:))**(XCPD/XRD)
904 !
905 !*  2.5.6 Compute the vapor mixing ratio and the final specific humdity
906 !
907 !  The vapor mixing ratio is calculated by an interating process on rv and
908 !  Thetav. Have a look to MODE_THERMO for further informations.
909 ALLOCATE (ZR_DUM(IIU,IJU,INLEVEL,1))
910 ALLOCATE (ZRV_LS(IIU,IJU,INLEVEL))
911 ALLOCATE (ZTEV_LS(IIU,IJU,INLEVEL))
912 ZTEV_LS(:,:,:) = ZTHV_LS(:,:,:) * ZEXNM_LS(:,:,:)
913 ZRV_LS(:,:,:) = SM_PMR_HU(CLUOUT0, ZPM_LS(:,:,:),      &
914 ZTEV_LS(:,:,:),ZH_LS(:,:,:),ZR_DUM(:,:,:,:),KITERMAX=100)
915 IF (HFILE(1:3)=='ATM') THEN
916   XQ_LS(:,:,:,1) = ZRV_LS(:,:,:) / (1. + ZRV_LS(:,:,:))
917 ELSE IF (HFILE=='CHEM') THEN
918   XQ_SV_LS(:,:,:,1) = ZRV_LS(:,:,:) / (1. + ZRV_LS(:,:,:))
919 ENDIF
920 !JUAN
921 CALL  MPPDB_CHECK3D(XQ_LS(:,:,:,1),"XQ_LS",PRECISION)
922 !JUAN
923 DEALLOCATE (ZTEV_LS)
924 DEALLOCATE (ZH_LS)
925 DEALLOCATE (ZR_DUM)
926 !
927 !*  2.5.7 Compute T from the interpolated Theta V
928 !
929 !   T = Thetav . (P/P0)^(Rd/Cpd) . ((1 + rv) / (1 + (Rv/Rd).rv))
930 !!
931 IF (HFILE(1:3)=='ATM') THEN
932   XT_LS(:,:,:) = ZTHV_LS(:,:,:) * ZEXNM_LS(:,:,:) &
933                               * ((1.+ZRV_LS(:,:,:))/(1.+(XRV/XRD)*ZRV_LS(:,:,:)))
934  CALL MPPDB_CHECK3D(XT_LS,"XT_LS",PRECISION)
935 ELSE IF (HFILE=='CHEM') THEN
936   XT_SV_LS(:,:,:) = ZTHV_LS(:,:,:) * ZEXNM_LS(:,:,:) &
937                               * ((1.+ZRV_LS(:,:,:))/(1.+(XRV/XRD)*ZRV_LS(:,:,:)))
938  CALL MPPDB_CHECK3D(XT_SV_LS,"XT_SV_LS",PRECISION)
939 ENDIF
940 !
941 DEALLOCATE (ZRV_LS)
942 DEALLOCATE (ZTHV_LS)
943 DEALLOCATE (ZPM_LS)
944 DEALLOCATE (ZEXNM_LS)
945 !
946 !*  2.5.8 Read the other specific ratios (if Arome model)
947 !
948 IF (NRR >1) THEN
949   WRITE (ILUOUT0,'(A)') ' | Reading Q fields (except humidity)'
950   DO JLOOP2=1,NRR-1
951     IPAR=246+JLOOP2-1
952     DO JLOOP1=1, INLEVEL
953       ILEV1 = JLOOP1-1+ISTARTLEVEL
954       CALL SEARCH_FIELD(IPAR,109,ILEV1,-1,IGRIB,INUM)
955
956       IF (INUM < 0) THEN
957         WRITE (ILUOUT0,'(A,I3,A,I2,A)') ' -> Specific ratio ',IPAR,' at level ',JLOOP1,' is missing - abort'
958         !callabortstop
959         CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
960         CALL ABORT
961         STOP
962       END IF
963       CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
964       ALLOCATE(ZVALUE(ISIZE))
965       CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
966       ALLOCATE(ZOUT(INO))
967       CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
968         ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
969       CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XQ_LS(:,:,INLEVEL-JLOOP1+1,1+JLOOP2))
970       DEALLOCATE (ZVALUE)
971       DEALLOCATE (ZOUT)
972     END DO
973   END DO
974 END IF
975 !
976 IF (CTURB=='TKEL') THEN
977   WRITE (ILUOUT0,'(A)') ' | Reading TKE field'
978   IPAR=251
979   DO JLOOP1=1, INLEVEL
980     ILEV1 = JLOOP1-1+ISTARTLEVEL
981     CALL SEARCH_FIELD(IPAR,109,ILEV1,-1,IGRIB,INUM)
982     IF (INUM <  0) THEN
983       WRITE (ILUOUT0,'(A,I3,A,I2,A)') ' -> TKE at level ',JLOOP1,' is missing - abort'
984       !callabortstop
985       CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
986       CALL ABORT
987       STOP
988     END IF
989     CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
990     ALLOCATE(ZVALUE(ISIZE))
991     CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
992     ALLOCATE(ZOUT(INO))
993     CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
994         ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
995     CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XTKE_LS(:,:,INLEVEL-JLOOP1+1))
996     DEALLOCATE (ZVALUE)
997     DEALLOCATE (ZOUT)
998   END DO
999 END IF
1000 DEALLOCATE(IINLO)
1001 !
1002 !---------------------------------------------------------------------------------------
1003 !* 2.6 Interpolation of MOCAGE variable
1004 !---------------------------------------------------------------------------------------
1005
1006 IF (IMODEL==5) THEN 
1007   ! Always initialize chemical scheme variables before INI_NSV call !
1008   CALL CH_INIT_SCHEME_n(IMI,LUSECHAQ,LUSECHIC,LCH_PH,ILUOUT0,KVERB)
1009   LUSECHEM = .TRUE.
1010   IF (LORILAM) THEN
1011     CORGANIC = "MPMPO"
1012     LVARSIGI = .TRUE.
1013     LVARSIGJ = .TRUE.
1014     CALL CH_AER_INIT_SOA(ILUOUT0, KVERB)
1015   END IF
1016   ! initialise NSV_* variables
1017   CALL INI_NSV(1)
1018   IF( HFILE=='ATM0' ) THEN
1019     ALLOCATE (XSV_LS(IIU,IJU,INLEVEL,NSV))
1020   ELSE IF (HFILE=='CHEM' ) THEN
1021     DEALLOCATE(XSV_LS)
1022     ALLOCATE (XSV_LS(IIU,IJU,INLEVEL,NSV))
1023   ELSE
1024     WRITE (ILUOUT0,'(A)') ' -> Mocage model: Bad input argument in read_all_data_grib_case - abort'
1025     !callabortstop
1026     CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
1027     CALL ABORT
1028     STOP
1029   END IF
1030   XSV_LS(:,:,:,:) = 0.
1031   ITYP=109
1032   ILEV1=-1
1033   ILEV2=-1
1034 !
1035   WRITE (ILUOUT0,'(A,A4,A)') ' | Reading Mocage species (ppp) from ',HFILE,' file'
1036 !
1037 !*       2.6.1  read mocage species
1038 !
1039 ! open input file
1040   CALL CH_OPEN_INPUT(YPRE_MOC, "MOC2MESONH", ICHANNEL, ILUOUT0, KVERB)
1041 !
1042 ! read number of mocage species to transfer into mesonh
1043   READ(ICHANNEL, *) IMOC
1044   IF (KVERB >= 5) WRITE(ILUOUT0,*) "number of mocage species to transfer into mesonh : ", IMOC
1045 !
1046 ! read data input format
1047   READ(ICHANNEL,"(A)") YFORMAT
1048   YFORMAT=UPCASE(YFORMAT)
1049   IF (KVERB >= 5) WRITE(ILUOUT0,*) "input format is: ", YFORMAT
1050 !
1051 ! allocate fields
1052   ALLOCATE(YMNHNAME(IMOC))
1053   ALLOCATE(INUMGRIB(IMOC))
1054 !
1055 ! read variables names and Grib code
1056   IF (INDEX(YFORMAT,'A') < INDEX(YFORMAT,'I')) THEN
1057     DO JI = 1, IMOC
1058       READ(ICHANNEL,YFORMAT) YMNHNAME(JI), INUMGRIB(JI)
1059       WRITE(ILUOUT0,YFORMAT) YMNHNAME(JI), INUMGRIB(JI)
1060     END DO
1061   ELSE
1062     DO JI = 1, IMOC
1063       READ(ICHANNEL,YFORMAT) INUMGRIB(JI), YMNHNAME(JI)
1064       WRITE(ILUOUT0,YFORMAT) INUMGRIB(JI), YMNHNAME(JI)
1065     END DO
1066   ENDIF
1067   !
1068   ! close file
1069   CALL CLOSE_ll(YPRE_MOC,IOSTAT=IRET)
1070   !
1071   !*  2.6.2   exchange mocage values onto prognostic variables XSV_LS
1072   !
1073   IF (KVERB >= 10) WRITE(ILUOUT0,'(A,I4)') ' NEQ=',NEQ
1074   !
1075   DO JNREAL = 1, NEQ
1076     INACT = 0
1077     search_loop2 : DO JN = 1, IMOC
1078       IF (CNAMES(JNREAL) .EQ. YMNHNAME(JN)) THEN
1079         INACT = JN
1080         EXIT search_loop2
1081       END IF
1082     END DO search_loop2
1083     IF (INACT .NE. 0) THEN
1084       DO JLOOP1=1, INLEVEL
1085         ILEV1 = JLOOP1
1086         CALL SEARCH_FIELD(INUMGRIB(JN),ITYP,ILEV1,ILEV2,IGRIB,INUM)
1087         IF (INUM <  0) THEN
1088           WRITE (ILUOUT0,'(A,I3,A,I2,A)') ' -> Atmospheric ',INUMGRIB(JN),&
1089           ' grib chemical species level ',JLOOP1,' is missing - abort'
1090           !callabortstop
1091           CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
1092           CALL ABORT
1093           STOP
1094         END IF
1095         CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
1096         ALLOCATE(IINLO(INJ))
1097         CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
1098           ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
1099         CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
1100         ALLOCATE(ZVALUE(ISIZE))
1101         CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
1102         ALLOCATE(ZOUT(INO))
1103         CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
1104            ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.TRUE. )
1105         CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU, &
1106                             XSV_LS(:,:,INLEVEL-JLOOP1+1,JNREAL))
1107         DEALLOCATE (ZVALUE)
1108         DEALLOCATE (ZOUT)
1109         DEALLOCATE(IINLO)
1110       END DO
1111     END IF
1112   END DO
1113   XSV_LS(:,:,:,:) = MAX(XSV_LS(:,:,:,:),0.)
1114 ELSE
1115   LORILAM  = .FALSE. 
1116   LUSECHEM = .FALSE. 
1117   ! initialise NSV_* variables
1118   CALL INI_NSV(1)
1119   IF (NSV > 0) THEN 
1120     ALLOCATE (XSV_LS(IIU,IJU,INLEVEL,NSV))
1121     XSV_LS(:,:,:,:) = 0.
1122   ELSE 
1123     ALLOCATE(XSV_LS(0,0,0,0))
1124   END IF
1125 END IF
1126 !
1127 !---------------------------------------------------------------------------------------
1128 !* 2.7 Search, read, interpolate and project winds
1129 !---------------------------------------------------------------------------------------
1130 !
1131 ! The way winds are processed depends upon the type of archive :
1132 !
1133 ! -> ECMWF
1134 !   Winds are projected from a standard lat,lon grid to MesoNH grid. This correcponds to
1135 !   a rotation of an angle :
1136 !    Alpha = k.(L-L0) - Beta
1137 !   k,L0 and Beta definiiton parameter of MesoNH grid
1138 !   L longitude of the point of the tangent coordinate system
1139 !
1140 ! -> Aladin
1141 !   The grid used by Aladin files is the same than the one of MesoNH. !
1142 ! So we have 2 sets  of parameters :
1143 !     k   L0   Beta      for MesoNH
1144 !     k'  L0'  Beta'     for Aladin (Beta'=0 for Aladin)
1145 !   We applied twice the formula seen for standard lat,lon grid and we get :
1146 !    Alpha = k.(L-L0) - Beta - k'.(L-L0')
1147 !
1148 ! -> Arpege
1149 ! Arpege winds are given on the tangent coordinate system of the stretched grid.
1150 ! Therefore they have first to be projected on a standard lat,lon grid. This is done
1151 ! before the interpolation. The projection formulas were given by Meteo France.
1152 ! After this projection, the file is simil
1153 !
1154 IF (HFILE(1:3)=='ATM') THEN
1155 ITYP  = 109
1156 ISTARTLEVEL = 1
1157 ILEV2 = -1
1158 ALLOCATE (XU_LS(IIU,IJU,INLEVEL))
1159 ALLOCATE (XV_LS(IIU,IJU,INLEVEL))
1160 ALLOCATE (ZTU_LS(INO))
1161 ALLOCATE (ZTV_LS(INO))
1162 !
1163 SELECT CASE (IMODEL)
1164   CASE (0)
1165     IPAR = 131
1166     CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
1167     IF (INUM< 0) THEN
1168       ISTARTLEVEL = 0
1169       CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
1170     END IF
1171   CASE (1,2,3)
1172     IPAR  = 33
1173     CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
1174     IF (INUM < 0) THEN
1175       ISTARTLEVEL = 0
1176       CALL SEARCH_FIELD(IPAR,ITYP,ISTARTLEVEL,ILEV2,IGRIB,INUM)
1177     END IF
1178 END SELECT
1179
1180 DO JLOOP1 = ISTARTLEVEL, ISTARTLEVEL+INLEVEL-1
1181   ILEV1 = JLOOP1
1182   ! read component u 
1183   CALL SEARCH_FIELD(IPAR,ITYP,ILEV1,ILEV2,IGRIB,INUM)
1184   IF (INUM < 0) THEN
1185     WRITE (ILUOUT0,'(A,I2,A)') ' -> Wind vector component "u", level ', &
1186       JLOOP1,' is missing - abort'
1187     !callabortstop
1188     CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
1189     CALL ABORT
1190     STOP
1191   END IF
1192   CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
1193   ALLOCATE(ZVALUE(ISIZE))
1194   CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
1195   IF (IMODEL==3) THEN
1196     ALLOCATE(ZTU0_LS(INI))
1197     ZTU0_LS(:) = ZVALUE(:)
1198   ELSE
1199     CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
1200     IF(ALLOCATED(IINLO)) DEALLOCATE (IINLO)
1201     ALLOCATE(IINLO(INJ))      
1202     CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
1203           ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
1204     ALLOCATE(ZOUT(INO))
1205     CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
1206            ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.TRUE.,PTIME_HORI,.FALSE. )
1207     ZTU_LS(:) = ZOUT(:)
1208     DEALLOCATE(IINLO)
1209     DEALLOCATE(ZOUT)
1210   END IF
1211   DEALLOCATE (ZVALUE)
1212   ! read component v and perform interpolation if not Arpege grid
1213   ILEV1 = JLOOP1
1214   CALL SEARCH_FIELD(IPAR+1,ITYP,ILEV1,ILEV2,IGRIB,INUM)
1215   IF (INUM < 0) THEN
1216     WRITE (ILUOUT0,'(A,I2,A)') ' -> Wind vector component "v", level ', &
1217       JLOOP1,' is missing - abort'
1218     !callabortstop
1219     CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
1220     CALL ABORT
1221     STOP
1222   END IF
1223   CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
1224   ALLOCATE(ZVALUE(ISIZE))
1225   CALL GRIB_GET(IGRIB(INUM),'values',ZVALUE)
1226   IF ((IMODEL==3)) THEN 
1227     CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
1228     ALLOCATE(IINLO(INJ))         
1229     CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
1230           ZXOUT,ZYOUT,INI,ZPARAM,IINLO)     
1231     ALLOCATE(ZTV0_LS(INI))
1232     ZTV0_LS(:) = ZVALUE(:)
1233   ELSE
1234      CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
1235     ALLOCATE(IINLO(INJ))          
1236     CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
1237           ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
1238     ALLOCATE(ZOUT(INO))
1239     CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
1240            ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.TRUE.,PTIME_HORI,.FALSE.)
1241     ZTV_LS(:) = ZOUT(:)
1242     DEALLOCATE(ZOUT)
1243   END IF
1244   DEALLOCATE (ZVALUE)
1245     ! interpolations for arpege grid
1246   IF ((IMODEL==3)) THEN
1247     ! Comes back to real winds instead of stretched winds
1248     ! (but still with components according to Arpege grid axes)
1249     ZLATPOLE = ZPARAM(7) * XPI/180.
1250     ZLONPOLE = ZPARAM(8) * XPI/180.
1251     ZC       = ZPARAM(9)
1252     ZD       = ZC * ZC
1253     JLOOP3 = 0
1254     JLOOP4 = 1
1255     ZLAT = ZPARAM(3) * XPI / 180.
1256     DO JLOOP2=1, INI
1257       ZLON = JLOOP3 * 2. * XPI / IINLO(JLOOP4)
1258       ! Compute the scale factor
1259       ZA = ((1.+ZD) - (1.-ZD)*SIN(ZLAT)) / (2. * ZC)
1260       ZTU0_LS(JLOOP2) = ZTU0_LS(JLOOP2) * ZA
1261       ZTV0_LS(JLOOP2) = ZTV0_LS(JLOOP2) * ZA
1262       ! next parallel
1263       JLOOP3 = JLOOP3 + 1
1264       IF (JLOOP3 == IINLO(JLOOP4)) THEN
1265         JLOOP3 = 0
1266         ZLAT = ZLAT + (((ZPARAM(5)-ZPARAM(3))/(ZPARAM(2)-1)) * XPI / 180.)
1267         JLOOP4 = JLOOP4 + 1
1268       END IF
1269     END DO
1270     !
1271     ! interpolation
1272     CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,&
1273                 INI,ZTU0_LS,INO,ZXOUT,ZYOUT,ZTU_LS,.TRUE.,PTIME_HORI,.FALSE.)
1274     CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,&
1275                 INI,ZTV0_LS,INO,ZXOUT,ZYOUT,ZTV_LS,.TRUE.,PTIME_HORI,.FALSE.)
1276     DEALLOCATE(IINLO)
1277     !
1278     ! Rotation of the components from Arpege grid axes to real sphere axes
1279     !
1280     DO JLOOP2=1, INO
1281       ZLAT = ZYOUT(JLOOP2) * XPI / 180.
1282       ZLON = ZXOUT(JLOOP2) * XPI / 180.
1283       ! Compute the rotation matrix
1284       ZA = (ZD+1.) + (ZD-1.)*SIN(ZLAT)
1285       ZB = (ZD-1.) + (ZD+1.)*SIN(ZLAT)
1286       ZE = 2.*ZC*COS(ZLATPOLE)*COS(ZLAT)*COS(ZLON) + ZB*SIN(ZLATPOLE)
1287       IF (ABS(ZE) .GE. ABS(ZA)) THEN
1288         ZF = -2.*ZC*COS(ZLATPOLE)/ ( COS(ZLAT)* ((ZD+1.)+(ZD-1.)*SIN(ZLATPOLE)) )
1289         ZSIN = -ZF*SIN(ZLONPOLE-ZLON)
1290         ZCOS =  ZF*COS(ZLONPOLE-ZLON)
1291       ELSE
1292         ZF = 1. / SQRT(ZA*ZA - ZE*ZE)
1293         ZSIN = -COS(ZLATPOLE)*SIN(ZLON)*ZA*ZF
1294         ZCOS = (2.*ZC*SIN(ZLATPOLE)*COS(ZLAT)-ZB*COS(ZLATPOLE)*COS(ZLON))*ZF
1295       ENDIF
1296       ZTEMP = ZTU_LS(JLOOP2)
1297       ZTU_LS(JLOOP2) = ZCOS*ZTEMP - ZSIN*ZTV_LS(JLOOP2)
1298       ZTV_LS(JLOOP2) = ZSIN*ZTEMP + ZCOS*ZTV_LS(JLOOP2)
1299     END DO
1300   END IF
1301   !
1302   ! Rotation of the components from the real sphere axes (Arpege, CEP)
1303   ! or model axes (Aladin) to MESO-NH axes
1304   !
1305   JLOOP4=0
1306   DO JJ=1,IJU
1307     DO JI=1,IIU
1308       JLOOP4=JLOOP4+1
1309       IF (IMODEL==2 .OR. IMODEL==1 ) THEN 
1310           IF (IMODEL==2) THEN          ! ALADIN  REUNION
1311               ZALPHA=0
1312           ELSE  !ALADIN
1313               ZALPHA = (XRPK*(ZLONOUT(JLOOP4)-XLON0)-XBETA) - &
1314               (SIN(ZPARAM(9)*XPI/180.)*(ZLONOUT(JLOOP4)-ZPARAM(10)))
1315           ENDIF
1316       ELSE  ! CEP, ARPEGE (after projection)
1317           ZALPHA = XRPK*(ZLONOUT(JLOOP4)-XLON0)-XBETA
1318       ENDIF
1319       ZALPHA = ZALPHA * XPI / 180.
1320       XU_LS(JI,JJ,INLEVEL+ISTARTLEVEL-JLOOP1)= &
1321           ZTU_LS(JLOOP4)*COS(ZALPHA) - ZTV_LS(JLOOP4)*SIN(ZALPHA)
1322       XV_LS(JI,JJ,INLEVEL+ISTARTLEVEL-JLOOP1)= &
1323           ZTU_LS(JLOOP4)*SIN(ZALPHA) + ZTV_LS(JLOOP4)*COS(ZALPHA)
1324     ENDDO
1325   ENDDO
1326   IF ((IMODEL==3)) THEN ! deallocation of Arpege arrays
1327     DEALLOCATE (ZTU0_LS)
1328     DEALLOCATE (ZTV0_LS)
1329   END IF
1330 END DO
1331 DEALLOCATE (ZTU_LS)
1332 DEALLOCATE (ZTV_LS)
1333 IF(ALLOCATED(IINLO)) DEALLOCATE (IINLO)
1334 END IF
1335 !
1336 !-------------------------------------------------------------------------------
1337 !* 2.8 Filter the characteristics of the large-scale vortex
1338 !-------------------------------------------------------------------------------
1339 IF (HFILE(1:3)=='ATM' .AND. LFILTERING) THEN
1340   WRITE (ILUOUT0,'(A)') ' | Starting the filtering of the fields to remove large-scale vortex'
1341   IF (INDEX(CFILTERING,'Q')/=0) THEN
1342     WRITE (ILUOUT0,'(A)') ' -> Filtering of Q is now available!'       
1343     WRITE (ILUOUT0,'(A,A5)') ' CFILTERING= ',CFILTERING
1344   ENDIF
1345   !
1346   IF (INDEX(CFILTERING,'P')/=0) THEN
1347     ! compute reduced surface pressure
1348     ALLOCATE(ZTVF_LS(IIU,IJU),ZMSLP_LS(IIU,IJU))
1349     ! compute pressure reduced to first level above mean sea level
1350     !                                        (rather than above ground level)
1351     ZGAMREF=-6.5E-3
1352     !virtual temperature at the first level above ground
1353     ZTVF_LS(:,:) = XT_LS(:,:,1)*(1.+XQ_LS(:,:,1,1)*(XRV/XRD-1))
1354     !virtual temperature averaged between first level above ground
1355     !                                 and first level above sea level
1356     ZTVF_LS(:,:) = ZTVF_LS(:,:)-0.5*ZGAMREF*XZS_LS(:,:)
1357     ZMSLP_LS(:,:)=XPS_LS(:,:)*EXP(XG*XZS_LS(:,:)/(XRD*ZTVF_LS(:,:)))
1358   ENDIF
1359   !
1360   IF (INDEX(CFILTERING,'P')==0) THEN
1361     IF (INDEX(CFILTERING,'Q')==0) THEN
1362       CALL REMOVAL_VORTEX(XZS_LS,XU_LS,XV_LS,XT_LS)
1363     ELSE
1364       CALL REMOVAL_VORTEX(XZS_LS,XU_LS,XV_LS,XT_LS,PQ_LS=XQ_LS(:,:,:,1))
1365     ENDIF
1366   ELSE
1367     IF (INDEX(CFILTERING,'Q')==0) THEN
1368       CALL REMOVAL_VORTEX(XZS_LS,XU_LS,XV_LS,XT_LS,PPS_LS=ZMSLP_LS)
1369     ELSE
1370       CALL REMOVAL_VORTEX(XZS_LS,XU_LS,XV_LS,XT_LS,PQ_LS=XQ_LS(:,:,:,1),PPS_LS=ZMSLP_LS)
1371     ENDIF
1372     XPS_LS(:,:)  = ZMSLP_LS(:,:)*EXP(-XG*XZS_LS(:,:)/(XRD*ZTVF_LS(:,:)))
1373     DEALLOCATE(ZTVF_LS,ZMSLP_LS)
1374   ENDIF
1375   !
1376 END IF
1377 !
1378 !---------------------------------------------------------------------------------------
1379 !* 2.9 Read date
1380 !---------------------------------------------------------------------------------------
1381 !
1382 WRITE (ILUOUT0,'(A)') ' | Reading date'
1383 CALL GRIB_GET(IGRIB(INUM),'dataDate',IDATE,IRET_GRIB)
1384 CALL GRIB_GET(IGRIB(INUM),'dataTime',ITIME,IRET_GRIB)
1385 TPTCUR%TIME=ITIME/100*3600+(ITIME-(ITIME/100)*100)*60
1386 TPTCUR%TDATE%YEAR=IDATE/10000
1387 TPTCUR%TDATE%MONTH=INT((IDATE-TPTCUR%TDATE%YEAR*10000)/100)
1388 TPTCUR%TDATE%DAY=IDATE-TPTCUR%TDATE%YEAR*10000-TPTCUR%TDATE%MONTH*100
1389 CALL GRIB_GET(IGRIB(INUM),'startStep',ITIMESTEP,IRET_GRIB)
1390 CALL GRIB_GET(IGRIB(INUM),'stepUnits',CSTEPUNIT,IRET_GRIB)
1391 IF (IMODEL==0) THEN
1392   ITWOZS=0
1393   IF ((TPTCUR%TDATE%YEAR ==2000).AND.(TPTCUR%TDATE%MONTH  >11)) ITWOZS=1
1394   IF ((TPTCUR%TDATE%YEAR ==2000).AND.(TPTCUR%TDATE%MONTH ==11)) THEN
1395     IF ( (TPTCUR%TDATE%DAY   >20 )   .OR.   &
1396         ((TPTCUR%TDATE%DAY  ==20 ).AND.(TPTCUR%TIME >=64800 ))) ITWOZS=1
1397   END IF
1398   IF ( TPTCUR%TDATE%YEAR ==2001 )  ITWOZS=1
1399   IF ((TPTCUR%TDATE%YEAR ==2002).AND.(TPTCUR%TDATE%MONTH  <11)) ITWOZS=1
1400   IF ((TPTCUR%TDATE%YEAR ==2002).AND.(TPTCUR%TDATE%MONTH ==11)) THEN
1401     IF ( (TPTCUR%TDATE%DAY   <24 )   .OR.   &
1402         ((TPTCUR%TDATE%DAY  ==25 ).AND.(TPTCUR%TIME  <64800 ))) ITWOZS=1
1403   END IF
1404   IF (ITWOZS==1) &
1405     WRITE(ILUOUT0,*) ' Check that both orography fields on 1st model level and on surface are used.'
1406 END IF
1407
1408 CALL MPPDB_CHECK3D(XU_LS,"XU_LS",PRECISION)
1409 CALL MPPDB_CHECK3D(XV_LS,"XV_LS",PRECISION)
1410
1411 SELECT CASE (CSTEPUNIT)       ! Time unit indicator
1412   CASE ('h')                    !hour
1413     TPTCUR%TIME   = TPTCUR%TIME + ITIMESTEP*3600.
1414   CASE ('m')                    !minute
1415     TPTCUR%TIME   = TPTCUR%TIME + ITIMESTEP*60.
1416   CASE ('s')                    !minute
1417     TPTCUR%TIME   = TPTCUR%TIME + ITIMESTEP
1418    CASE DEFAULT
1419     WRITE (ILUOUT0,'(A,A,A)') ' | error CSTEPUNIT=',CSTEPUNIT, ' is different of s,m or h'
1420 END SELECT
1421 CALL ADD_FORECAST_TO_DATE(TPTCUR%TDATE%YEAR, &
1422                           TPTCUR%TDATE%MONTH,&
1423                           TPTCUR%TDATE%DAY,  &
1424                           TPTCUR%TIME        )
1425 IF (HFILE(1:3)=='ATM') THEN
1426   CALL SM_PRINT_TIME(TPTCUR,CLUOUT0,'MESONH current date')
1427   TDTCUR = TPTCUR
1428   TDTMOD = TPTCUR
1429   TDTSEG = TPTCUR
1430   TDTEXP = TPTCUR
1431 ELSE IF (HFILE=='CHEM') THEN
1432   CALL SM_PRINT_TIME(TPTCUR,CLUOUT0,'current date in MesoNH format')
1433 ENDIF
1434 !
1435 !-------------------------------------------------------------------------------
1436 !* 2.10 Read and interpolate dummy fields listed in free-format part of nml file
1437 !-------------------------------------------------------------------------------
1438 IF (ODUMMY_REAL) THEN
1439   !
1440   WRITE (ILUOUT0,'(A)') ' | Try to read 2D dummy fields'
1441   !
1442   !*       2.10.1   read 2D dummy fields
1443   !
1444   ! close file
1445   CALL CLOSE_ll(HPRE_REAL1,IOSTAT=IRET)
1446   ! open input file
1447   CALL CH_OPEN_INPUT(HPRE_REAL1, "DUMMY_2D", ICHANNEL, ILUOUT0, KVERB)
1448   !
1449   ! read number of dummy 2D fields to transfer into mesonh
1450   READ(ICHANNEL, *) IMOC
1451   IF (KVERB >= 5) WRITE(ILUOUT0,*) "number of dummy fields to transfer into Mesonh : ", IMOC
1452   ALLOCATE(XDUMMY_2D(IIU,IJU,IMOC),CDUMMY_2D(IMOC))
1453   ALLOCATE(INUMGRIB(IMOC),INUMLEV(IMOC),INUMLEV1(IMOC),INUMLEV2(IMOC))
1454   INUMLEV(:)=-1 ; INUMLEV1(:)=-1 ; INUMLEV2(:)=-1
1455   !
1456   IVAR=0
1457   ! read variables names and Grib codes
1458   DO JI = 1, IMOC
1459     READ(ICHANNEL,'(A)') YINPLINE
1460     YINPLINE= TRIM(ADJUSTL(YINPLINE))
1461     IF (LEN_TRIM(YINPLINE) == 0) CYCLE ! skip blank line
1462     ! transform tab and comma character into blank
1463     DO JJ=1,LEN_TRIM(YINPLINE)
1464       IF (YINPLINE(JJ:JJ)==YPTAB .OR. YINPLINE(JJ:JJ)==YPCOM) YINPLINE(JJ:JJ)= ' '
1465     END DO
1466     IF (KVERB >= 10) WRITE(ILUOUT0,*) 'Line read : ', YINPLINE
1467     ! extract field name
1468     INDX= INDEX(YINPLINE,' ')
1469     YFIELD= YINPLINE(1:INDX-1)
1470     IF (KVERB >= 5) WRITE(ILUOUT0,*) 'Field being treated : ', YFIELD
1471     ITYP=105
1472     ILEV1=-1
1473     ILEV2=-1
1474     ! extract the parameter indicator
1475     YINPLINE= ADJUSTL(YINPLINE(INDX:))
1476     INDX= INDEX(YINPLINE,' ')
1477     IF (INDX == 1) THEN
1478       WRITE(ILUOUT0,*) ' Parameter indicator is missing. ',YFIELD,' not treated.'
1479       CYCLE
1480     END IF
1481     IVAR=IVAR+1
1482     READ(YINPLINE(1:INDX-1),*) IPAR
1483     IF (NVERB>=5) WRITE(ILUOUT0,*) ' Parameter indicator: ',IPAR
1484     ! extract the level indicator (optional)
1485     YINPLINE= ADJUSTL(YINPLINE(INDX:))
1486     INDX= INDEX(YINPLINE,' ')
1487     IF (INDX /= 1) THEN
1488       READ(YINPLINE(1:INDX-1),*) ITYP
1489       IF (NVERB>=5) WRITE(ILUOUT0,*) ' Level indicator is indicated: ',ITYP 
1490     END IF
1491     ! extract the first level value (optional)
1492     YINPLINE= ADJUSTL(YINPLINE(INDX:))
1493     INDX= INDEX(YINPLINE,' ')
1494     IF (INDX /= 1) THEN
1495       READ(YINPLINE(1:INDX-1),*) ILEV1
1496       IF (NVERB>=5) WRITE(ILUOUT0,*) ' Level1 value is indicated: ',ILEV1
1497     END IF
1498     ! extract the second level value (optional)
1499     YINPLINE= ADJUSTL(YINPLINE(INDX:))
1500     INDX= INDEX(YINPLINE,' ')
1501     IF (INDX /= 1) THEN
1502       READ(YINPLINE(1:INDX-1),*) ILEV2
1503       IF (NVERB>=5) WRITE(ILUOUT0,*) ' Level2 value is indicated: ',ILEV2
1504     END IF
1505     !
1506     CDUMMY_2D(IVAR)=YFIELD ; INUMGRIB(IVAR)=IPAR
1507     INUMLEV(IVAR)=ITYP     ; INUMLEV1(IVAR)=ILEV1 ; INUMLEV2(IVAR)=ILEV2
1508     !
1509   END DO
1510   !
1511   IF (NVERB>=10) THEN
1512     WRITE(ILUOUT0,*) CDUMMY_2D(1:IVAR)
1513     WRITE(ILUOUT0,*) INUMGRIB(1:IVAR)
1514     WRITE(ILUOUT0,*) INUMLEV(1:IVAR)
1515     WRITE(ILUOUT0,*) INUMLEV1(1:IVAR)
1516     WRITE(ILUOUT0,*) INUMLEV2(1:IVAR)
1517   END IF
1518   !
1519   IF (IVAR /= IMOC) THEN
1520     WRITE (ILUOUT0,'(A,I3,A,I3,A)') ' -> Number of correct lines (',IVAR,') is different of ',IMOC,' - abort'
1521    !callabortstop
1522     CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
1523     CALL ABORT
1524     STOP
1525   END IF
1526   !
1527   !* 2.10.2 read and interpolate variables onto dummy variables XDUMMY_2D
1528   !
1529   DO JI = 1, IMOC
1530     WRITE(ILUOUT0,'(A,4(I3,1X))') CDUMMY_2D(JI),INUMGRIB(JI),INUMLEV(JI),INUMLEV1(JI),INUMLEV2(JI)
1531     CALL SEARCH_FIELD(IPAR,ITYP,ILEV1,ILEV2,IGRIB,INUM)
1532     IF (INUM < 0) THEN
1533       WRITE (ILUOUT0,'(A,I3,A,I2,A)') ' -> 2D field ',INUMGRIB(JI),' is missing - abort'
1534       !callabortstop
1535       CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
1536       CALL ABORT
1537       STOP
1538     END IF
1539     CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
1540     ALLOCATE(IINLO(INJ))
1541     CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM),IIU,IJU,ZLONOUT,ZLATOUT,&
1542         ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
1543     CALL GRIB_GET_SIZE(IGRIB(INUM),'values',ISIZE)
1544     ALLOCATE(ZVALUE(ISIZE))
1545     CALL GRIB_GET(IGRIB(INUM_ZS),'values',ZVALUE)
1546     ALLOCATE(ZOUT(INO))
1547     CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
1548          ZVALUE,INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE. )
1549     DEALLOCATE(IINLO)
1550     DEALLOCATE(ZVALUE)
1551     CALL ARRAY_1D_TO_2D(INO,ZOUT,IIU,IJU,XDUMMY_2D(:,:,JI))
1552     DEALLOCATE (ZOUT)
1553   END DO
1554 !
1555 ENDIF
1556 !
1557 !---------------------------------------------------------------------------------------
1558 !
1559 !* 3. VERTICAL GRID
1560 !
1561 IF (HFILE(1:3)=='ATM') THEN
1562   WRITE (ILUOUT0,'(A)') ' | Reading of vertical grid in progress'
1563   CALL READ_VER_GRID(HPRE_REAL1)
1564 END IF
1565
1566 !
1567 !---------------------------------------------------------------------------------------
1568 !
1569 !* 4. Free all temporary allocations
1570 !
1571 DEALLOCATE (ZLATOUT)
1572 DEALLOCATE (ZLONOUT)
1573 DEALLOCATE (ZXOUT)
1574 DEALLOCATE (ZYOUT)
1575 DEALLOCATE(ZPARAM)
1576 DEALLOCATE(ZPARAM_ZS)
1577 DEALLOCATE(IINLO_ZS)
1578 DO JLOOP=1,ICOUNT
1579   CALL GRIB_RELEASE(IGRIB(JLOOP))
1580 ENDDO
1581 DEALLOCATE(IGRIB)
1582
1583 WRITE (ILUOUT0,'(A,A4,A)') ' -- Grib decoder for ',HFILE,' file ended successfully'
1584 !
1585 !---------------------------------------------------------------------------------------
1586 !---------------------------------------------------------------------------------------
1587 !---------------------------------------------------------------------------------------
1588 !
1589 !
1590
1591 !
1592 CONTAINS
1593 !
1594 !
1595 !     ##########################################################################
1596       SUBROUTINE ARRAY_1D_TO_2D (KN1,P1,KL1,KL2,P2)
1597 !     ##########################################################################
1598 !
1599 !       Small routine used to store a linear array into a 2 dimension array
1600 !
1601 IMPLICIT NONE
1602 INTEGER,                INTENT(IN)  :: KN1
1603 REAL,DIMENSION(KN1),    INTENT(IN)  :: P1
1604 INTEGER,                INTENT(IN)  :: KL1
1605 INTEGER,                INTENT(IN)  :: KL2
1606 REAL,DIMENSION(KL1,KL2),INTENT(OUT) :: P2
1607 INTEGER                 :: JLOOP1_A1T2
1608 INTEGER                 :: JLOOP2_A1T2
1609 INTEGER                 :: JPOS_A1T2
1610 !
1611 IF (KN1 < KL1*KL2) THEN
1612   WRITE (ILUOUT0,'(A)') ' | Error in "ARRAY_1D_TO_2D", sizes do not match - abort'
1613  !callabortstop
1614   CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
1615   CALL ABORT
1616   STOP
1617 END IF
1618 JPOS_A1T2 = 1
1619 DO JLOOP2_A1T2 = 1, KL2
1620   DO JLOOP1_A1T2 = 1, KL1
1621     P2(JLOOP1_A1T2,JLOOP2_A1T2) = P1(JPOS_A1T2)
1622     JPOS_A1T2 = JPOS_A1T2 + 1
1623   END DO
1624 END DO
1625 END SUBROUTINE ARRAY_1D_TO_2D
1626 !
1627 !
1628 !---------------------------------------------------------------------------------------
1629 !---------------------------------------------------------------------------------------
1630 !---------------------------------------------------------------------------------------
1631 !#################################################################################
1632 SUBROUTINE SEARCH_FIELD(KPARAM,KLTYPE,KLEV1,KLEV2,KGRIB,KNUM)
1633 !#################################################################################
1634 ! search the grib message corresponding to KPARAM,KLTYPE,KLEV1,KLEV2 in all 
1635 ! the KGIRB messages
1636 !
1637 USE MODD_LUNIT
1638 USE GRIB_API
1639 !
1640 USE MODE_FM
1641 USE MODE_IO_ll
1642 !
1643 IMPLICIT NONE
1644 !
1645 !
1646 INTEGER,INTENT(IN)              :: KPARAM ! Parameter to read
1647 INTEGER,INTENT(IN)              :: KLTYPE ! Level type
1648 INTEGER,INTENT(IN)              :: KLEV1  ! Level parameter 1
1649 INTEGER,INTENT(IN)              :: KLEV2  ! Level parameter 2
1650 INTEGER(KIND=kindOfInt),DIMENSION(:),INTENT(IN) :: KGRIB ! number of grib messages
1651 INTEGER,INTENT(OUT)             :: KNUM  ! number of the message researched
1652 !
1653 ! Declaration of local variables
1654 !
1655 INTEGER :: IFOUND  ! Number of correct parameters
1656 INTEGER :: IRET    ! error code 
1657 INTEGER :: IPARAM  ! Parameter read
1658 INTEGER :: ILEV1   ! Level parameter 1
1659 INTEGER :: ILEV2   ! Level parameter 2
1660 INTEGER :: JLOOP   ! Dummy counter
1661 INTEGER :: IVERSION
1662 CHARACTER(LEN=20) :: YLTYPELU
1663 CHARACTER(LEN=20) :: YLTYPE
1664 !
1665 ! Variables used to display messages
1666 INTEGER :: ILUOUT0   ! Logical unit number of the listing
1667 !
1668 CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRET)
1669 !
1670 SELECT CASE (KLTYPE) 
1671 CASE(109)
1672   YLTYPE='hybrid'
1673 CASE(1)
1674   YLTYPE='surface'
1675 CASE DEFAULT
1676   YLTYPE='unknown'
1677 END SELECT
1678 DO JLOOP=1,SIZE(KGRIB)
1679       IFOUND = 0
1680       ! 
1681       CALL GRIB_GET(KGRIB(JLOOP),'editionNumber',IVERSION,IRET_GRIB)
1682       IF (IRET_GRIB >   0) THEN
1683         WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
1684         CYCLE
1685       ELSE IF (IRET_GRIB == -6) THEN
1686         WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
1687         CYCLE
1688       ENDIF
1689       IF (IVERSION == 2) THEN
1690         CALL GRIB_GET(KGRIB(JLOOP),'paramId',IPARAM,IRET_GRIB)
1691       ELSE
1692         CALL GRIB_GET(KGRIB(JLOOP),'indicatorOfParameter',IPARAM,IRET_GRIB)
1693       ENDIF
1694       IF (IRET_GRIB >   0) THEN
1695         WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
1696         CYCLE
1697       ELSE IF (IRET_GRIB == -6) THEN
1698         WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
1699         CYCLE
1700       ENDIF
1701       !
1702       CALL GRIB_GET(KGRIB(JLOOP),'typeOfLevel',YLTYPELU,IRET_GRIB)
1703       IF (IRET_GRIB >   0) THEN
1704         WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
1705         CYCLE
1706       ELSE IF (IRET_GRIB == -6) THEN
1707         WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
1708         CYCLE
1709       ENDIF
1710       !
1711       CALL GRIB_GET(KGRIB(JLOOP),'topLevel',ILEV1,IRET_GRIB)
1712       IF (IRET_GRIB >   0) THEN
1713         WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
1714         CYCLE
1715       ELSE IF (IRET_GRIB == -6) THEN
1716         WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
1717         CYCLE
1718       ENDIF
1719       ! 
1720       CALL GRIB_GET(KGRIB(JLOOP),'bottomLevel',ILEV2,IRET_GRIB)
1721       IF (IRET_GRIB >   0) THEN
1722         WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
1723         CYCLE
1724       ELSE IF (IRET_GRIB == -6) THEN
1725         WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
1726         CYCLE
1727       ENDIF
1728       IF (IPARAM==KPARAM) IFOUND = 1
1729       IF (YLTYPELU==YLTYPE .OR. -1==KLTYPE) IFOUND = IFOUND + 1
1730       IF (ILEV1== KLEV1  .OR. -1== KLEV1) IFOUND = IFOUND + 1
1731       IF (ILEV2== KLEV2 .OR. -1== KLEV2) IFOUND = IFOUND + 1
1732       IF (IFOUND == 4) THEN
1733           KNUM=JLOOP
1734           EXIT
1735       ELSE  ! field not found
1736           KNUM=-1
1737       END IF
1738  
1739 END DO
1740
1741 END SUBROUTINE SEARCH_FIELD
1742
1743 !#################################################################################
1744 SUBROUTINE COORDINATE_CONVERSION(KMODEL,KGRIB,KNOLON,KNOLARG,&
1745                            PLONOUT,PLATOUT,PLXOUT,PLYOUT,KNI,PPARAM,KINLO)
1746 !#################################################################################
1747 !perform coordinate conversion from lat/lon system to x,y (depends on the grib
1748 ! type)
1749 !!   AUTHOR
1750 !!   ------
1751 !!
1752 !!   G. Tanguy
1753 !!
1754 !!   MODIFICATIONS
1755 !!   -------------
1756 !!
1757 !!   Original       08/06/2010
1758
1759 USE MODD_CST
1760 USE MODI_LATLONTOXY
1761 USE GRIB_API
1762 !
1763 IMPLICIT NONE
1764 !
1765 !
1766 INTEGER(KIND=kindOfInt),INTENT(IN)                            :: KGRIB  ! number of the grib message
1767 INTEGER,INTENT(IN)                            :: KMODEL ! number of the model
1768 INTEGER,INTENT(OUT)                           :: KNI    ! number of points
1769 INTEGER,INTENT(IN)                            :: KNOLON,KNOLARG  ! Number of output points
1770 REAL,DIMENSION( KNOLON*KNOLARG),INTENT(IN)    :: PLONOUT ! Output coordinate,
1771 REAL,DIMENSION( KNOLON*KNOLARG),INTENT(IN)    :: PLATOUT ! lat/lon system
1772 REAL,DIMENSION( KNOLON*KNOLARG),INTENT(INOUT) :: PLXOUT  ! Converted output coordinates
1773 REAL,DIMENSION( KNOLON*KNOLARG),INTENT(INOUT) :: PLYOUT  ! (depends on Grib Grid type)
1774 REAL,DIMENSION(:),INTENT(INOUT)               :: PPARAM  ! output parameters of
1775 !                                                the grid to avoid many calculations
1776 INTEGER,DIMENSION(:),INTENT(INOUT)            :: KINLO   ! Number of points along a parallel
1777 !===============================
1778 INTEGER                                       :: IINLA ! Number of points along a meridian
1779 INTEGER                                       :: JLOOP1,JLOOP2 ! Dummy counter
1780 INTEGER                                       :: INO ! Number of output points
1781 REAL                                          :: ZILA1  ! Grib first point latitude
1782 REAL                                          :: ZILO1  ! Grib first point longitude
1783 REAL                                          :: ZILA2  ! Grib last point latitude
1784 REAL                                          :: ZILO2  ! Grib last point longitude
1785 REAL                                          :: ZILASP ! Grib streching pole lat
1786 REAL                                          :: ZILOSP ! Grib streching pole lon
1787 LOGICAL                                       :: GREADY ! Used to test if projection is needed
1788 INTEGER                                       :: ILENX ! nb points in X
1789 INTEGER                                       :: ILENY ! nb points in Y
1790 INTEGER                                       :: IEARTH  ! 
1791 REAL                                          :: ZSTRECH ! streching of arpege grid
1792 CHARACTER(LEN=20)                             :: CGRID   ! type of the grid
1793 INTEGER(KIND=kindOfInt)                                       :: IMISSING ! dummy variable
1794 ! Aladin projection
1795 REAL                                          :: ZALALAT0  ! Grid definition parameters
1796 REAL                                          :: ZALALON0  !  |
1797 REAL                                          :: ZALALATOR !  |
1798 REAL                                          :: ZALALONOR !  |
1799 REAL                                          :: ZALARPK   !  |
1800 REAL, DIMENSION(:,:), ALLOCATABLE             :: ZXM       ! Intermediate arrays
1801 REAL, DIMENSION(:,:), ALLOCATABLE             :: ZYM       !  |
1802 REAL, DIMENSION(:,:), ALLOCATABLE             :: ZLONM     !  |
1803 REAL, DIMENSION(:,:), ALLOCATABLE             :: ZLATM     !  |
1804 ! CEP projection
1805 REAL, DIMENSION(:), ALLOCATABLE               :: ZLATGRIB
1806 REAL, DIMENSION(:), ALLOCATABLE               :: ZLONGRIB
1807 INTEGER                                       :: INBLATGRIB,INBLONGRIB
1808 !JUAN
1809 INTEGER(KIND=kindOfInt),DIMENSION(:),ALLOCATABLE            :: INLO_GRIB   ! Number of points along a parallel
1810 !JUAN
1811 !
1812 !--------------------------------------------------------------------------------
1813 !
1814 !JUAN
1815 ALLOCATE(INLO_GRIB(SIZE(KINLO)))
1816 !JUAN
1817 INO= KNOLON*KNOLARG
1818 CALL GRIB_GET(KGRIB,'typeOfGrid',CGRID)
1819 SELECT CASE (KMODEL)
1820 !
1821 CASE(0,5) ! CEP/MOCAGE
1822 ! en theorie il faut ces 4 lignes
1823 !  CALL GRIB_GET(KGRIB,'latitudeOfFirstGridPointInDegrees',ZILA1)
1824 !  CALL GRIB_GET(KGRIB,'longitudeOfFirstGridPointInDegrees',ZILO1)
1825 !  CALL GRIB_GET(KGRIB,'latitudeOfLastGridPointInDegrees',ZILA2)
1826 !  CALL GRIB_GET(KGRIB,'longitudeOfLastGridPointInDegrees',ZILO2)
1827 ! pourtant au passage de GRIB1 a GRIB2 les arrondi etait fait differement
1828 ! et on n'obtenais pas les meme resultat entre un fichier grib1 et le meme 
1829 ! convertit en GRIB2
1830 ! Du coup en faisant ce qui suit on prend une valeur recalculee par grib_api
1831 ! suivant l'ordre N de la gausienne donc plus precise et donc la meme entre le
1832 ! GRIB1 et le GRIB2
1833   CALL GRIB_GET(KGRIB,'Nj',IINLA,IRET_GRIB)
1834   CALL GRIB_GET_SIZE(KGRIB,'latitudes',INBLATGRIB)
1835   CALL GRIB_GET_SIZE(KGRIB,'longitudes',INBLONGRIB)
1836   ALLOCATE(ZLATGRIB(INBLATGRIB))
1837   ALLOCATE(ZLONGRIB(INBLONGRIB))
1838   CALL GRIB_GET(KGRIB,'latitudes',ZLATGRIB)
1839   CALL GRIB_GET(KGRIB,'longitudes',ZLONGRIB)
1840   ZILA1=MAXVAL(ZLATGRIB)
1841   ZILO1=MINVAL(ZLONGRIB)
1842   ZILA2=MINVAL(ZLATGRIB)
1843   ZILO2=MAXVAL(ZLONGRIB)
1844   KNI=0
1845   CALL GRIB_IS_MISSING(KGRIB,'pl',IMISSING,IRET_GRIB)
1846   IF (IRET_GRIB /= 0 .OR. IMISSING==1)  THEN   ! pl not present
1847     CALL GRIB_GET(KGRIB,'Ni',INLO_GRIB(1),IRET_GRIB)
1848     INLO_GRIB(2:)=INLO_GRIB(1)
1849     KNI=IINLA*INLO_GRIB(1)
1850     GREADY= (PPARAM(1)==INLO_GRIB(1) .AND. PPARAM(2)==IINLA .AND.&
1851       PPARAM(3)==ZILA1 .AND. PPARAM(4)==ZILO1 .AND.&
1852       PPARAM(5)==ZILA2 .AND. PPARAM(6)==ZILO2)
1853    PPARAM(1)=INLO_GRIB(1)
1854    PPARAM(2)=IINLA
1855    PPARAM(3)=ZILA1
1856    PPARAM(4)=ZILO1 
1857    PPARAM(5)=ZILA2
1858    PPARAM(6)=ZILO2
1859   ELSE ! pl present in the grib
1860     CALL GRIB_GET(KGRIB,'pl',INLO_GRIB)
1861     DO JLOOP1=1 ,IINLA
1862       KNI = KNI + INLO_GRIB(JLOOP1)
1863     ENDDO
1864     GREADY= (PPARAM(1)==INLO_GRIB(1) .AND.&
1865       PPARAM(3)==ZILA1 .AND. PPARAM(4)==ZILO1 .AND.&
1866       PPARAM(5)==ZILA2 .AND. PPARAM(6)==ZILO2)
1867     PPARAM(1)=INLO_GRIB(1)
1868     PPARAM(2)=IINLA
1869     PPARAM(3)=ZILA1
1870     PPARAM(4)=ZILO1 
1871     PPARAM(5)=ZILA2
1872     PPARAM(6)=ZILO2    
1873   END IF
1874   IF (.NOT. GREADY) THEN
1875     PLXOUT=PLONOUT
1876     PLYOUT=PLATOUT
1877   ENDIF
1878 !
1879 CASE(1) ! ALADIN
1880 !
1881   CALL GRIB_GET(KGRIB,'Nj',IINLA,IRET_GRIB)
1882   CALL GRIB_GET(KGRIB,'Ni',INLO_GRIB(1),IRET_GRIB)
1883   INLO_GRIB(2:)=INLO_GRIB(1)
1884   CALL GRIB_GET(KGRIB,'DxInMetres',ILENX)
1885   CALL GRIB_GET(KGRIB,'DyInMetres',ILENY)
1886   CALL GRIB_GET(KGRIB,'LoVInDegrees',ZALALON0)
1887   CALL GRIB_GET(KGRIB,'Latin1InDegrees',ZALALAT0)
1888   KNI = IINLA*INLO_GRIB(1)
1889   ZILA1 = 0.
1890   ZILO1 = 0.
1891   ZILA2 = ZILA1 + (IINLA   -1)*ILENY
1892   ZILO2 = ZILO1 + (INLO_GRIB(1)-1)*ILENX
1893   GREADY= (PPARAM(1)==INLO_GRIB(1) .AND. PPARAM(2)==IINLA .AND.&
1894       PPARAM(3)==ZILA1 .AND. PPARAM(4)==ZILO1 .AND.&
1895       PPARAM(5)==ZILA2 .AND. PPARAM(6)==ZILO2.AND.&
1896       PPARAM(7)==ILENX .AND. PPARAM(8)==ILENY.AND.&
1897       PPARAM(9)==ZALALAT0 .AND. PPARAM(10)==ZALALON0)
1898   IF(.NOT. GREADY) THEN
1899      PPARAM(1)=INLO_GRIB(1)
1900      PPARAM(2)=IINLA
1901      PPARAM(3)=ZILA1
1902      PPARAM(4)=ZILO1 
1903      PPARAM(5)=ZILA2
1904      PPARAM(6)=ZILO2    
1905      PPARAM(7)=ILENX
1906      PPARAM(8)=ILENY 
1907      PPARAM(9)=ZALALAT0
1908      PPARAM(10)=ZALALON0
1909 !        
1910      IF (ZALALON0 > 180.) ZALALON0 = ZALALON0 - 360.
1911      CALL GRIB_GET(KGRIB,'latitudeOfFirstGridPointInDegrees',ZALALATOR)
1912     CALL GRIB_GET(KGRIB,'longitudeOfFirstGridPointInDegrees',ZALALONOR)
1913      IF (ZALALONOR > 180.) ZALALONOR = ZALALONOR - 360.
1914      ZALARPK = SIN(ZALALAT0/180.*XPI)
1915      ALLOCATE (ZXM(KNOLON,KNOLARG))
1916      ALLOCATE (ZYM(KNOLON,KNOLARG))
1917      ALLOCATE (ZLONM(KNOLON,KNOLARG))
1918      ALLOCATE (ZLATM(KNOLON,KNOLARG))
1919      JLOOP1=0
1920      DO JLOOP2=1, KNOLARG
1921        ZLONM(1:KNOLON,JLOOP2) = PLONOUT(1+JLOOP1:KNOLON+JLOOP1)
1922        ZLATM(1:KNOLON,JLOOP2) = PLATOUT(1+JLOOP1:KNOLON+JLOOP1)
1923        JLOOP1 = JLOOP1+KNOLON
1924      END DO
1925      CALL SM_LATLONTOXY_A (ZALALAT0,ZALALON0,ZALARPK,ZALALATOR,ZALALONOR, &
1926                            ZXM,ZYM,ZLATM,ZLONM,KNOLON,KNOLARG,6367470.)
1927      JLOOP1=0
1928      DO JLOOP2=1, KNOLARG
1929         PLXOUT(1+JLOOP1:KNOLON+JLOOP1)=ZXM(1:KNOLON,JLOOP2)
1930         PLYOUT(1+JLOOP1:KNOLON+JLOOP1)=ZYM(1:KNOLON,JLOOP2)
1931         JLOOP1 = JLOOP1+KNOLON
1932      ENDDO
1933      DEALLOCATE (ZLATM)
1934      DEALLOCATE (ZLONM)
1935      DEALLOCATE (ZYM)
1936      DEALLOCATE (ZXM)
1937   END IF
1938 !
1939 CASE(2) ! ALADIN REUNION
1940 !
1941   CALL GRIB_GET(KGRIB,'Nj',IINLA,IRET_GRIB)
1942   CALL GRIB_GET(KGRIB,'Ni',INLO_GRIB(1),IRET_GRIB)
1943   INLO_GRIB(2:)=INLO_GRIB(1)
1944   CALL GRIB_GET(KGRIB,'DiInMetres',ILENX)
1945   CALL GRIB_GET(KGRIB,'DjInMetres',ILENY)
1946   CALL GRIB_GET(KGRIB,'LaDInDegrees',ZALALAT0)
1947   KNI = IINLA*INLO_GRIB(1)
1948   ZILA1 = 0.
1949   ZILO1 = 0.
1950   ZILA2 = ZILA1 + (IINLA   -1)*ILENY
1951   ZILO2 = ZILO1 + (INLO_GRIB(1)-1)*ILENX
1952   GREADY= (PPARAM(1)==INLO_GRIB(1) .AND. PPARAM(2)==IINLA .AND.&
1953       PPARAM(3)==ZILA1 .AND. PPARAM(4)==ZILO1 .AND.&
1954       PPARAM(5)==ZILA2 .AND. PPARAM(6)==ZILO2.AND.&
1955       PPARAM(7)==ILENX .AND. PPARAM(8)==ILENY.AND.&
1956       PPARAM(9)==ZALALAT0)
1957   IF(.NOT. GREADY) THEN
1958     PPARAM(1)=INLO_GRIB(1)
1959     PPARAM(2)=IINLA
1960     PPARAM(3)=ZILA1
1961     PPARAM(4)=ZILO1 
1962     PPARAM(5)=ZILA2
1963     PPARAM(6)=ZILO2    
1964     PPARAM(7)=ILENX
1965     PPARAM(8)=ILENY 
1966     PPARAM(9)=ZALALAT0
1967     ZALALON0 = 0.
1968     CALL GRIB_GET(KGRIB,'latitudeOfFirstGridPointInDegrees',ZALALATOR)
1969     CALL GRIB_GET(KGRIB,'longitudeOfFirstGridPointInDegrees',ZALALONOR)
1970     IF (ZALALONOR > 180.) ZALALONOR = ZALALONOR - 360.
1971     ZALARPK = 0
1972     ALLOCATE (ZXM(KNOLON,KNOLARG))
1973     ALLOCATE (ZYM(KNOLON,KNOLARG))
1974     ALLOCATE (ZLONM(KNOLON,KNOLARG))
1975     ALLOCATE (ZLATM(KNOLON,KNOLARG))
1976     JLOOP1=0
1977     DO JLOOP2=1, KNOLARG
1978        ZLONM(1:KNOLON,JLOOP2) = PLONOUT(1+JLOOP1:KNOLON+JLOOP1)
1979        ZLATM(1:KNOLON,JLOOP2) = PLATOUT(1+JLOOP1:KNOLON+JLOOP1)
1980        JLOOP1 = JLOOP1+KNOLON
1981     END DO
1982     CALL GRIB_GET(KGRIB,'earthIsOblate',IEARTH)
1983     IF (IEARTH==0) THEN
1984         CALL SM_LATLONTOXY_A (ZALALAT0,ZALALON0,ZALARPK,ZALALATOR,ZALALONOR, &
1985                               ZXM,ZYM,ZLATM,ZLONM,KNOLON,KNOLARG,6367470.)
1986     ELSE
1987         CALL SM_LATLONTOXY_A (ZALALAT0,ZALALON0,ZALARPK,ZALALATOR,ZALALONOR, &
1988                               ZXM,ZYM,ZLATM,ZLONM,KNOLON,KNOLARG)
1989     END IF
1990     JLOOP1=0
1991     DO JLOOP2=1, KNOLARG
1992        PLXOUT(1+JLOOP1:KNOLON+JLOOP1)=ZXM(1:KNOLON,JLOOP2)
1993        PLYOUT(1+JLOOP1:KNOLON+JLOOP1)=ZYM(1:KNOLON,JLOOP2)
1994        JLOOP1 = JLOOP1+KNOLON
1995     ENDDO
1996     DEALLOCATE (ZLATM)
1997     DEALLOCATE (ZLONM)
1998     DEALLOCATE (ZYM)
1999     DEALLOCATE (ZXM)       
2000   END IF
2001 !
2002 CASE(3,4) ! ARPEGE
2003 !
2004   CALL GRIB_GET(KGRIB,'latitudeOfFirstGridPointInDegrees',ZILA1)
2005   CALL GRIB_GET(KGRIB,'longitudeOfFirstGridPointInDegrees',ZILO1)
2006   CALL GRIB_GET(KGRIB,'latitudeOfLastGridPointInDegrees',ZILA2)
2007   CALL GRIB_GET(KGRIB,'longitudeOfLastGridPointInDegrees',ZILO2)
2008   CALL GRIB_GET(KGRIB,'latitudeOfStretchingPoleInDegrees',ZILASP)
2009   CALL GRIB_GET(KGRIB,'longitudeOfStretchingPoleInDegrees',ZILOSP)
2010   CALL GRIB_GET(KGRIB,'stretchingFactor',ZSTRECH)
2011   CALL GRIB_GET(KGRIB,'Nj',IINLA,IRET_GRIB)
2012 !
2013   KNI=0
2014   CALL GRIB_IS_MISSING(KGRIB,'pl',IRET_GRIB)
2015   IF (IRET_GRIB == 1)  THEN !  regular
2016      CALL GRIB_GET(KGRIB,'Ni',INLO_GRIB(1),IRET_GRIB)
2017      INLO_GRIB(2:)=INLO_GRIB(1)
2018      KNI=IINLA*INLO_GRIB(1)
2019      GREADY= (PPARAM(1)==INLO_GRIB(1) .AND. PPARAM(2)==IINLA .AND.&
2020        PPARAM(3)==ZILA1 .AND. PPARAM(4)==ZILO1 .AND.&
2021        PPARAM(5)==ZILA2 .AND. PPARAM(6)==ZILO2 .AND.&
2022        PPARAM(7)==ZILASP .AND. PPARAM(8)==ZILOSP .AND.&
2023        PPARAM(9)==ZSTRECH)
2024   ELSE !  quasi-regular
2025      CALL GRIB_GET(KGRIB,'pl',INLO_GRIB)
2026      DO JLOOP1=1 ,IINLA
2027         KNI = KNI + INLO_GRIB(JLOOP1)
2028      ENDDO
2029      GREADY= (PPARAM(1)==INLO_GRIB(1) .AND.&
2030        PPARAM(3)==ZILA1 .AND. PPARAM(4)==ZILO1 .AND.&
2031        PPARAM(5)==ZILA2 .AND. PPARAM(6)==ZILO2 .AND.&
2032        PPARAM(7)==ZILASP .AND. PPARAM(8)==ZILOSP .AND.&
2033        PPARAM(9)==ZSTRECH)
2034   END IF
2035 !
2036   IF (.NOT. GREADY) THEN
2037     CALL ARPEGE_STRETCH_A(INO,ZILASP,ZILOSP, &
2038         ZSTRECH,PLATOUT,PLONOUT,PLYOUT,PLXOUT)
2039      PPARAM(1)=INLO_GRIB(1)
2040      PPARAM(2)=IINLA
2041      PPARAM(3)=ZILA1
2042      PPARAM(4)=ZILO1 
2043      PPARAM(5)=ZILA2
2044      PPARAM(6)=ZILO2
2045      PPARAM(7)=ZILASP
2046      PPARAM(8)=ZILOSP
2047      PPARAM(9)=ZSTRECH    
2048   END IF
2049 END SELECT
2050 !JUAN
2051 KINLO=INLO_GRIB
2052 !JUAN
2053 END SUBROUTINE COORDINATE_CONVERSION
2054 !
2055 !     ###################################################################
2056       SUBROUTINE ARPEGE_STRETCH_A(KN,PLAP,PLOP,PCOEF,PLAR,PLOR,PLAC,PLOC)
2057 !     ###################################################################
2058 !!****  *ARPEGE_STRETCH_A* - Projection to Arpege stretched grid
2059 !!
2060 !!   PURPOSE
2061 !!   -------
2062 !!
2063 !!   Projection from standard Lat,Lon grid to Arpege stretched grid
2064 !!
2065 !!   METHOD
2066 !!   ------
2067 !!
2068 !!   The projection is defined in two steps :
2069 !!    1. A rotation to place the stretching pole at the north pole
2070 !!    2. The stretching
2071 !!   This routine is a basic implementation of the informations founded in
2072 !!     'Note de travail Arpege nr.3'
2073 !!     'Transformation de coordonnees'
2074 !!     J.F.Geleyn 1988
2075 !!   This document describes a slightly different transformation in 3 steps. Only the
2076 !!   two first steps are to be taken in account (at the time of writing this paper has
2077 !!   not been updated).
2078 !!
2079 !!   EXTERNAL
2080 !!   --------
2081 !!
2082 !!   Module MODD_CST
2083 !!     XPI
2084 !!
2085 !!   IMPLICIT ARGUMENTS
2086 !!   ------------------
2087 !!
2088 !!   REFERENCE
2089 !!   ---------
2090 !!
2091 !!   This routine is based on :
2092 !!     'Note de travail ARPEGE' number 3
2093 !!     by J.F. GELEYN (may 1988)
2094 !!
2095 !!   AUTHOR
2096 !!   ------
2097 !!
2098 !!   V.Bousquet
2099 !!
2100 !!   MODIFICATIONS
2101 !!   -------------
2102 !!
2103 !!   Original       07/01/1999
2104 !!
2105 !----------------------------------------------------------------------------
2106 !
2107 !*   0. DECLARATIONS
2108 !    ---------------
2109 !
2110 USE MODD_CST
2111 !
2112 IMPLICIT NONE
2113 !
2114 !*   0.1. Declaration of arguments
2115 !    -----------------------------
2116 !
2117 INTEGER,             INTENT(IN)  :: KN            ! Number of points to convert
2118 REAL,                INTENT(IN)  :: PLAP          ! Latitude of stretching pole
2119 REAL,                INTENT(IN)  :: PLOP          ! Longitude of stretching pole
2120 REAL,                INTENT(IN)  :: PCOEF         ! Stretching coefficient
2121 REAL, DIMENSION(KN), INTENT(IN)  :: PLAR          ! Lat. of points
2122 REAL, DIMENSION(KN), INTENT(IN)  :: PLOR          ! Lon. of points
2123 REAL, DIMENSION(KN), INTENT(OUT) :: PLAC          ! Computed pseudo-lat. of points
2124 REAL, DIMENSION(KN), INTENT(OUT) :: PLOC          ! Computed pseudo-lon. of points
2125 !
2126 !*   0.2. Declaration of local variables
2127 !    -----------------------------------
2128 !
2129 REAL                       :: ZSINSTRETCHLA ! Sine of stretching point lat.
2130 REAL                       :: ZSINSTRETCHLO ! Sine of stretching point lon.
2131 REAL                       :: ZCOSSTRETCHLA ! Cosine of stretching point lat.
2132 REAL                       :: ZCOSSTRETCHLO ! Cosine of stretching point lon.
2133 REAL                       :: ZSINLA        ! Sine of computed point latitude
2134 REAL                       :: ZSINLO        ! Sine of computed point longitude
2135 REAL                       :: ZCOSLA        ! Cosine of computed point latitude
2136 REAL                       :: ZCOSLO        ! Cosine of computed point longitude
2137 REAL                       :: ZSINLAS       ! Sine of point's pseudo-latitude
2138 REAL                       :: ZSINLOS       ! Sine of point's pseudo-longitude
2139 REAL                       :: ZCOSLOS       ! Cosine of point's pseudo-lon.
2140 REAL                       :: ZA,ZB,ZD      ! Dummy variables used for
2141 REAL                       :: ZX,ZY         ! computations
2142 !
2143 INTEGER                    :: JLOOP1        ! Dummy loop counter
2144 !
2145 !----------------------------------------------------------------------------
2146 !
2147 ZSINSTRETCHLA = SIN(PLAP*XPI/180.)
2148 ZCOSSTRETCHLA = COS(PLAP*XPI/180.)
2149 ZSINSTRETCHLO = SIN(PLOP*XPI/180.)
2150 ZCOSSTRETCHLO = COS(PLOP*XPI/180.)
2151 ! L = longitude (0 = Greenwich, + toward east)
2152 ! l = latitude (90 = N.P., -90 = S.P.)
2153 ! p stands for stretching pole
2154 PLAC(:) = PLAR(:) * XPI / 180.
2155 PLOC(:) = PLOR(:) * XPI / 180.
2156 ! A = 1 + c.c
2157 ZA = 1. + PCOEF*PCOEF
2158 ! B = 1 - c.c
2159 ZB = 1. - PCOEF*PCOEF
2160 DO JLOOP1=1, KN
2161   ZSINLA = SIN(PLAC(JLOOP1))
2162   ZCOSLA = COS(PLAC(JLOOP1))
2163   ZSINLO = SIN(PLOC(JLOOP1))
2164   ZCOSLO = COS(PLOC(JLOOP1))
2165   ! X = cos(Lp-L)
2166   ZX = ZCOSLO*ZCOSSTRETCHLO + ZSINLO*ZSINSTRETCHLO
2167   ! Y = sin(Lp-L)
2168   ZY = ZSINSTRETCHLO*ZCOSLO - ZSINLO*ZCOSSTRETCHLO
2169   ! D = (1+c.c) + (1-c.c)(sin lp.sin l + cos lp.cos l.cos(Lp-L))
2170   ZD = ZA + ZB*(ZSINSTRETCHLA*ZSINLA+ZCOSSTRETCHLA*ZCOSLA*ZX)
2171   !          (1-c.c)+(1+c.c)((sin lp.sin l + cos lp.cos l.cos(Lp-L))
2172   ! sin lr = -------------------------------------------------------
2173   !                                  D
2174   ZSINLAS = (ZB + ZA*(ZSINSTRETCHLA*ZSINLA+ZCOSSTRETCHLA*ZCOSLA*ZX)) / ZD
2175   ! D' = D * cos lr
2176   ZD = ZD * (AMAX1(1e-6,SQRT(1.-ZSINLAS*ZSINLAS)))
2177   !          2.c.(cos lp.sin l - sin lp.cos l.cos(Lp-L))
2178   ! cos Lr = -------------------------------------------
2179   !                              D'
2180   ZCOSLOS = 2.*PCOEF*(ZCOSSTRETCHLA*ZSINLA-ZSINSTRETCHLA*ZCOSLA*ZX) / ZD
2181   !          2.c.cos l.cos(Lp-L)
2182   ! sin Lr = -------------------
2183   !                  D'
2184   ZSINLOS = 2.*PCOEF*(ZCOSLA*ZY) / ZD
2185   ! saturations (corrects calculation errors)
2186   ZSINLAS = MAX(ZSINLAS,-1.)
2187   ZSINLAS = MIN(ZSINLAS, 1.)
2188   ZCOSLOS = MAX(ZCOSLOS,-1.)
2189   ZCOSLOS = MIN(ZCOSLOS, 1.)
2190   ! back from sine & cosine
2191   PLAC(JLOOP1) = ASIN(ZSINLAS)
2192   IF (ZSINLOS>0) THEN
2193     PLOC(JLOOP1) =  ACOS(ZCOSLOS)
2194   ELSE
2195     PLOC(JLOOP1) = -ACOS(ZCOSLOS)
2196   ENDIF
2197 ENDDO
2198 PLOC(:) = PLOC(:) * 180. / XPI
2199 PLAC(:) = PLAC(:) * 180. / XPI
2200 RETURN
2201 END SUBROUTINE ARPEGE_STRETCH_A
2202 !
2203 !
2204 END SUBROUTINE READ_ALL_DATA_GRIB_CASE