SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
prep_grib_grid.F90
Go to the documentation of this file.
1 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SFX_LIC for details. version 1.
5 ! #########
6  SUBROUTINE prep_grib_grid(HGRIB,KLUOUT,HINMODEL,HGRIDTYPE,TPTIME_GRIB)
7 ! ##########################################################################
8 !
9 !!**** *PREP_GRIB_GRID* - reads GRIB grid.
10 !!
11 !! PURPOSE
12 !! -------
13 !!
14 !!** METHOD
15 !! ------
16 !!
17 !! EXTERNAL
18 !! --------
19 !!
20 !! IMPLICIT ARGUMENTS
21 !! ------------------
22 !!
23 !!
24 !! REFERENCE
25 !! ---------
26 !!
27 !!
28 !! AUTHOR
29 !! ------
30 !!
31 !! V. Masson (from read_all_data_grib_case)
32 !!
33 !! MODIFICATIONS
34 !! -------------
35 !! Original 06/2003
36 !! S. Faroux 01/2011 : to use library GRIB_API instead of GRIBEX (from
37 !! read_all_data_grib_case)
38 !-------------------------------------------------------------------------------
39 !
40 !* 0. DECLARATIONS
41 ! ------------
42 !
43 USE modd_type_date_surf
44 !
45 USE mode_read_grib
46 !
47 USE modd_grid_rotlatlon
48 USE modd_grid_gauss, ONLY : xila1, xilo1, xila2, xilo2, ninla, ninlo, nilen, lrotpole, xcoef, xlap, xlop
49 USE modd_grid_arome, ONLY : xx, xy, nx, ny, xlat0, xlon0, xlator, xlonor, xrpk, xbeta
50 USE modd_grid_grib, ONLY : nni
51 USE modd_surf_par, ONLY : xundef, nundef
52 USE modd_csts, ONLY : xpi
53 !
54 USE yomhook ,ONLY : lhook, dr_hook
55 USE parkind1 ,ONLY : jprb
56 !
57 USE modi_abor1_sfx
58 !
59 IMPLICIT NONE
60 !
61 !* 0.1. Declaration of arguments
62 ! ------------------------
63 !
64  CHARACTER(LEN=*), INTENT(IN) :: hgrib ! Grib file name
65 INTEGER, INTENT(IN) :: kluout ! logical unit of output listing
66  CHARACTER(LEN=6), INTENT(OUT) :: hinmodel ! Grib originating model
67  CHARACTER(LEN=10), INTENT(OUT) :: hgridtype ! Grid type
68 TYPE (date_time) :: tptime_grib ! current date and time
69 
70 !
71 !* 0.2 Declaration of local variables
72 ! ------------------------------
73 ! General purpose variables
74 INTEGER(KIND=kindOfInt) :: iret ! Return code from subroutines
75 !
76 ! Variable involved in the task of reading the grib file
77 INTEGER(KIND=kindOfInt) :: imissing
78 INTEGER(KIND=kindOfInt) :: iunit
79 INTEGER(KIND=kindOfInt) :: igrib
80 INTEGER :: icenter ! number of center
81  CHARACTER(LEN=20) :: hgrid ! type of grid
82 INTEGER :: iscan, jscan
83 INTEGER :: ilenx ! nb points in X
84 INTEGER :: ileny ! nb points in Y
85 INTEGER :: itime
86 INTEGER :: iunittime,ip1
87 !
88 ! Grib Grid definition variables
89 INTEGER :: jloop1 ! Dummy counter
90 !JUAN
91 INTEGER(KIND=kindOfInt), DIMENSION(:), ALLOCATABLE :: ininlo_grib
92 !JUAN
93 REAL(KIND=JPRB) :: zhook_handle
94 !
95 !---------------------------------------------------------------------------------------
96 !
97 IF (lhook) CALL dr_hook('PREP_GRIB_GRID',0,zhook_handle)
98 WRITE (kluout,'(A)') ' -- Grib reader started'
99 !
100 ! open grib file
101  CALL grib_open_file(iunit,hgrib,'R',iret)
102 IF (iret /= 0) THEN
103  CALL abor1_sfx('PREP_GRIB_GRID: Error opening the grib file '//hgrib)
104 END IF
105 !
106  CALL grib_new_from_file(iunit,igrib,iret)
107 IF (iret /= 0) THEN
108  CALL abor1_sfx('PREP_GRIB_GRID: Error in reading the grib file')
109 END IF
110 !
111 ! close the grib file
112  CALL grib_close_file(iunit)
113 !
114 !---------------------------------------------------------------------------------------
115 !* 2. Fix originating center
116 !---------------------------------------------------------------------------------------
117 !
118  CALL grib_get(igrib,'centre',icenter,iret)
119 IF (iret /= 0) THEN
120  CALL abor1_sfx('PREP_GRIB_GRID: Error in reading center')
121 END IF
122 !
123  CALL grib_get(igrib,'typeOfGrid',hgrid,iret)
124 IF (iret /= 0) THEN
125  CALL abor1_sfx('PREP_GRIB_GRID: Error in reading type of grid')
126 END IF
127 !
128 SELECT CASE (icenter)
129 
130  CASE (96)
131  WRITE (kluout,'(A)') ' | Grib file from HARMONY'
132  hinmodel='ALADIN'
133  hgridtype='AROME '
134 
135  CASE (82)
136  WRITE (kluout,'(A)') ' | Grib file from HIRLAM'
137  hinmodel='HIRLAM'
138  hgridtype='ROTLATLON '
139 
140  CASE (98)
141  WRITE (kluout,'(A)') ' | Grib file from European Center for Medium-range Weather Forecast'
142  hinmodel = 'ECMWF '
143  hgridtype= 'GAUSS '
144 
145  CASE (85)
146  SELECT CASE (hgrid)
147 
148  CASE('regular_gg')
149  WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Arpege model'
150  WRITE (kluout,'(A)') 'but same grid as ECMWF model (unstretched)'
151  hinmodel = 'ARPEGE'
152  hgridtype= 'GAUSS '
153 
154  CASE('reduced_gg')
155  WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Arpege model'
156  WRITE (kluout,'(A)') 'but reduced grid'
157  hinmodel = 'ARPEGE'
158  hgridtype= 'GAUSS '
159 
160  CASE('regular_ll')
161  WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Mocage model'
162  hinmodel = 'MOCAGE'
163  hgridtype= 'LATLON '
164 
165  CASE('unknown_PLPresent')
166  WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Arpege model'
167  hinmodel = 'ARPEGE'
168  hgridtype= 'ROTGAUSS '
169 
170  CASE('lambert')
171  WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Aladin france model'
172  hinmodel = 'ALADIN'
173  hgridtype= 'AROME '
174 
175  CASE('mercator')
176  WRITE (kluout,'(A)') ' | Grib file from French Weather Service - Aladin reunion model'
177  hinmodel = 'ALADIN'
178  hgridtype= 'MERCATOR '
179 
180  END SELECT
181 
182  CASE default
183  CALL abor1_sfx('PREP_GRIB_GRID: GRIB FILE FORMAT NOT SUPPORTED')
184 
185 END SELECT
186 
187 !---------------------------------------------------------------------------------------
188 !* 3. Number of points
189 !---------------------------------------------------------------------------------------
190 !
191 nx = nundef
192 ny = nundef
193 ninla = nundef
194 nilen = nundef
195 IF (ALLOCATED(ninlo)) DEALLOCATE(ninlo)
196 !
197 SELECT CASE (hgridtype)
198 
199  CASE ('AROME ','MERCATOR ')
200  ! 3.1 lambert conformal projection(aladin files)
201  ! or mercator projection(aladin reunion files)
202  CALL grib_get(igrib,'Nj',ny,iret)
203  CALL grib_get(igrib,'Ni',nx,iret)
204  nni= nx * ny
205  !
206  !
207  CASE ('GAUSS ','ROTGAUSS ','LATLON ')
208  ! 3.2 usual or gaussian lat,lon grid(ecmwf files)
209  !
210  CALL grib_get(igrib,'Nj',ninla,iret)
211  ALLOCATE (ninlo(ninla))
212  ALLOCATE (ininlo_grib(ninla))
213  nilen = 0
214  CALL grib_is_missing(igrib,'pl',imissing,iret)
215  IF (iret /= 0 .OR. imissing==1) THEN ! regular
216  CALL grib_get(igrib,'Ni',ininlo_grib(1),iret)
217  ininlo_grib(2:ninla)=ininlo_grib(1)
218  nilen=ninla*ininlo_grib(1)
219  ELSE ! quasi-regular
220  CALL grib_get(igrib,'pl',ininlo_grib)
221  DO jloop1=1 ,ninla
222  nilen = nilen + ininlo_grib(jloop1)
223  ENDDO
224  ENDIF
225  nni = nilen
226  ninlo = ininlo_grib !JUAN
227  CASE ('ROTLATLON')
228  CALL grib_get(igrib,'Nj',nry,iret)
229  CALL grib_get(igrib,'Ni',nrx,iret)
230  nni = nrx * nry
231 
232  CASE default
233  CALL abor1_sfx('PREP_GRIB_GRID: GRID PROJECTION NOT SUPPORTED')
234  !
235 END SELECT
236 !
237 !---------------------------------------------------------------------------------------
238 !* 4. Updates grid information
239 !---------------------------------------------------------------------------------------
240 !
241 xx = xundef
242 xy = xundef
243 xila1 = xundef
244 xilo1 = xundef
245 xila2 = xundef
246 xilo2 = xundef
247 lrotpole = .false.
248 xcoef = xundef
249 xlap = xundef
250 xlop = xundef
251 
252 SELECT CASE (hgridtype)
253 
254  CASE ('AROME ')
255  ! 4.1 lambert conformal projection(aladin files)
256  !
257  CALL grib_get(igrib,'xDirectionGridLength',ilenx)
258  CALL grib_get(igrib,'yDirectionGridLength',ileny)
259  xy = (ny-1)*ileny
260  xx = (nx-1)*ilenx
261 
262  CALL grib_get(igrib,'Latin1InDegrees',xlat0)
263  CALL grib_get(igrib,'LoVInDegrees',xlon0)
264  IF (xlon0 > 180.) xlon0 = xlon0 - 360.
265 
266  CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xlator)
267  CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xlonor)
268  IF (xlonor > 180.) xlonor = xlonor - 360.
269 
270  xrpk = sin(xlat0/180.*xpi)
271  xbeta = 0.
272  !
273  CASE ('GAUSS ','LATLON ')
274  hgridtype = 'GAUSS '
275  ! 4.2 usual or gaussian lat,lon grid(ecmwf files)
276  ! no projection - just stores the grid definition
277  !
278  CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xila1)
279  CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xilo1)
280  CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xila2)
281  CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xilo2)
282 
283  lrotpole = .false.
284  !
285  CASE ('ROTLATLON ')
286  !
287  ! 4.2.5 rotated lat/lon grid(hirlam)
288  !
289  CALL grib_get(igrib,'iScansNegatively',iscan)
290  CALL grib_get(igrib,'jScansNegatively',jscan)
291 
292  IF (iscan+jscan == 0 ) THEN
293  CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xrila2)
294  CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xrilo1)
295  CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xrila1)
296  CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xrilo2)
297  ELSEIF (iscan+jscan == 2) THEN
298  CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xrila1)
299  CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xrilo2)
300  CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xrila2)
301  CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xrilo1)
302  ELSEIF (iscan == 1) THEN
303  CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xrila2)
304  CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xrilo2)
305  CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xrila1)
306  CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xrilo1)
307  ELSEIF (jscan == 1) THEN
308  CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xrila1)
309  CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xrilo1)
310  CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xrila2)
311  CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xrilo2)
312  ENDIF
313 
314  CALL grib_get(igrib,'latitudeOfSouthernPoleInDegrees',xrlap)
315  CALL grib_get(igrib,'longitudeOfSouthernPoleInDegrees',xrlop)
316 
317  CALL grib_get(igrib,'iDirectionIncrementInDegrees',xrdx)
318  CALL grib_get(igrib,'jDirectionIncrementInDegrees',xrdy)
319 
320  WRITE(kluout,*)'XRILA1,XRILO1',xrila1,xrilo1
321  WRITE(kluout,*)'XRILA2,XRILO2',xrila2,xrilo2
322  WRITE(kluout,*)'XRLAP,XRLOP',xrlap,xrlop
323  WRITE(kluout,*)'XRDX,XRDY',xrdx,xrdy
324  !
325  CASE ('ROTGAUSS ')
326  ! 4.3 stretched lat,lon grid(arpege files)
327  !
328  hgridtype = 'GAUSS '
329  CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xila1)
330  CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xilo1)
331  CALL grib_get(igrib,'latitudeOfLastGridPointInDegrees',xila2)
332  CALL grib_get(igrib,'longitudeOfLastGridPointInDegrees',xilo2)
333 
334  lrotpole = .true.
335  CALL grib_get(igrib,'stretchingFactor',xcoef)
336  CALL grib_get(igrib,'latitudeOfStretchingPoleInDegrees',xlap)
337  CALL grib_get(igrib,'longitudeOfStretchingPoleInDegrees',xlop)
338  !
339 
340  CASE ('MERCATOR ')
341  ! 4.4 mercator projection(aladin reunion files)
342  !
343  hgridtype = 'AROME '
344  CALL grib_get(igrib,'Dj',ileny)
345  CALL grib_get(igrib,'Di',ilenx)
346  xy = (ny-1)*ileny
347  xx = (nx-1)*ilenx
348 
349  CALL grib_get(igrib,'LaDInDegrees',xlat0)
350  xlon0 = 0.
351 
352  CALL grib_get(igrib,'latitudeOfFirstGridPointInDegrees',xlator)
353  CALL grib_get(igrib,'longitudeOfFirstGridPointInDegrees',xlonor)
354  IF (xlonor > 180.) xlonor = xlonor - 360.
355 
356  xrpk = 0.
357  xbeta = 0.
358 
359  CASE default
360  WRITE (kluout,'(A)') 'No such projection implemented in prep_grib_grid ',hgrid
361  CALL abor1_sfx('PREP_GRIB_GRID: UNKNOWN PROJECTION')
362  !
363 END SELECT
364 !---------------------------------------------------------------------------------------
365 !* 2.4 Read date
366 !---------------------------------------------------------------------------------------
367 !
368 WRITE (kluout,'(A)') ' | Reading date'
369 !
370  CALL grib_get(igrib,'year',tptime_grib%TDATE%YEAR,iret)
371  CALL grib_get(igrib,'month',tptime_grib%TDATE%MONTH,iret)
372  CALL grib_get(igrib,'day',tptime_grib%TDATE%DAY,iret)
373  CALL grib_get(igrib,'time',itime,iret)
374 tptime_grib%TIME=int(itime/100)*3600+(itime-int(itime/100)*100)*60
375 !
376  CALL grib_get(igrib,'P1',ip1,iret)
377 IF ( ip1>0 ) THEN
378  CALL grib_get(igrib,'unitOfTimeRange',iunittime,iret)
379  SELECT CASE (iunittime) ! Time unit indicator
380  CASE (1) !hour
381  tptime_grib%TIME = tptime_grib%TIME + ip1*3600.
382  CASE (0) !minute
383  tptime_grib%TIME = tptime_grib%TIME + ip1*60.
384  END SELECT
385 ENDIF
386 !
387 !---------------------------------------------------------------------------------------
388 !
389  CALL grib_release(igrib,iret)
390 IF (iret /= 0) THEN
391  CALL abor1_sfx('PREP_GRIB_GRID: Error in releasing the grib message memory')
392 END IF
393 !
394 IF (lhook) CALL dr_hook('PREP_GRIB_GRID',1,zhook_handle)
395 !
396 END SUBROUTINE prep_grib_grid
modd_grid_grib
Definition: modd_grid_grib.F90:6
abor1_sfx
subroutine abor1_sfx(YTEXT)
Definition: abor1_sfx.F90:6
modd_grid_rotlatlon
Definition: modd_grid_rotlatlon.F90:6
modd_type_date_surf::date_time
Definition: modd_type_date_surf.F90:46
modd_grid_gauss
Definition: modd_grid_gauss.F90:6
modd_surf_par
Definition: modd_surf_par.F90:6
modd_type_date_surf
Definition: modd_type_date_surf.F90:6
mode_read_grib
Definition: mode_read_grib.F90:6
prep_grib_grid
subroutine prep_grib_grid(HGRIB, KLUOUT, HINMODEL, HGRIDTYPE, TPTIME_GRIB)
Definition: prep_grib_grid.F90:6
modd_csts
Definition: modd_csts.F90:6
modd_grid_arome
Definition: modd_grid_arome.F90:6
Generated on Fri May 13 2016 09:31:05 for SURFEX by   doxygen 1.8.5