SURFEX v8.1
General documentation of Surfex
mode_read_grib.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 MODULE mode_read_grib
7 ! #####################
8 !-------------------------------------------------------------------
9 !
10 USE modi_abor1_sfx
11 USE grib_api
12 !
13 USE yomhook ,ONLY : lhook, dr_hook
14 USE parkind1 ,ONLY : jprb
15 !
16 CONTAINS
17 !-------------------------------------------------------------------
18 ! ####################
19  SUBROUTINE make_grib_index(HGRIB)
20 ! ####################
21 !
22 USE modd_grid_grib, ONLY : cgrib_file, nidx
23 !
24 IMPLICIT NONE
25 !
26  CHARACTER(LEN=*), INTENT(IN) :: HGRIB
27 !
28 INTEGER(KIND=kindOfInt) :: IRET
29 REAL(KIND=JPRB) :: ZHOOK_HANDLE
30 !
31 IF (lhook) CALL dr_hook('MODE_READ_GRIB:MAKE_GRIB_INDEX',0,zhook_handle)
32 !
33 IF (cgrib_file==hgrib) CALL dr_hook('MODE_READ_GRIB:MAKE_GRIB_INDEX',1,zhook_handle)
34 IF (cgrib_file==hgrib) RETURN
35 !
36 cgrib_file=hgrib
37 !
38  CALL grib_index_create(nidx,hgrib,'indicatorOfParameter',iret)
39 IF (iret/=0) CALL abor1_sfx("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while creating the grib index")
40 !
41 IF (lhook) CALL dr_hook('MODE_READ_GRIB:MAKE_GRIB_INDEX',1,zhook_handle)
42 !
43 END SUBROUTINE make_grib_index
44 !-------------------------------------------------------------------
45 ! ####################
46  SUBROUTINE clear_grib_index
47 ! ####################
48 !
49 USE modd_grid_grib, ONLY : cgrib_file, nidx
50 !
51 IMPLICIT NONE
52 !
53 INTEGER(KIND=kindOfInt) :: IRET
54 !
55 REAL(KIND=JPRB) :: ZHOOK_HANDLE
56 !
57 IF (lhook) CALL dr_hook('MODE_READ_GRIB:CLEAR_GRIB_INDEX',0,zhook_handle)
58 !
59 IF (cgrib_file.NE."") THEN
60  cgrib_file=""
61  CALL grib_index_release(nidx,iret)
62  IF (iret/=0) CALL abor1_sfx("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while deleting the grib index")
63 ENDIF
64 !
65 IF (lhook) CALL dr_hook('MODE_READ_GRIB:CLEAR_GRIB_INDEX',1,zhook_handle)
66 !
67 END SUBROUTINE clear_grib_index
68 !-------------------------------------------------------------------
69 ! ####################
70  SUBROUTINE get_grib_message(KLUOUT,KLTYPE,KLEV1,KLEV2,KGRIB,KFOUND)
71 ! ####################
72 !
73 USE modd_grid_grib, ONLY : nidx
74 !
75 IMPLICIT NONE
76 !
77 INTEGER, INTENT(IN) :: KLUOUT
78 INTEGER, INTENT(INOUT) :: KLTYPE
79 INTEGER, INTENT(INOUT) :: KLEV1
80 INTEGER, INTENT(INOUT) :: KLEV2
81 INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KGRIB
82 INTEGER, INTENT(OUT) :: KFOUND
83 !
84 INTEGER :: ILTYPE
85 INTEGER :: ILEV1
86 INTEGER :: ILEV2
87 INTEGER(KIND=kindOfInt) :: IRET
88 !
89 REAL(KIND=JPRB) :: ZHOOK_HANDLE
90 !
91 IF (lhook) CALL dr_hook('MODE_READ_GRIB:GET_GRIB_MESSAGE',0,zhook_handle)
92 !
93 iret = 0
94 kfound=0
95 !
96 DO WHILE (iret /= grib_end_of_index .AND. kfound/=3)
97  !
98  iret = 0
99  kfound=0
100  !
101  IF (kltype/=-2) THEN
102  CALL grib_get(kgrib,'indicatorOfTypeOfLevel',iltype,iret)
103  CALL test_iret(kluout,iltype,kltype,iret)
104  ENDIF
105  !
106  IF (iret.EQ.0) THEN
107  !
108  kfound = kfound + 1
109  !
110  IF (klev1/=-2) THEN
111  CALL grib_get(kgrib,'topLevel',ilev1,iret)
112  CALL test_iret(kluout,ilev1,klev1,iret)
113  ENDIF
114  !
115  IF (iret.EQ.0) THEN
116  !
117  kfound = kfound + 1
118  !
119  IF (klev2/=-2) THEN
120  CALL grib_get(kgrib,'bottomLevel',ilev2,iret)
121  CALL test_iret(kluout,ilev2,klev2,iret)
122  ENDIF
123  !
124  IF (iret.EQ.0) kfound = kfound + 1
125  !
126  ENDIF
127  !
128  ENDIF
129  !
130  IF (kfound.NE.3) THEN
131  CALL grib_release(kgrib)
132  CALL grib_new_from_index(nidx,kgrib,iret)
133  ENDIF
134  !
135 ENDDO
136 !
137 IF (lhook) CALL dr_hook('MODE_READ_GRIB:GET_GRIB_MESSAGE',1,zhook_handle)
138 !
139 CONTAINS
140 !
141 ! ##############
142  SUBROUTINE test_iret(KLUOUT,VAL1,VAL0,KRET)
143 ! ##############
144 !
145 IMPLICIT NONE
146 !
147 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
148 INTEGER, INTENT(IN) :: VAL1
149 INTEGER, INTENT(INOUT) :: VAL0
150 INTEGER(KIND=kindOfInt), INTENT(INOUT) :: KRET ! number of the message researched
151 !
152 REAL(KIND=JPRB) :: ZHOOK_HANDLE
153 !
154 IF (lhook) CALL dr_hook('MODE_READ_GRIB:TEST_IRET',0,zhook_handle)
155 !
156 IF (kret > 0) THEN
157  WRITE (kluout,'(A)')' | Error encountered in the Grib file, skipping field'
158 ELSE IF (kret == -6) THEN
159  WRITE (kluout,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
160 ELSEIF (val1 /= val0) THEN
161  IF (val0 == -1) THEN
162  val0 = val1
163  ELSE
164  kret=1
165  ENDIF
166 ENDIF
167 !
168 IF (lhook) CALL dr_hook('MODE_READ_GRIB:TEST_IRET',1,zhook_handle)
169 END SUBROUTINE test_iret
170 !
171 END SUBROUTINE get_grib_message
172 !-------------------------------------------------------------------
173 ! ####################
174  SUBROUTINE read_grib(HGRIB,KLUOUT,KPARAM,KRET,PFIELD,KLTYPE,KLEV1,KLEV2,KPARAM2)
175 ! ####################
176 !
177 USE modd_grid_grib, ONLY : nidx
178 !
179 IMPLICIT NONE
180 !
181  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! name of the GRIB file
182 INTEGER, INTENT(IN) :: KLUOUT
183 INTEGER,INTENT(IN) :: KPARAM ! Parameter to read
184 INTEGER(KIND=kindOfInt), INTENT(OUT) :: KRET
185 REAL, DIMENSION(:), POINTER :: PFIELD
186 INTEGER,INTENT(INOUT), OPTIONAL :: KLTYPE ! Level type
187 INTEGER,INTENT(INOUT), OPTIONAL :: KLEV1 ! Level parameter 1
188 INTEGER,INTENT(INOUT), OPTIONAL :: KLEV2 ! Level parameter 2
189 INTEGER, INTENT(INOUT), OPTIONAL :: KPARAM2
190 !
191 INTEGER :: ILTYPE, ILEV1, ILEV2
192 INTEGER(KIND=kindOfInt) :: IGRIB
193 INTEGER :: ISIZE, IFOUND
194 REAL(KIND=JPRB) :: ZHOOK_HANDLE
195 !
196 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB',0,zhook_handle)
197 !
198 iltype=-2
199 IF (PRESENT(kltype)) iltype=kltype
200 ilev1=-2
201 IF (PRESENT(klev1)) ilev1=klev1
202 ilev2=-2
203 IF (PRESENT(klev2)) ilev2=klev2
204 !
205  CALL make_grib_index(hgrib)
206 !
207 ifound=0
208 kret=0
209 !
210  CALL grib_index_select(nidx,'indicatorOfParameter',kparam,kret)
211  CALL grib_new_from_index(nidx,igrib,kret)
212 IF (kret.EQ.0) CALL get_grib_message(kluout,iltype,ilev1,ilev2,igrib,ifound)
213 !
214 IF (PRESENT(kparam2)) THEN
215  IF (ifound/=3) THEN
216  CALL grib_index_select(nidx,'indicatorOfParameter',kparam2,kret)
217  CALL grib_new_from_index(nidx,igrib,kret)
218  IF (kret.EQ.0) THEN
219  iltype=-2
220  IF (PRESENT(kltype)) iltype=kltype
221  CALL get_grib_message(kluout,iltype,ilev1,ilev2,igrib,ifound)
222  ENDIF
223  ELSE
224  kparam2 = kparam
225  ENDIF
226 ENDIF
227 !
228 IF (ifound==3) THEN
229  !
230  IF (PRESENT(kltype)) kltype = iltype
231  IF (PRESENT(klev1)) klev1 = ilev1
232  IF (PRESENT(klev2)) klev2 = ilev2
233  !
234  IF (.NOT.ASSOCIATED(pfield)) THEN
235  CALL grib_get_size(igrib,'values',isize,kret)
236  IF (kret.NE.0) CALL abor1_sfx("MODE_READ_GRIB:READ_GRIB: Problem getting size of values")
237  ALLOCATE(pfield(isize))
238  ENDIF
239  !
240  CALL grib_get(igrib,'values',pfield,kret)
241  IF (kret.NE.0) CALL abor1_sfx("MODE_READ_GRIB:READ_GRIB: Problem getting values")
242  CALL grib_release(igrib,kret)
243  IF (kret.NE.0) CALL abor1_sfx("MODE_READ_GRIB:READ_GRIB: Problem releasing memory")
244  !
245 ELSE
246  !
247  kret = 1
248  !
249 ENDIF
250 !
251 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB',1,zhook_handle)
252 END SUBROUTINE read_grib
253 !
254 !-------------------------------------------------------------------
255 !-------------------------------------------------------------------
256 ! ####################
257  SUBROUTINE read_grib_land_mask(HGRIB,KLUOUT,HINMODEL,PMASK)
258 ! ####################
259 !
260 IMPLICIT NONE
261 !
262  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
263 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
264  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
265 REAL, DIMENSION(:), POINTER :: PMASK ! Land mask
266 !
267 INTEGER(KIND=kindOfInt) :: IRET ! return code
268 INTEGER :: ILTYPE ! leveltype
269 INTEGER :: ILEV ! level
270 REAL(KIND=JPRB) :: ZHOOK_HANDLE
271 !-------------------------------------------------------------------
272 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_LAND_MASK',0,zhook_handle)
273 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_LAND_MASK: | Reading land mask from ',hinmodel
274 !
275 SELECT CASE (hinmodel)
276  CASE ('ECMWF ')
277  CALL read_grib(hgrib,kluout,172,iret,pmask)
278  CASE ('ARPEGE','ALADIN','MOCAGE')
279  CALL read_grib(hgrib,kluout,81,iret,pmask)
280  CASE ('HIRLAM')
281  iltype=105
282  ilev =0
283  CALL read_grib(hgrib,kluout,81,iret,pmask,kltype=iltype,klev1=ilev)
284  CASE DEFAULT
285  CALL abor1_sfx('MODE_READ_GRIB:READ_GRIB_LAND_MASK: OPTION NOT SUPPORTED '//hinmodel)
286 END SELECT
287 !
288 IF (iret /= 0) THEN
289  CALL abor1_sfx('MODE_READ_GRIB: LAND SEA MASK MISSING (READ_GRIB_LAND_MASK)')
290 END IF
291 !
292 WHERE (pmask>0.5)
293  pmask = 1.
294 ELSEWHERE
295  pmask = 0.
296 END WHERE
297 !
298 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_LAND_MASK',1,zhook_handle)
299 END SUBROUTINE read_grib_land_mask
300 !-------------------------------------------------------------------
301 ! ############################
302  SUBROUTINE read_grib_zs(HGRIB,KLUOUT,HINMODEL,PZS)
303 ! ############################
304 !
305 USE modd_csts, ONLY : xg
306 !
307 IMPLICIT NONE
308 !
309  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
310 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
311  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
312 REAL, DIMENSION(:), POINTER :: PZS !
313 !
314 INTEGER(KIND=kindOfInt) :: IRET ! return code
315 REAL(KIND=JPRB) :: ZHOOK_HANDLE
316 !-------------------------------------------------------------------
317 !* Read orography
318 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_ZS',0,zhook_handle)
319 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_ZS: | Reading orography from ',hinmodel
320 !
321 SELECT CASE (hinmodel)
322  CASE ('ECMWF ')
323  CALL read_grib(hgrib,kluout,129,iret,pzs)
324  CASE ('ARPEGE','MOCAGE')
325  CALL read_grib(hgrib,kluout,8,iret,pzs)
326  CASE ('HIRLAM','ALADIN')
327  CALL read_grib(hgrib,kluout,6,iret,pzs)
328  CASE DEFAULT
329  CALL abor1_sfx('MODE_READ_GRIB:READ_GRIB_ZS:OPTION NOT SUPPORTED '//hinmodel)
330 END SELECT
331 !
332 IF (iret /= 0) THEN
333  CALL abor1_sfx('MODE_READ_GRIB: OROGRAPHY MISSING (READ_GRIB_ZS_LAND)')
334 END IF
335 !
336 ! Datas given in archives are multiplied by the gravity acceleration
337 pzs = pzs / xg
338 !
339 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_ZS',1,zhook_handle)
340 END SUBROUTINE read_grib_zs
341 !-------------------------------------------------------------------
342 ! ############################
343  SUBROUTINE read_grib_zs_land(HGRIB,KLUOUT,HINMODEL,PMASK,PZSL)
344 ! ############################
345 !
346 IMPLICIT NONE
347 !
348  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
349 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
350  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
351 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
352 REAL, DIMENSION(:), POINTER :: PZSL !
353 !
354 REAL(KIND=JPRB) :: ZHOOK_HANDLE
355 !-------------------------------------------------------------------
356 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_ZS_LAND',0,zhook_handle)
357 !
358  CALL read_grib_zs(hgrib,kluout,hinmodel,pzsl)
359 !
360 IF (SIZE(pmask)==SIZE(pzsl)) &
361  WHERE (pmask(:)/=1.) pzsl = 0.
362 !
363 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_ZS_LAND',1,zhook_handle)
364 END SUBROUTINE read_grib_zs_land
365 !-------------------------------------------------------------------
366 ! ############################
367  SUBROUTINE read_grib_zs_sea(HGRIB,KLUOUT,HINMODEL,PMASK,PZSS)
368 ! ############################
369 !
370 IMPLICIT NONE
371 !
372  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
373 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
374  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
375 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
376 REAL, DIMENSION(:), POINTER :: PZSS !
377 !
378 REAL(KIND=JPRB) :: ZHOOK_HANDLE
379 !-------------------------------------------------------------------
380 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_ZS_SEA',0,zhook_handle)
381 !
382  CALL read_grib_zs(hgrib,kluout,hinmodel,pzss)
383 !
384 IF (SIZE(pmask)==SIZE(pzss)) &
385  WHERE (pmask(:)/=0.) pzss = 0.
386 !
387 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_ZS_SEA',1,zhook_handle)
388 END SUBROUTINE read_grib_zs_sea
389 !-------------------------------------------------------------------
390 ! ###########################
391  SUBROUTINE read_grib_t(HGRIB,KLUOUT,HINMODEL,PT)
392 ! ###########################
393 !
394 IMPLICIT NONE
395 !
396  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
397 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
398  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
399 REAL, DIMENSION(:), POINTER :: PT !
400 !
401 INTEGER(KIND=kindOfInt) :: IRET ! return code
402 INTEGER :: ILTYPE ! type of level (Grib code table 3)
403 INTEGER :: ILEV ! level definition
404 REAL(KIND=JPRB) :: ZHOOK_HANDLE
405 !-------------------------------------------------------------------
406 !* Read surface temperature
407 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_T',0,zhook_handle)
408 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_T: | Reading surface temperature'
409 !
410 SELECT CASE (hinmodel)
411  CASE ('ECMWF ')
412  CALL read_grib(hgrib,kluout,139,iret,pt)
413 
414  CASE ('ARPEGE','ALADIN','MOCAGE')
415  ilev=0
416  iltype=111
417  CALL read_grib(hgrib,kluout,11,iret,pt,kltype=iltype,klev1=ilev)
418  IF (iret /= 0) THEN
419  iltype=1
420  CALL read_grib(hgrib,kluout,11,iret,pt,kltype=iltype)
421  IF (iret /= 0) THEN
422  iltype=105
423  CALL read_grib(hgrib,kluout,11,iret,pt,kltype=iltype,klev1=ilev)
424  ENDIF
425  END IF
426 
427  CASE ('HIRLAM ')
428  WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_T: | Reading surface temperature tile 4'
429  iltype=105
430  ilev=904
431  CALL read_grib(hgrib,kluout,11,iret,pt,kltype=iltype,klev1=ilev)
432 
433  CASE DEFAULT
434  CALL abor1_sfx('MODE_READ_GRIB:READ_GRIB_T:OPTION NOT SUPPORTED '//hinmodel)
435 END SELECT
436 !
437 IF (iret /= 0) THEN
438  CALL abor1_sfx('MODE_READ_GRIB: SURFACE TEMPERATURE MISSING (READ_GRIB_T)')
439 END IF
440 !
441 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_T',1,zhook_handle)
442 END SUBROUTINE read_grib_t
443 !-------------------------------------------------------------------
444 ! ###########################
445  SUBROUTINE read_grib_ts(HGRIB,KLUOUT,HINMODEL,PMASK,PTS)
446 ! ###########################
447 !
448 USE modd_surf_par, ONLY : xundef
449 !
450 IMPLICIT NONE
451 !
452  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
453 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
454  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
455 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
456 REAL, DIMENSION(:), POINTER :: PTS !
457 !
458 REAL(KIND=JPRB) :: ZHOOK_HANDLE
459 !-------------------------------------------------------------------
460 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TS',0,zhook_handle)
461 !
462  CALL read_grib_t(hgrib,kluout,hinmodel,pts)
463 !
464 IF (SIZE(pmask)==SIZE(pts)) WHERE (pmask(:)/=1.) pts = xundef
465 !
466 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TS',1,zhook_handle)
467 END SUBROUTINE read_grib_ts
468 !-------------------------------------------------------------------
469 ! ###########################
470  SUBROUTINE read_grib_sst(HGRIB,KLUOUT,HINMODEL,PMASK,PSST)
471 ! ###########################
472 !
473 USE modd_surf_par, ONLY : xundef
474 !
475 IMPLICIT NONE
476 !
477  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
478 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
479  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
480 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
481 REAL, DIMENSION(:), POINTER :: PSST !
482 !
483 INTEGER :: IRET
484 REAL(KIND=JPRB) :: ZHOOK_HANDLE
485 !-------------------------------------------------------------------
486 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SST',0,zhook_handle)
487 !
488 SELECT CASE (hinmodel)
489  CASE ('ECMWF ')
490  CALL read_grib(hgrib,kluout,34,iret,psst)
491  IF (iret /= 0) CALL read_grib_t(hgrib,kluout,hinmodel,psst)
492  CASE ('ARPEGE','ALADIN','MOCAGE','HIRLAM')
493  CALL read_grib_t(hgrib,kluout,hinmodel,psst)
494  CASE DEFAULT
495  CALL abor1_sfx('MODE_READ_GRIB:READ_GRIB_SST:OPTION NOT SUPPORTED '//hinmodel)
496 END SELECT
497 !
498 IF (SIZE(pmask)==SIZE(psst)) WHERE (pmask(:)/=0.) psst = xundef
499 !
500 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SST',1,zhook_handle)
501 END SUBROUTINE read_grib_sst
502 !-------------------------------------------------------------------
503 ! ###########################
504  SUBROUTINE read_grib_tswater(HGRIB,KLUOUT,HINMODEL,PMASK,PTS)
505 ! ###########################
506 !
507 USE modd_surf_par, ONLY : xundef
508 !
509 IMPLICIT NONE
510 !
511  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
512 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
513  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
514 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
515 REAL, DIMENSION(:), POINTER :: PTS !
516 !
517 INTEGER :: IRET
518 REAL(KIND=JPRB) :: ZHOOK_HANDLE
519 !-------------------------------------------------------------------
520 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TSWATER',0,zhook_handle)
521 !
522 SELECT CASE (hinmodel)
523  CASE ('ECMWF ')
524  CALL read_grib(hgrib,kluout,3080,iret,pts)
525  IF (iret /= 0) CALL read_grib_t2(hgrib,kluout,hinmodel,pmask,pts)
526  CASE ('ARPEGE','ALADIN','MOCAGE','HIRLAM')
527  CALL read_grib_t2(hgrib,kluout,hinmodel,pmask,pts)
528  CASE DEFAULT
529  CALL abor1_sfx('MODE_READ_GRIB:READ_GRIB_TSWATER:OPTION NOT SUPPORTED '//hinmodel)
530 END SELECT
531 !
532 IF (SIZE(pmask)==SIZE(pts)) WHERE (pmask(:)/=0.) pts = xundef
533 !
534 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TSWATER',1,zhook_handle)
535 END SUBROUTINE read_grib_tswater
536 !-------------------------------------------------------------------
537 ! ###########################
538  SUBROUTINE read_grib_t2(HGRIB,KLUOUT,HINMODEL,PMASK,PT2)
539 ! ###########################
540 !
541 USE modd_surf_par, ONLY : xundef
542 !
543 IMPLICIT NONE
544 !
545  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
546 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
547  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
548 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
549 REAL, DIMENSION(:), POINTER :: PT2 !
550 !
551 INTEGER(KIND=kindOfInt) :: IRET
552 INTEGER :: ILTYPE, ILEV ! type of level (Grib code table 3)
553 REAL(KIND=JPRB) :: ZHOOK_HANDLE
554 !-------------------------------------------------------------------
555 !* Read deep soil temperature
556 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_T2',0,zhook_handle)
557 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_T2: | Reading deep soil temperature'
558 !
559 SELECT CASE (hinmodel)
560  CASE ('ECMWF ')
561  CALL read_grib(hgrib,kluout,170,iret,pt2)
562  CASE ('ARPEGE','ALADIN','MOCAGE')
563  iltype=111
564  CALL read_grib(hgrib,kluout,11,iret,pt2,kltype=iltype)
565  IF (iret /= 0) THEN
566  iltype=105
567  CALL read_grib(hgrib,kluout,11,iret,pt2,kltype=iltype)
568  ENDIF
569  CASE ('HIRLAM ')
570  iltype=105
571  ilev=954
572  CALL read_grib(hgrib,kluout,11,iret,pt2,kltype=iltype,klev1=ilev)
573  CASE DEFAULT
574  CALL abor1_sfx('MODE_READ_GRIB:READ_GRIB_T2:OPTION NOT SUPPORTED '//hinmodel)
575 END SELECT
576 !
577 IF (iret /= 0) THEN
578  CALL abor1_sfx('MODE_READ_GRIB: DEEP SOIL TEMPERATURE MISSING (READ_GRIB_T2)')
579 END IF
580 !
581 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_T2',1,zhook_handle)
582 END SUBROUTINE read_grib_t2
583 !-------------------------------------------------------------------
584 !-------------------------------------------------------------------
585 ! ###########################
586  SUBROUTINE read_grib_t2_land(HGRIB,KLUOUT,HINMODEL,PMASK,ZFIELD)
587 ! ###########################
588 !
589 USE modd_surf_par, ONLY : xundef
590 !
591 IMPLICIT NONE
592 !
593  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
594 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
595  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
596 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
597 REAL, DIMENSION(:), POINTER :: ZFIELD !
598 !
599 REAL(KIND=JPRB) :: ZHOOK_HANDLE
600 !-------------------------------------------------------------------
601 !* Read deep soil temperature
602 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_T2_LAND',0,zhook_handle)
603 !
604  CALL read_grib_t2(hgrib,kluout,hinmodel,pmask,zfield)
605 !
606 IF (SIZE(pmask)==SIZE(zfield)) WHERE (pmask(:)/=1.) zfield = xundef
607 !
608 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_T2_LAND',1,zhook_handle)
609 END SUBROUTINE read_grib_t2_land
610 !-------------------------------------------------------------------
611 !-------------------------------------------------------------------
612 SUBROUTINE put_layer_depth(KLUOUT,KLEV,HROUT,KLTYPE,KLEV1,KLEV2,KNLAYERDEEP,PV4,PV,PD)
613 !
614 IMPLICIT NONE
615 !
616 INTEGER, INTENT(IN) :: KLUOUT
617 INTEGER, INTENT(IN) :: KLEV
618  CHARACTER(LEN=*), INTENT(IN) :: HROUT
619 INTEGER, INTENT(INOUT) :: KLTYPE
620 INTEGER, INTENT(IN) :: KLEV1
621 INTEGER, INTENT(IN) :: KLEV2
622 INTEGER, INTENT(IN) :: KNLAYERDEEP
623 REAL, INTENT(IN) :: PV4
624 REAL, INTENT(IN) :: PV
625 REAL, INTENT(OUT) :: PD
626 REAL(KIND=JPRB) :: ZHOOK_HANDLE
627 !
628 IF (lhook) CALL dr_hook('MODE_READ_GRIB:PUT_LAYER_DEPTH',0,zhook_handle)
629 !
630 IF (klev2 == -1) kltype = 0
631 IF (kltype==112) THEN
632  pd = (klev2 - klev1) / 100.
633 ELSE
634  IF (knlayerdeep == 4) THEN
635  pd = pv4
636  ELSE
637  pd = pv
638  END IF
639  WRITE (kluout,'(A,i1,A,f5.2,A)') 'MODE_READ_GRIB:'//trim(hrout)//': | Level ',&
640  klev,' height not available, assuming ',pd,' m'
641 END IF
642 !
643 IF (lhook) CALL dr_hook('MODE_READ_GRIB:PUT_LAYER_DEPTH',1,zhook_handle)
644 END SUBROUTINE put_layer_depth
645 !-------------------------------------------------------------------
646 ! #######################
647  SUBROUTINE fill_pfield(KLUOUT,HROUT,KNLAYERDEEP,PDIN,PFIELDIN,PMASK,PFIELDOUT,PDOUT)
648 ! #######################
649 !
650 USE modd_surf_par, ONLY : xundef
651 !
652 IMPLICIT NONE
653 !
654 INTEGER, INTENT(IN) :: KLUOUT
655  CHARACTER(LEN=*), INTENT(IN) :: HROUT
656 INTEGER, INTENT(IN) :: KNLAYERDEEP
657 REAL, DIMENSION(:), INTENT(IN) :: PDIN
658 REAL, DIMENSION(:,:), INTENT(IN) :: PFIELDIN
659 REAL, DIMENSION(:), INTENT(IN) :: PMASK
660 REAL, DIMENSION(:,:), POINTER :: PFIELDOUT
661 REAL, DIMENSION(:,:), POINTER :: PDOUT
662 !
663  CHARACTER(LEN=20) :: FMT0
664 INTEGER :: JL
665 REAL(KIND=JPRB) :: ZHOOK_HANDLE
666 !
667 IF (lhook) CALL dr_hook('MODE_READ_GRIB:FILL_PFIELD',0,zhook_handle)
668 !--------------------------------------------------------------------------------
669 ! 1. Display the number of layer found
670 ! -----------------------
671 WRITE(fmt0,fmt='(A8,I1,A11)') '(A,I1,A,',knlayerdeep,'(F5.2,","))'
672 WRITE (kluout,fmt=fmt0) 'MODE_READ_GRIB:'//trim(hrout)//': | ',knlayerdeep,&
673  ' deep layers, heights are : ',pdin(1:knlayerdeep)
674 !--------------------------------------------------------------------------------
675 ! 2. Set temperature profile and layer thicknesses
676 ! -----------------------------------------------
677 ALLOCATE(pfieldout(SIZE(pfieldin,1),knlayerdeep))
678 ALLOCATE(pdout(SIZE(pfieldin,1),knlayerdeep))
679 !
680 DO jl=1,knlayerdeep
681  pdout(:,jl)=sum(pdin(1:jl))
682  pfieldout(:,jl)=pfieldin(:,jl)
683  IF (SIZE(pmask)==SIZE(pfieldout,1)) &
684  WHERE (pmask(:)/=1.) pfieldout(:,jl) = xundef
685 ENDDO
686 !
687 IF (lhook) CALL dr_hook('MODE_READ_GRIB:FILL_PFIELD',1,zhook_handle)
688 END SUBROUTINE fill_pfield
689 !-------------------------------------------------------------------
690 ! #######################
691  SUBROUTINE read_grib_tg_ecmwf(HGRIB,KLUOUT,HINMODEL,PMASK,PTG,PD)
692 ! #######################
693 !
694 IMPLICIT NONE
695 !
696 !* dummy arguments
697 ! ---------------
698  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
699 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
700  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
701 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
702 REAL, DIMENSION(:,:), POINTER :: PTG ! field to initialize
703 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
704 !
705 !* local variables
706 ! ---------------
707 INTEGER(KIND=kindOfInt) :: IRET ! return code
708 INTEGER :: ILTYPE ! type of level (Grib code table 3)
709 INTEGER :: ILEV1 ! level definition
710 INTEGER :: ILEV2 ! level definition
711 INTEGER :: JL ! layer loop counter
712 INTEGER :: INLAYERDEEP! number of deep moisture layers
713 REAL, DIMENSION(:), POINTER :: ZFIELD => null() ! first layer temperature
714 REAL, DIMENSION(:,:), ALLOCATABLE:: ZTG ! first layer temperature
715 REAL, DIMENSION(:) , ALLOCATABLE:: ZD
716 REAL(KIND=JPRB) :: ZHOOK_HANDLE
717 !--------------------------------------------------------------------------------
718 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TG_ECMWF',0,zhook_handle)
719 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_TG_ECMWF: | Reading soil temperature'
720 !
721 ALLOCATE(zd(5))
722 !
723 ! 1. Search and read level 1 (and its depth)
724 ! --------------------------------------
725 iltype= -1
726 ilev1 = -1
727 ilev2 = -1
728  CALL read_grib(hgrib,kluout,139,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
729 !
730 IF (iret== 0) THEN
731  CALL put_layer_depth(kluout,1,'READ_GRIB_TG_ECMWF',iltype,ilev1,ilev2,4,0.07,0.07,zd(1))
732  ALLOCATE(ztg(SIZE(zfield),5))
733  ztg(:,1)=zfield
734 ELSE
735  CALL abor1_sfx('MODE_READ_GRIB: SOIL TEMPERATURE LEVEL 1 MISSING (READ_GRIB_TG_ECMWF)')
736 ENDIF
737 !
738 ! 2. Search and read level 4 (and its depth) This level is optionnal
739 ! ---------------------------------------------------------------
740 iltype= -1
741 ilev1 = -1
742 ilev2 = -1
743  CALL read_grib(hgrib,kluout,236,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
744 !
745 IF (iret == 0) THEN
746  inlayerdeep = 4
747  CALL put_layer_depth(kluout,4,'READ_GRIB_TG_ECMWF',iltype,ilev1,ilev2,inlayerdeep,1.89,1.89,zd(4))
748  ztg(:,4)=zfield
749 ELSE
750  inlayerdeep = 3
751  zd(4) = 0.
752 ENDIF
753 !
754 ! 3. Search and read level 3 (and its depth) This level is optionnal
755 ! ---------------------------------------------------------------
756 iltype= -1
757 ilev1 = -1
758 ilev2 = -1
759  CALL read_grib(hgrib,kluout,183,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
760 !
761 IF (iret == 0) THEN
762  CALL put_layer_depth(kluout,3,'READ_GRIB_TG_ECMWF',iltype,ilev1,ilev2,inlayerdeep,0.72,0.42,zd(3))
763  ztg(:,3)=zfield
764 ELSE
765  inlayerdeep = 2
766  zd(3) = 0.
767 ENDIF
768 !
769 ! 4. Search and read level 2 (and its depth)
770 ! ---------------------------------------
771 iltype= -1
772 ilev1 = -1
773 ilev2 = -1
774  CALL read_grib(hgrib,kluout,170,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
775 !
776 IF (iret== 0) THEN
777  CALL put_layer_depth(kluout,2,'READ_GRIB_TG_ECMWF',iltype,ilev1,ilev2,inlayerdeep,0.21,0.42,zd(2))
778  ztg(:,2)=zfield
779  DEALLOCATE(zfield)
780 ELSE
781  CALL abor1_sfx('MODE_READ_GRIB: SOIL TEMPERATURE LEVEL 2 MISSING (READ_GRIB_TG_ECMWF)')
782 ENDIF
783 !--------------------------------------------------------------------------------
784 ! 5. Assumes uniform temperature profile up to 3m depth
785 ! -------------------------------------------------
786 !
787 IF(sum(zd(1:inlayerdeep)) < 3.) THEN
788  !We add a temperature layer
789  inlayerdeep=inlayerdeep+1
790  zd(inlayerdeep)=3.-sum(zd(1:inlayerdeep-1))
791  ztg(:,inlayerdeep)=ztg(:,inlayerdeep-1)
792 ENDIF
793 !
794 !--------------------------------------------------------------------------------
795 ! 6. Set temperature profile and layer thicknesses
796 ! ----------------------------------------------
797  CALL fill_pfield(kluout,'READ_GRIB_TG_ECMWF',inlayerdeep,zd,ztg,pmask,ptg,pd)
798 DEALLOCATE(zd)
799 DEALLOCATE(ztg)
800 !
801 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TG_ECMWF',1,zhook_handle)
802 END SUBROUTINE read_grib_tg_ecmwf
803 !-------------------------------------------------------------------
804 ! #######################
805  SUBROUTINE read_grib_wg_ecmwf_1(HGRIB,KLUOUT,HINMODEL,PMASK,PWG,PD)
806 ! #######################
807 !
808 ! This tasks is divided in the following steps :
809 ! - computing the MesoNH constants
810 ! - reading the grib datas according to the type of file (ECMWF/Arpege/Aladin)
811 ! - converting from specific humidity to relative humidity
812 ! - interpolation with land mask
813 ! - converting back from relative humidity to specific humidity with MesoNH constants
814 ! Five different models are supported :
815 ! - ECMWF with 2 layers (untested)
816 ! - ECMWF with 3 layers (archive before 1991 - Blondin model)
817 ! - ECMWF with 4 layers (archive after 1991 - Viterbo model)
818 ! - Arpege/Aladin before ISBA (I don't know the name of this model)
819 ! - Arpege/Aladin with ISBA model
820 ! The available model is detect according to the fields which are presents :
821 ! - ECMWF archive : loads as many layers as possible
822 ! - Arpege/Aladin archive : ISBA model needs Clay and Sans fraction fields, if they
823 ! are present, they are used and the model is declared to be ISBA.
824 ! To detect the height of the layers, two methods are used :
825 ! - if level type is not 112, a default value is assumed and a warning message is
826 ! displayed
827 ! - if level type is ID 112, then the position of the top and bottom surface may be
828 ! given. If they are present, they are used, if not the default value is taken and
829 ! a warning message is issued.
830 !
831 IMPLICIT NONE
832 !
833 !* dummy arguments
834 ! ---------------
835  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
836 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
837  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
838 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
839 REAL, DIMENSION(:,:), POINTER :: PWG ! field to initialize
840 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
841 !
842 !* local variables
843 ! ---------------
844 INTEGER(KIND=kindOfInt) :: IRET ! return code
845 INTEGER :: IPAR ! parameter number for field reading
846 INTEGER :: ILTYPE ! type of level (Grib code table 3)
847 INTEGER :: ILEV1 ! level definition
848 INTEGER :: ILEV2 ! level definition
849 INTEGER :: INLAYERDEEP! number of deep moisture layers
850 REAL, DIMENSION(:), POINTER :: ZFIELD => null() ! first layer temperature
851 REAL, DIMENSION(:,:), ALLOCATABLE:: ZWG ! first layer temperature
852 REAL, DIMENSION(:) , ALLOCATABLE:: ZD
853 REAL(KIND=JPRB) :: ZHOOK_HANDLE
854 !--------------------------------------------------------------------------------
855 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WG_ECMWF_1',0,zhook_handle)
856 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_WG_ECMWF_1: | Reading soil moisture'
857 !
858 ALLOCATE(zd(4))
859 !
860 ! 1. Search and read level 1 (and its depth)
861 ! --------------------------------------
862 iltype= -1
863 ilev1 = -1
864 ilev2 = -1
865 ipar=39
866  CALL read_grib(hgrib,kluout,140,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2,kparam2=ipar)
867 !
868 IF (iret == 0) THEN
869  CALL put_layer_depth(kluout,1,'READ_GRIB_WG_ECMWF_1',iltype,ilev1,ilev2,4,0.07,0.07,zd(1))
870  ALLOCATE(zwg(SIZE(zfield,1),4))
871  zwg(:,1)=zfield
872  !
873  IF (ipar==140) zwg(:,1)=zwg(:,1) / zd(1)
874 ELSE
875  CALL abor1_sfx('MODE_READ_GRIB: SOIL MOISTURE LEVEL 1 MISSING (READ_GRIB_WG_ECMWF_1)')
876 ENDIF
877 !
878 ! 2. Search and read level 4 (and its depth) This level is optionnal
879 ! ---------------------------------------------------------------
880 iltype= -1
881 ilev1 = -1
882 ilev2 = -1
883 ipar=42
884  CALL read_grib(hgrib,kluout,237,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2,kparam2=ipar)
885 !
886 IF (iret == 0) THEN
887  inlayerdeep = 4
888  CALL put_layer_depth(kluout,4,'READ_GRIB_WG_ECMWF_1',iltype,ilev1,ilev2,inlayerdeep,1.89,1.89,zd(4))
889  zwg(:,4)=zfield
890  !
891  IF (ipar==237) zwg(:,4)=zwg(:,4) / zd(1)
892 ELSE
893  inlayerdeep = 3
894  zd(4) = 0.
895 ENDIF
896 !
897 ! 3. Search and read level 3 (and its depth) This level is optionnal
898 ! ---------------------------------------------------------------
899 iltype= -1
900 ilev1 = -1
901 ilev2 = -1
902 ipar=41
903  CALL read_grib(hgrib,kluout,184,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2,kparam2=ipar)
904 !
905 IF (iret == 0) THEN
906  CALL put_layer_depth(kluout,3,'READ_GRIB_WG_ECMWF_1',iltype,ilev1,ilev2,inlayerdeep,0.72,0.42,zd(3))
907  zwg(:,3)=zfield
908  !
909  IF (ipar==184) zwg(:,3)=zwg(:,3) / zd(1)
910 ELSE
911  inlayerdeep = 2
912  zd(3) = 0.
913 ENDIF
914 !
915 ! 4. Search and read level 2 (and its depth)
916 ! ---------------------------------------
917 iltype= -1
918 ilev1 = -1
919 ilev2 = -1
920 ipar=40
921  CALL read_grib(hgrib,kluout,171,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2,kparam2=ipar)
922 !
923 IF (iret == 0) THEN
924  CALL put_layer_depth(kluout,2,'READ_GRIB_WG_ECMWF_1',iltype,ilev1,ilev2,inlayerdeep,0.21,0.42,zd(2))
925  zwg(:,2)=zfield
926  DEALLOCATE(zfield)
927  !
928  IF (ipar==171) zwg(:,2)=zwg(:,2) / zd(1)
929 ELSE
930  CALL abor1_sfx('MODE_READ_GRIB: SOIL MOISTURE LEVEL 2 MISSING (READ_GRIB_WG_ECMWF_1)')
931 ENDIF
932 !
933 !--------------------------------------------------------------------------------
934 !
935  CALL fill_pfield(kluout,'READ_GRIB_WG_ECMWF_1',inlayerdeep,zd,zwg,pmask,pwg,pd)
936 DEALLOCATE(zd)
937 DEALLOCATE(zwg)
938 !
939 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WG_ECMWF_1',1,zhook_handle)
940 END SUBROUTINE read_grib_wg_ecmwf_1
941 !-------------------------------------------------------------------
942 ! #######################
943  SUBROUTINE ecmwf_wgi(PTG,PWG,PWGI)
944 ! #######################
945 !
946 ! ECMWF grib only contain (ice+water) content.
947 ! This routine computes iced part and water part according to the formula
948 ! given in ECMWF documentation. But we use real water content instead of
949 ! (CL+CH) times saturation capacity.
950 !
951 USE modd_csts, ONLY : xtt, xpi
952 !
953 IMPLICIT NONE
954 !
955 !* dummy arguments
956 ! ---------------
957 REAL, DIMENSION(:,:), INTENT(IN) :: PTG ! Temperature profil
958 REAL, DIMENSION(:,:), INTENT(INOUT) :: PWG ! INPUT contains (water+ice) profil
959  ! OUTPUT contains only water profil
960 REAL, DIMENSION(:,:), INTENT(OUT) :: PWGI ! ice profil
961 !
962 !* local variables
963 ! ---------------
964 REAL :: ZT1, ZT2 ! Temperature threshold
965 REAL(KIND=JPRB) :: ZHOOK_HANDLE
966 !--------------------------------------------------------------------------------
967 IF (lhook) CALL dr_hook('MODE_READ_GRIB:ECMWF_WGI',0,zhook_handle)
968 !
969 zt1=xtt + 1.
970 zt2=xtt - 3.
971 !
972 WHERE(ptg(:,:) > zt1)
973  pwgi(:,:) = 0.
974 ELSEWHERE(ptg(:,:) < zt2)
975  pwgi(:,:) = pwg(:,:)
976  pwg(:,:) = 0.
977 ELSEWHERE
978  pwgi(:,:)=pwg(:,:) * 0.5* (1 - sin(xpi * (ptg(:,:) - 0.5*zt1 - 0.5*zt2) / &
979  (zt1 - zt2 ) ))
980  pwg(:,:) = pwg(:,:) - pwgi(:,:)
981 ENDWHERE
982 !
983 IF (lhook) CALL dr_hook('MODE_READ_GRIB:ECMWF_WGI',1,zhook_handle)
984 END SUBROUTINE ecmwf_wgi
985 !--------------------------------------------------------------------------------
986 ! #######################
987  SUBROUTINE harmonize_grib_wg_wgi_ecmwf(HGRIB,KLUOUT,HINMODEL,PMASK,PWG,PD,PWGI)
988 ! #######################
989 !
990 ! ECMWF grib only contain (ice+water) content.
991 ! This routine computes iced part and water part according to the formula
992 ! given in ECMWF documentation. But we use real water content instead of
993 ! (CL+CH) times saturation capacity.
994 !
995 USE modd_csts, ONLY : xtt, xpi
996 !
997 IMPLICIT NONE
998 !
999 !* dummy arguments
1000 ! ---------------
1001  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1002 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1003  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1004 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1005 REAL, DIMENSION(:,:), OPTIONAL, POINTER :: PWG ! INPUT contains (water+ice) profil
1006  ! OUTPUT contains only water profil
1007 REAL, DIMENSION(:,:), OPTIONAL, POINTER :: PD ! thickness of each layer
1008 REAL, DIMENSION(:,:), OPTIONAL, POINTER :: PWGI ! ice profil
1009 !* local variables
1010 ! ---------------
1011 REAL, DIMENSION(:,:), POINTER :: ZWG => null() ! profile of soil water contents
1012 REAL, DIMENSION(:,:), POINTER :: ZD => null() ! thickness of each layer
1013 REAL, DIMENSION(:,:), POINTER :: ZTG => null() ! profile of temperature
1014 REAL, DIMENSION(:,:), POINTER :: ZDT => null() ! thickness of each temperature layer
1015 REAL, DIMENSION(:,:), ALLOCATABLE:: ZWGI ! profile of soil ice contents
1016 REAL :: ZT1, ZT2 ! Temperature threshold
1017 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1018 !--------------------------------------------------------------------------------
1019 IF (lhook) CALL dr_hook('MODE_READ_GRIB:HARMONIZE_GRIB_WG_WGI_ECMWF',0,zhook_handle)
1020 !
1021  CALL read_grib_tg_ecmwf(hgrib,kluout,hinmodel,pmask,ztg,zdt)
1022  CALL read_grib_wg_ecmwf_1(hgrib,kluout,hinmodel,pmask,zwg,zd)
1023 !
1024 IF (SIZE(ztg,2) .LT. SIZE(zwg,2)) THEN
1025  WRITE (kluout,'(A)') 'MODE_READ_GRIB:HARMONIZE_GRIB_WG_WGI_ECMWF: '
1026  WRITE (kluout,'(A)') 'ERROR, YOU HAVE NOT THE SAME NUMBER OF LEVELS '
1027  WRITE (kluout,'(A)') 'IN SOIL FOR TEMPERATURE AND HUMIDITY '
1028  WRITE (kluout,'(A)') 'VERIFY GRIB FILE '
1029  CALL abor1_sfx("MODE_READ_GRIB:HARMONIZE_GRIB_WG_WGI_ECMWF: VERIFY NUMBER OF LEVELS IN GRIB FILE")
1030 ENDIF
1031 !
1032 IF (PRESENT(pd)) THEN
1033  ALLOCATE(pd(SIZE(zd,1),SIZE(zd,2)))
1034  pd(:,:)=zd(:,:)
1035 ENDIF
1036 IF (PRESENT(pwgi)) THEN
1037  ALLOCATE(pwgi(SIZE(zwg,1),SIZE(zwg,2)))
1038  pwgi(:,:)=0.
1039 ENDIF
1040 !
1041 !If same vertical grids are taken into account for WG and TG we can
1042 !compute ice content and new water content
1043 IF(all(zdt(:,1:SIZE(zwg,2))==zd(:,1:SIZE(zwg,2)))) THEN
1044  ALLOCATE(zwgi(SIZE(zwg,1),SIZE(zwg,2)))
1045  CALL ecmwf_wgi(ztg(:,1:SIZE(zwg,2)),zwg,zwgi)
1046  IF (PRESENT(pwgi)) pwgi(:,:)=zwgi(:,:)
1047  DEALLOCATE(zwgi)
1048 ENDIF
1049 !
1050 IF (PRESENT(pwg)) THEN
1051  ALLOCATE(pwg(SIZE(zwg,1),SIZE(zwg,2)))
1052  pwg(:,:)=zwg(:,:)
1053 ENDIF
1054 !
1055 DEALLOCATE(zwg)
1056 DEALLOCATE(zd)
1057 DEALLOCATE(ztg)
1058 DEALLOCATE(zdt)
1059 !
1060 IF (lhook) CALL dr_hook('MODE_READ_GRIB:HARMONIZE_GRIB_WG_WGI_ECMWF',1,zhook_handle)
1061 END SUBROUTINE harmonize_grib_wg_wgi_ecmwf
1062 !--------------------------------------------------------------------------------
1063 ! #######################
1064  SUBROUTINE read_grib_wg_ecmwf(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
1065 ! #######################
1066 !
1067 USE modd_grid_grib, ONLY : nni
1068 USE modd_surf_par, ONLY : xundef
1069 !
1070 IMPLICIT NONE
1071 !
1072 !* dummy arguments
1073 ! ---------------
1074  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1075 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1076  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1077 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1078 REAL, DIMENSION(:,:), POINTER :: PFIELD ! field to initialize
1079 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
1080 !* local variables
1081 ! ---------------
1082 INTEGER(KIND=kindOfInt) :: IRET ! return code
1083 REAL, DIMENSION(:), POINTER :: ZSLT => null() ! soil type
1084 REAL, DIMENSION(:), ALLOCATABLE:: ZWWILT ! ECMWF wilting point
1085 REAL, DIMENSION(:), ALLOCATABLE:: ZWFC ! ECMWF field capacity
1086 INTEGER :: JL ! loop counter on layers
1087 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1088 !
1089 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WG_ECMWF',0,zhook_handle)
1090 !
1091  CALL harmonize_grib_wg_wgi_ecmwf(hgrib,kluout,hinmodel,pmask,pwg=pfield,pd=pd)
1092 !
1093 ! 1. Get soil type to compute SWI
1094 ! ----------------------------
1095  CALL read_grib(hgrib,kluout,43,iret,zslt)
1096 !--------------------------------------------------------------------------------
1097 ALLOCATE (zwfc(SIZE(pfield,1)))
1098 ALLOCATE (zwwilt(SIZE(pfield,1)))
1099 zwfc(:) = 0.
1100 zwwilt(:) = 0.
1101 !
1102 IF (iret == 0) THEN
1103 !
1104 ! 2.1 Convert from specific humidity to relative humidity using soil types
1105 ! --------------------------------------------------------------------
1106  WHERE (zslt(:)==1.)
1107  zwfc(:) = 0.242
1108  zwwilt(:) = 0.059
1109  ELSEWHERE (zslt(:)==2.)
1110  zwfc(:) = 0.346
1111  zwwilt(:) = 0.151
1112  ELSEWHERE (zslt(:)==3.)
1113  zwfc(:) = 0.382
1114  zwwilt(:) = 0.133
1115  ELSEWHERE (zslt(:)==4.)
1116  zwfc(:) = 0.448
1117  zwwilt(:) = 0.279
1118  ELSEWHERE (zslt(:)==5.)
1119  zwfc(:) = 0.541
1120  zwwilt(:) = 0.335
1121  ELSEWHERE (zslt(:)==6.)
1122  zwfc(:) = 0.662
1123  zwwilt(:) = 0.267
1124  ENDWHERE
1125  !
1126 ELSE
1127 !
1128 ! 2.2 Convert from specific humidity to relative humidity single soil type
1129 ! --------------------------------------------------------------------
1130  ! Compute model's constants
1131  IF (SIZE(pfield,2)==4) THEN
1132  zwfc(:) = 0.323
1133  zwwilt(:) = 0.171
1134  ELSE
1135  zwfc(:) = 0.171
1136  zwwilt(:) = 0.086
1137  END IF
1138  !
1139 ENDIF
1140 !
1141 DO jl=1,SIZE(pfield,2)
1142  WHERE ( pfield(:,jl).NE.xundef .AND. zwfc(:).NE.0. )
1143  pfield(:,jl) = (pfield(:,jl) - zwwilt(:)) / (zwfc(:) - zwwilt(:))
1144  ELSEWHERE
1145  pfield(:,jl) = 0.
1146  ENDWHERE
1147 ENDDO
1148 !
1149 IF (ASSOCIATED(zslt)) DEALLOCATE(zslt)
1150 DEALLOCATE(zwfc)
1151 DEALLOCATE(zwwilt)
1152 !
1153 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WG_ECMWF',1,zhook_handle)
1154 END SUBROUTINE read_grib_wg_ecmwf
1155 !----------------------------------------------------------------------------
1156 ! #######################
1157  SUBROUTINE read_grib_wgi_ecmwf(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
1158 ! #######################
1159 !
1160 USE modd_grid_grib, ONLY : nni
1161 USE modd_surf_par, ONLY : xundef
1162 !
1163 IMPLICIT NONE
1164 !
1165 !* dummy arguments
1166 ! ---------------
1167  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1168 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1169  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1170 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1171 REAL, DIMENSION(:,:), POINTER :: PFIELD ! field to initialize
1172 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
1173 !* local variables
1174 ! ---------------
1175 INTEGER(KIND=kindOfInt) :: IRET ! return code
1176 REAL, DIMENSION(:), POINTER :: ZSLT => null() ! soil type
1177 REAL, DIMENSION(:) , ALLOCATABLE:: ZWSAT ! ECMWF saturation
1178 INTEGER :: JL ! loop counter on layers
1179 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1180 !--------------------------------------------------------------------------------
1181 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WGI_ECMWF',0,zhook_handle)
1182 !
1183  CALL harmonize_grib_wg_wgi_ecmwf(hgrib,kluout,hinmodel,pmask,pd=pd,pwgi=pfield)
1184 !
1185 ! 1. Get soil type to compute WSAT
1186 ! ----------------------------
1187  CALL read_grib(hgrib,kluout,43,iret,zslt)
1188 !--------------------------------------------------------------------------------
1189 ALLOCATE (zwsat(SIZE(pfield,1)))
1190 zwsat(:)=0.
1191 !
1192 IF (iret == 0) THEN
1193 !
1194 ! 2.1 Convert from specific humidity to relative humidity using soil types
1195 ! --------------------------------------------------------------------
1196  WHERE (zslt(:)==1.)
1197  zwsat(:) = 0.403
1198  ELSEWHERE (zslt(:)==2.)
1199  zwsat(:) = 0.439
1200  ELSEWHERE (zslt(:)==3.)
1201  zwsat(:) = 0.430
1202  ELSEWHERE (zslt(:)==4.)
1203  zwsat(:) = 0.520
1204  ELSEWHERE (zslt(:)==5.)
1205  zwsat(:) = 0.614
1206  ELSEWHERE (zslt(:)==6.)
1207  zwsat(:) = 0.766
1208  ENDWHERE
1209 !
1210 ELSE
1211 !
1212 ! 2.2 Convert from specific humidity to relative humidity single soil type
1213 ! --------------------------------------------------------------------
1214  ! Compute model's constants
1215  IF (SIZE(pfield,2)==4) THEN
1216  zwsat(:) = 0.472
1217  ELSE
1218  zwsat(:) = 0.286
1219  END IF
1220  !
1221 ENDIF
1222 !
1223 ! Then perform conversion
1224 DO jl=1,SIZE(pfield,2)
1225  WHERE ( pfield(:,jl).NE.xundef .AND. zwsat(:).NE.0. )
1226  pfield(:,jl) = pfield(:,jl) / zwsat(:)
1227  ELSEWHERE
1228  pfield(:,jl) = 0.
1229  ENDWHERE
1230 ENDDO
1231 !
1232 IF (ASSOCIATED(zslt)) DEALLOCATE(zslt)
1233 DEALLOCATE(zwsat)
1234 !
1235 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WGI_ECMWF',1,zhook_handle)
1236 END SUBROUTINE read_grib_wgi_ecmwf
1237 !-------------------------------------------------------------------
1238 !-------------------------------------------------------------------
1239 ! #######################
1240  SUBROUTINE read_grib_tg_meteo_france(HGRIB,KLUOUT,HINMODEL,PMASK,PTG,PDT)
1241 ! #######################
1242 !
1243 USE modd_surf_par, ONLY : xundef
1244 !
1245 IMPLICIT NONE
1246 !
1247 !* dummy arguments
1248 ! ---------------
1249  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1250 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1251  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1252 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1253 REAL, DIMENSION(:,:), POINTER :: PTG ! field to initialize
1254 REAL, DIMENSION(:,:), POINTER :: PDT ! thickness of each layer
1255 !* local variables
1256 ! ---------------
1257 REAL, DIMENSION(:), POINTER :: ZFIELD => null() ! field to read
1258 INTEGER :: JL ! layer loop counter
1259 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1260 !--------------------------------------------------------------------------------
1261 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TG_METEO_FRANCE',0,zhook_handle)
1262 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_TG_METEO_FRANCE: | Reading soil temperature'
1263 !--------------------------------------------------------------------------------
1264 ! 1. Allocate soil temperature profile
1265 ! ---------------------------------
1266 !--------------------------------------------------------------------------------
1267 ! 2. Search and read level 1 (and its depth)
1268 ! ---------------------------------------
1269  CALL read_grib_ts(hgrib,kluout,hinmodel,pmask,zfield)
1270 !
1271 ALLOCATE(ptg(SIZE(zfield),3))
1272 ALLOCATE(pdt(SIZE(zfield),3))
1273 !
1274 ptg(:,1) = zfield(:)
1275 pdt(:,1) = 0.01
1276 !--------------------------------------------------------------------------------
1277 ! 3. Deep soil temperature
1278 ! ---------------------
1279  CALL read_grib_t2_land(hgrib,kluout,hinmodel,pmask,zfield)
1280 !
1281 ptg(:,2) = zfield(:)
1282 pdt(:,2) = 0.4 ! deep temperature layer depth assumed equal to 0.4m
1283 DEALLOCATE(zfield)
1284 !--------------------------------------------------------------------------------
1285 ! 4. Assumes uniform temperature profile below
1286 ! -----------------------------------------
1287 ptg(:,3) = ptg(:,2)
1288 pdt(:,3) = 5. ! temperature profile down to 5m
1289 !
1290 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TG_METEO_FRANCE',1,zhook_handle)
1291 END SUBROUTINE read_grib_tg_meteo_france
1292 !-------------------------------------------------------------------
1293 ! #######################
1294  SUBROUTINE read_grib_sand_clay_meteo_france(HGRIB,KLUOUT,HINMODEL,PSAND,PCLAY,GISBA)
1295 ! ######################
1296 !
1297 IMPLICIT NONE
1298 !
1299 !* dummy arguments
1300 ! ---------------
1301  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1302 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1303  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1304 REAL, DIMENSION(:), POINTER :: PSAND ! field to initialize
1305 REAL, DIMENSION(:), POINTER :: PCLAY ! thickness of each layer
1306 LOGICAL, INTENT(OUT) :: GISBA ! T: surface scheme in file is ISBA
1307 !* local variables
1308 ! ---------------
1309 INTEGER(KIND=kindOfInt) :: IRET ! return code
1310 INTEGER :: IPAR ! parameter number for field reading
1311 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1312 !-------------------------------------------------------------------------------
1313 ! 1. Search and read clay fraction if available
1314 ! ------------------------------------------
1315 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SAND_CLAY_METEO_FRANCE',0,zhook_handle)
1316 !
1317 IF (hinmodel == 'ARPEGE' .OR. hinmodel == 'MOCAGE') ipar=171
1318 IF (hinmodel == 'ALADIN') ipar=128
1319  CALL read_grib(hgrib,kluout,ipar,iret,pclay)
1320 !
1321 ! if not available, the model is not ISBA (IWMODE=1)
1322 IF (iret /= 0) THEN
1323  gisba = .false.
1324 ELSE
1325  gisba = .true.
1326  pclay(:) = pclay(:) / 100. ! this field is given in percent
1327  WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_SAND_CLAY_METEO_FRANCE: | The soil model is ISBA'
1328 END IF
1329 !-------------------------------------------------------------------------------
1330 ! 2. Search and read sand fraction if available
1331 ! ------------------------------------------
1332 IF (hinmodel == 'ARPEGE' .OR. hinmodel == 'MOCAGE') ipar=172
1333 IF (hinmodel == 'ALADIN') ipar=129
1334  CALL read_grib(hgrib,kluout,ipar,iret,psand)
1335 !
1336 ! if not available, the model is not ISBA (IWMODE=1)
1337 IF (gisba) THEN
1338  IF (iret /= 0) THEN
1339  CALL abor1_sfx('MODE_READ_GRIB: SAND FRACTION MISSING (READ_GRIB_SAND_CLAY_METEO_FRANCE)')
1340  ELSE
1341  psand(:) = psand(:) / 100. ! this field is given in percent
1342  END IF
1343 END IF
1344 !
1345 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SAND_CLAY_METEO_FRANCE',1,zhook_handle)
1346 END SUBROUTINE read_grib_sand_clay_meteo_france
1347 !-----------------------------------------------------------------------------------
1348 ! #######################
1349  SUBROUTINE read_grib_wg_meteo_france(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
1350 ! #######################
1351 !
1352 USE modd_grid_grib, ONLY : nni
1353 USE modd_surf_par, ONLY : xundef
1354 !
1355 IMPLICIT NONE
1356 !
1357 !* dummy arguments
1358 ! ---------------
1359  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1360 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1361  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1362 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1363 REAL, DIMENSION(:,:), POINTER :: PFIELD ! field to initialize
1364 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
1365 !
1366 !* local variables
1367 ! ---------------
1368 LOGICAL :: GISBA ! T: surface scheme in file is ISBA
1369 INTEGER(KIND=kindOfInt) :: IRET ! return code
1370 INTEGER :: ILTYPE ! type of level (Grib code table 3)
1371 INTEGER :: ILEV1 ! level definition
1372 INTEGER :: ILEV2 ! level definition
1373 REAL, DIMENSION(:), POINTER :: ZCLAY => null() ! clay fraction
1374 REAL, DIMENSION(:), POINTER :: ZSAND => null() ! sand fraction
1375 REAL, DIMENSION(:), POINTER :: ZFIELD => null()
1376 REAL, DIMENSION(:), ALLOCATABLE:: ZWWILT ! wilting point
1377 REAL, DIMENSION(:), ALLOCATABLE:: ZWFC ! field capacity
1378 REAL, DIMENSION(:), ALLOCATABLE:: ZWSAT ! saturation
1379 INTEGER :: JL ! layer loop counter
1380 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1381 !-------------------------------------------------------------------------------
1382 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WG_METEO_FRANCE',0,zhook_handle)
1383 !
1384 ! 1. Search and read clay and sand fractions if available
1385 ! ----------------------------------------------------
1386  CALL read_grib_sand_clay_meteo_france(hgrib,kluout,hinmodel,zsand,zclay,gisba)
1387 !-------------------------------------------------------------------------------
1388 IF (gisba) THEN
1389  ALLOCATE(pfield(SIZE(zsand),3))
1390  ALLOCATE(pd(SIZE(zsand),3))
1391 ELSE
1392  ALLOCATE(pfield(nni,3))
1393  ALLOCATE(pd(nni,3))
1394 ENDIF
1395 !
1396 pd(:,1) = 0.01
1397 pd(:,2) = 0.20
1398 !-------------------------------------------------------------------------------
1399 ! 2. Read layer 1 moisture
1400 ! ---------------------
1401 ilev1 = 0
1402 IF (hinmodel == 'ARPEGE' .OR. hinmodel=='MOCAGE') THEN
1403  iltype = 112
1404  ilev2 = 1
1405  CALL read_grib(hgrib,kluout,153,iret,zfield,klev1=ilev1,klev2=ilev2)
1406 ELSE
1407  iltype = 105
1408  ilev2 = 0
1409  CALL read_grib(hgrib,kluout,86,iret,zfield,klev1=ilev1,klev2=ilev2)
1410 ENDIF
1411 IF (iret /= 0) THEN
1412  CALL abor1_sfx('MODE_READ_GRIB: SOIL MOISTURE LEVEL 1 MISSING (READ_GRIB_WG_METEO_FRANCE)')
1413 END IF
1414 !
1415 pfield(:,1) = zfield(:)
1416 !-------------------------------------------------------------------------------
1417 ! 3. Read layer 2 moisture
1418 ! ---------------------
1419 IF (hinmodel == 'ARPEGE' .OR. hinmodel=='MOCAGE') THEN
1420  iltype = 112
1421  ilev1 = 0
1422  ilev2 = 250
1423  CALL read_grib(hgrib,kluout,153,iret,zfield,klev1=ilev1,klev2=ilev2)
1424 ELSE
1425  iltype = 111
1426  ilev1 = -1
1427  ilev2 = -1
1428  CALL read_grib(hgrib,kluout,86,iret,zfield,klev1=ilev1,klev2=ilev2)
1429 ENDIF
1430 IF (iret /= 0) THEN
1431  CALL abor1_sfx('MODE_READ_GRIB: SOIL MOISTURE LEVEL 2 MISSING (READ_GRIB_WG_METEO_FRANCE)')
1432 END IF
1433 !
1434 pfield(:,2) = zfield(:)
1435 !-------------------------------------------------------------------------------
1436 ! 4. Read layer 2 depth (ISBA only)
1437 ! -----------------------------
1438 !* note that soil water reservoir is considered uniform between 0.2m and GRIB soil depth
1439 IF (gisba) THEN
1440  IF (hinmodel == 'ARPEGE' .OR. hinmodel == 'MOCAGE') THEN
1441  CALL read_grib(hgrib,kluout,173,iret,zfield)
1442  ELSE
1443  CALL read_grib(hgrib,kluout,130,iret,zfield)
1444  ENDIF
1445  IF (iret /= 0) THEN
1446  CALL abor1_sfx('MODE_READ_GRIB: SOIL MOISTURE LEVEL 2 DEPTH MISSING (READ_GRIB_WG_METEO_FRANCE)')
1447  END IF
1448  pd(:,3) = zfield(:)
1449  DEALLOCATE(zfield)
1450 ELSE
1451  pd(:,3) = 2.
1452 END IF
1453 !-------------------------------------------------------------------------------
1454 ! 5. Compute relative humidity from units kg/m^2
1455 ! -------------------------------------------
1456 ! Compute ISBA model constants (if needed)
1457 IF (gisba) THEN
1458  !
1459  !* updates Wg in m3/m3
1460  pfield(:,1) = pfield(:,1) / 10.
1461  pfield(:,2) = pfield(:,2) /(1000. * pd(:,3))
1462  !
1463  ALLOCATE (zwsat(SIZE(zsand)))
1464  zwsat(:) = (-1.08*100. * zsand(:) + 494.305) * 0.001
1465  pfield(:,1) = max(min(pfield(:,1),zwsat(:)),0.)
1466  pfield(:,2) = max(min(pfield(:,2),zwsat(:)),0.)
1467  DEALLOCATE(zwsat)
1468  DEALLOCATE (zsand)
1469  !
1470  ALLOCATE (zwwilt(SIZE(zclay)))
1471  ALLOCATE (zwfc(SIZE(zclay)))
1472  zwwilt(:) = 37.1342e-3 * sqrt( 100. * zclay(:) )
1473  zwfc(:) = 89.0467e-3 * (100. * zclay(:) )**0.3496
1474  pfield(:,1) = (pfield(:,1) - zwwilt(:)) / (zwfc(:) - zwwilt(:))
1475  pfield(:,2) = (pfield(:,2) - zwwilt(:)) / (zwfc(:) - zwwilt(:))
1476  DEALLOCATE (zwwilt)
1477  DEALLOCATE (zwfc)
1478  DEALLOCATE (zclay)
1479  !
1480 ELSE ! Non ISBA
1481  !
1482  pfield(:,2) = (pfield(:,1)+pfield(:,2)) / (20. + 100.)
1483  pfield(:,1) = pfield(:,1) / 20.
1484  !
1485 END IF
1486 !
1487 pfield(:,3) = pfield(:,2)
1488 !--------------------------------------------------------------------------------
1489 ! 6. Apply land mask
1490 ! ---------------
1491 IF (SIZE(pmask)==SIZE(pfield,1)) THEN
1492  DO jl=1,SIZE(pfield,2)
1493  WHERE (pmask(:)/=1.) pfield(:,jl) = xundef
1494  END DO
1495 ENDIF
1496 !
1497 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WG_METEO_FRANCE',1,zhook_handle)
1498 END SUBROUTINE read_grib_wg_meteo_france
1499 !--------------------------------------------------------------------------------
1500 ! #######################
1501  SUBROUTINE read_grib_wgi_meteo_france(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
1502 ! #######################
1503 !
1504 USE modd_grid_grib, ONLY : nni
1505 USE modd_surf_par, ONLY : xundef
1506 !
1507 IMPLICIT NONE
1508 !
1509 !* dummy arguments
1510 ! ---------------
1511  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1512 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1513  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1514 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1515 REAL, DIMENSION(:,:), POINTER :: PFIELD ! field to initialize
1516 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
1517 !
1518 !* local variables
1519 ! ---------------
1520 LOGICAL :: GISBA ! T: surface scheme in file is ISBA
1521 INTEGER(KIND=kindOfInt) :: IRET ! return code
1522 INTEGER :: ILTYPE ! type of level (Grib code table 3)
1523 INTEGER :: ILEV1 ! level definition
1524 INTEGER :: ILEV2 ! level definition
1525 REAL, DIMENSION(:), POINTER :: ZCLAY => null() ! clay fraction
1526 REAL, DIMENSION(:), POINTER :: ZSAND => null() ! sand fraction
1527 REAL, DIMENSION(:), POINTER :: ZFIELD => null()
1528 REAL, DIMENSION(:), ALLOCATABLE:: ZWSAT ! saturation
1529 INTEGER :: JL ! layer loop counter
1530 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1531 !-------------------------------------------------------------------------------
1532 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WGI_METEO_FRANCE',0,zhook_handle)
1533 !
1534 ! 1. Search and read clay fraction if available
1535 ! ------------------------------------------
1536  CALL read_grib_sand_clay_meteo_france(hgrib,kluout,hinmodel,zsand,zclay,gisba)
1537 !-------------------------------------------------------------------------------
1538 IF (gisba) THEN
1539  ALLOCATE(pfield(SIZE(zsand),2))
1540  ALLOCATE(pd(SIZE(zsand),2))
1541 ELSE
1542  ALLOCATE(pfield(nni,2))
1543  ALLOCATE(pd(nni,2))
1544 ENDIF
1545 !
1546 pd(:,1) = 0.01
1547 !-------------------------------------------------------------------------------
1548 ! 2. Read layer 1 soil ice
1549 ! ---------------------
1550 ilev1 = 0
1551 IF (hinmodel == 'ARPEGE' .OR. hinmodel=='MOCAGE') THEN
1552  iltype = 112
1553  ilev2 = 1
1554  CALL read_grib(hgrib,kluout,152,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
1555 ELSE
1556  iltype = 105
1557  ilev2 = 0
1558  CALL read_grib(hgrib,kluout,139,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
1559 END IF
1560 !
1561 IF (iret == 0) THEN
1562  pfield(:,1) = zfield(:)
1563  WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_WGI_METEO_FRANCE: -> Soil ice level 1 is present'
1564 ELSE
1565  pfield(:,1) = 0.
1566 END IF
1567 !-------------------------------------------------------------------------------
1568 ! 3. Read layer 2 soil ice
1569 ! ---------------------
1570 IF (hinmodel == 'ARPEGE' .OR. hinmodel=='MOCAGE') THEN
1571  iltype = 112
1572  ilev1 = 0
1573  ilev2 = 250
1574  CALL read_grib(hgrib,kluout,152,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
1575 ELSE
1576  iltype = 111
1577  ilev1 = -1
1578  ilev2 = -1
1579  CALL read_grib(hgrib,kluout,139,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
1580 END IF
1581 !
1582 IF (iret == 0) THEN
1583  pfield(:,2) = zfield(:)
1584  WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_WGI_METEO_FRANCE: -> Soil ice level 2 is present'
1585 ELSE
1586  pfield(:,2) = 0.
1587 END IF
1588 !-------------------------------------------------------------------------------
1589 ! 4. Read layer 2 depth (ISBA only)
1590 ! ------------------------------
1591 IF (gisba) THEN
1592  IF (hinmodel == 'ARPEGE' .OR. hinmodel=='MOCAGE') THEN
1593  CALL read_grib(hgrib,kluout,173,iret,zfield)
1594  ELSE
1595  CALL read_grib(hgrib,kluout,130,iret,zfield)
1596  ENDIF
1597  IF (iret /= 0) THEN
1598  CALL abor1_sfx('MODE_READ_GRIB: SOIL ICE LEVEL 2 MISSING (READ_GRIB_WGI_METEO_FRANCE)')
1599  END IF
1600  pd(:,2) = zfield(:)
1601  DEALLOCATE(zfield)
1602 ELSE
1603  pd(:,2) = 2.
1604 END IF
1605 !-------------------------------------------------------------------------------
1606 ! 5. Compute relative humidity from units kg/m^2
1607 ! -------------------------------------------
1608 IF (gisba) THEN
1609  !
1610  !* updates Wgi in m3/m3
1611  pfield(:,1) = pfield(:,1) / 10.
1612  pfield(:,2) = pfield(:,2) /(1000. * pd(:,2))
1613  !
1614  ALLOCATE (zwsat(nni))
1615  zwsat(:) = (-1.08*100. * zsand(:) + 494.305) * 0.001
1616  pfield(:,1) = pfield(:,1) / zwsat(:)
1617  pfield(:,2) = pfield(:,2) / zwsat(:)
1618  DEALLOCATE (zwsat)
1619  DEALLOCATE (zsand)
1620  DEALLOCATE (zclay)
1621  !
1622 ELSE ! Non ISBA
1623  !
1624  pfield(:,1) = 0.
1625  pfield(:,2) = 0.
1626  !
1627 END IF
1628 !--------------------------------------------------------------------------------
1629 ! 6. Apply land mask
1630 ! ---------------
1631 IF (SIZE(pmask)==SIZE(pfield,1)) THEN
1632  DO jl=1,SIZE(pfield,2)
1633  WHERE (pmask(:)/=1.) pfield(:,jl) = xundef
1634  END DO
1635 ENDIF
1636 !
1637 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WGI_METEO_FRANCE',1,zhook_handle)
1638 END SUBROUTINE read_grib_wgi_meteo_france
1639 !-------------------------------------------------------------------
1640 !-------------------------------------------------------------------
1641 ! #######################
1642  SUBROUTINE read_grib_tg_hirlam(HGRIB,KLUOUT,HINMODEL,PMASK,PTG,PDT)
1643 ! #######################
1644 !
1645 USE modd_grid_grib, ONLY : nni
1646 USE modd_surf_par, ONLY : xundef
1647 !
1648 IMPLICIT NONE
1649 !
1650 !* dummy arguments
1651 ! ---------------
1652  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1653 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1654  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1655 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1656 REAL, DIMENSION(:,:), POINTER :: PTG ! field to initialize
1657 REAL, DIMENSION(:,:), POINTER :: PDT ! thickness of each layer
1658 !
1659 !* local variables
1660 ! ---------------
1661 INTEGER(KIND=kindOfInt) :: IRET ! return code
1662 INTEGER :: ILEV1 ! level definition
1663 INTEGER :: ILEV2 ! level definition
1664 REAL, DIMENSION(:), POINTER :: ZFIELD => null()
1665 INTEGER :: JL ! layer loop counter
1666 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1667 !--------------------------------------------------------------------------------
1668 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TG_HIRLAM',0,zhook_handle)
1669 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_TG_HIRLAM: | Reading soil temperature'
1670 !--------------------------------------------------------------------------------
1671 ! 1. Search and read level 1 (and its depth)
1672 ! -----------------------
1673 ilev1 = 904
1674 ilev2 = -1
1675  CALL read_grib(hgrib,kluout,11,iret,zfield,klev1=ilev1,klev2=ilev2)
1676 IF (iret /= 0 ) THEN
1677  CALL abor1_sfx('MODE_READ_GRIB: SOIL TEMPERATURE LEVEL 1 MISSING (READ_GRIB_TG_HIRLAM)')
1678 END IF
1679 !
1680 ALLOCATE(ptg(SIZE(zfield),3))
1681 ALLOCATE(pdt(SIZE(zfield),3))
1682 ptg(:,1)= zfield(:)
1683 pdt(:,1) = 0.01
1684 !--------------------------------------------------------------------------------
1685 ! 2. Deep soil temperature
1686 ! ---------------------
1687 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_TG_HIRLAM: | Reading deep soil temperature'
1688 !
1689 ilev1 = 954
1690 ilev2 = -1
1691  CALL read_grib(hgrib,kluout,11,iret,zfield,klev1=ilev1,klev2=ilev2)
1692 IF (iret /= 0) THEN
1693  CALL abor1_sfx('MODE_READ_GRIB: DEEP SOIL TEMPERATURE MISSING (READ_GRIB_TG_HIRLAM)')
1694 END IF
1695 !
1696 ptg(:,2)= zfield(:)
1697 DEALLOCATE(zfield)
1698 pdt(:,2) = 0.4 ! deep temperature layer depth assumed equal to 0.40m
1699 !--------------------------------------------------------------------------------
1700 ! 4. Assumes uniform temperature profile below
1701 ! -----------------------------------------
1702 ptg(:,3) = ptg(:,2)
1703 pdt(:,3) = 5. ! temperature profile down to 5m
1704 !--------------------------------------------------------------------------------
1705 ! 5. Apply land mask
1706 ! ---------------
1707 IF (SIZE(pmask)==SIZE(ptg,1)) THEN
1708  DO jl=1,SIZE(ptg,2)
1709  WHERE (pmask(:)/=1.) ptg(:,jl) = xundef
1710  END DO
1711 ENDIF
1712 !
1713 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TG_HIRLAM',1,zhook_handle)
1714 END SUBROUTINE read_grib_tg_hirlam
1715 !-------------------------------------------------------------------
1716 ! #######################
1717  SUBROUTINE read_grib_wg_hirlam(HGRIB,KLUOUT,HINMODEL,PMASK,PFIELD,PD)
1718 ! #######################
1719 !
1720 USE modd_grid_grib, ONLY : nni
1721 USE modd_surf_par, ONLY : xundef
1722 !
1723 IMPLICIT NONE
1724 !
1725 !* dummy arguments
1726 ! ---------------
1727  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1728 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1729  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1730 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1731 REAL, DIMENSION(:,:), POINTER :: PFIELD ! field to initialize
1732 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
1733 !
1734 !* local variables
1735 ! ---------------
1736 INTEGER(KIND=kindOfInt) :: IRET ! return code
1737 INTEGER :: ILTYPE ! type of level (Grib code table 3)
1738 INTEGER :: ILEV1 ! level definition
1739 INTEGER :: ILEV2 ! level definition
1740 REAL, DIMENSION(:), POINTER :: ZFIELD => null()
1741 REAL, DIMENSION(:,:), ALLOCATABLE :: ZWG ! first water reservoir
1742 REAL, DIMENSION(:), ALLOCATABLE :: ZD ! Height of each layer
1743 INTEGER :: INLAYERDEEP! number of deep moisture layers
1744 REAL :: ZWWILT ! ECMWF wilting point
1745 REAL :: ZWFC ! ECMWF field capacity
1746 REAL :: ZWSAT ! ECMWF saturation
1747 INTEGER :: JL ! loop counter on layers
1748 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1749 !--------------------------------------------------------------------------------
1750 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WG_HIRLAM',0,zhook_handle)
1751 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_WG_HIRLAM: | Reading soil moisture'
1752 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_WG_HIRLAM: | WARNING READING LOW VEGETATION TILE (NR 4) ONLY'
1753 !
1754 ALLOCATE(zd(2))
1755 !
1756 ! 1. Search and read level 1 (and its depth)
1757 ! -----------------------
1758 iltype=105
1759 ilev1=904
1760 ilev2=-1
1761  CALL read_grib(hgrib,kluout,86,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
1762 IF (iret /= 0 ) THEN
1763  CALL abor1_sfx('MODE_READ_GRIB: SOIL MOISTURE LEVEL 1 MISSING (READ_GRIB_WG_HIRLAM)')
1764 END IF
1765 !
1766 ALLOCATE(zwg(SIZE(zfield),2))
1767 zwg(:,1)=zfield
1768 !
1769 zd(1) = 0.01
1770 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_WG_HIRLAM: | Level 1 height set to 0.01 m '
1771 !
1772 zwg(:,1) = zwg(:,1) / zd(1) ! convert units to m3/m3
1773 !--------------------------------------------------------------------------------
1774 ! 2. Search and read level 2 (and its depth) This level is optionnal
1775 ! -----------------------
1776 iltype=105
1777 ilev1=954
1778 ilev2=-1
1779  CALL read_grib(hgrib,kluout,86,iret,zfield,kltype=iltype,klev1=ilev1,klev2=ilev2)
1780 IF (iret /= 0 ) THEN
1781  CALL abor1_sfx('MODE_READ_GRIB: SOIL MOISTURE LEVEL 2 MISSING (READ_GRIB_WG_HIRLAM)')
1782 END IF
1783 !
1784 zwg(:,2)=zfield ! already units m3/m3
1785 DEALLOCATE(zfield)
1786 !
1787 inlayerdeep = 2
1788 zd(2) = 0.42
1789 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_WG_HIRLAM: | Level 2 height set to 0.42 m '
1790 !
1791 WRITE (kluout,'(A)') 'WARNING MODE_READ_GRIB: ZWG3 AND ZWG4 SET TO 0. (READ_GRIB_WG_HIRLAM)'
1792 !--------------------------------------------------------------------------------
1793 ! 3. Set water content profile and layer thicknesses
1794 ! -----------------------------------------------
1795  CALL fill_pfield(kluout,'READ_GRIB_WG_HIRLAM',inlayerdeep,zd,zwg,pmask,pfield,pd)
1796 DEALLOCATE(zd)
1797 DEALLOCATE(zwg)
1798 !--------------------------------------------------------------------------------
1799 ! 4. Convert from specific humidity to relative humidity
1800 ! ---------------------------------------------------
1801 ! Compute model's constants
1802 zwfc = 0.171
1803 zwwilt = 0.086
1804 !
1805 ! Then perform conversion
1806 DO jl=1,inlayerdeep
1807  WHERE (pfield(:,jl).NE.xundef) pfield(:,jl) = (pfield(:,jl) - zwwilt) / (zwfc - zwwilt)
1808 ENDDO
1809 !
1810 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WG_HIRLAM',1,zhook_handle)
1811 END SUBROUTINE read_grib_wg_hirlam
1812 !-------------------------------------------------------------------
1813 ! #######################
1814  SUBROUTINE read_grib_wgi_hirlam(HGRIB,KLUOUT,PFIELD,PD)
1815 ! #######################
1816 !
1817 USE modd_grid_grib, ONLY : nni
1818 !
1819 IMPLICIT NONE
1820 !
1821 !* dummy arguments
1822 ! ---------------
1823  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1824 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1825 REAL, DIMENSION(:,:), POINTER :: PFIELD ! field to initialize
1826 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
1827 !
1828 !* local variables
1829 ! ---------------
1830 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1831 !--------------------------------------------------------------------------------
1832 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WGI_HIRLAM',0,zhook_handle)
1833 !
1834 ALLOCATE (pfield(nni,2))
1835 ALLOCATE (pd(nni,2))
1836 pfield(:,:) = 0.
1837 !
1838 pd(:,1) = 0.01
1839 pd(:,2) = 1.
1840 !
1841 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_WGI_HIRLAM',1,zhook_handle)
1842 END SUBROUTINE read_grib_wgi_hirlam
1843 !-------------------------------------------------------------------
1844 !-------------------------------------------------------------------
1845 ! #######################
1846  SUBROUTINE read_grib_snow_veg_and_depth(HGRIB,KLUOUT,HINMODEL,PMASK,PSNV,PSNVD)
1847 ! #######################
1848 !!
1849 !! MODIFICATIONS
1850 !! -------------
1851 !! C Ardilouze 07/2013 : possibility to read snow density (ERAI-land)
1852 !!
1853 !-------------------------------------------------------------------------------
1854 !
1855 USE modd_grid_grib, ONLY : nni
1856 USE modd_surf_par, ONLY : xundef
1857 USE modd_snow_par, ONLY : xrhosmax
1858 !
1859 IMPLICIT NONE
1860 !
1861 !* dummy arguments
1862 ! ---------------
1863  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1864 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1865  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1866 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1867 REAL, DIMENSION(:), OPTIONAL, POINTER :: PSNV ! field to initialize
1868 REAL, DIMENSION(:), OPTIONAL, POINTER :: PSNVD ! field to initialize
1869 !
1870 !* local variables
1871 ! ---------------
1872 INTEGER(KIND=kindOfInt) :: IRET ! return code
1873 REAL, DIMENSION(:), POINTER :: ZFIELD => null() ! field to initialize
1874 REAL, DIMENSION(:), POINTER :: ZFIELD2 => null() ! field to initialize
1875 REAL, DIMENSION(:), POINTER :: ZRHO => null()
1876 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1877 !--------------------------------------------------------------------------------
1878 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SNOW_VEG_AND_DEPTH',0,zhook_handle)
1879 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_SNOW_VEG_AND_DEPTH: | Reading snow depth and density (if present)'
1880 !
1881 SELECT CASE(hinmodel)
1882  CASE('ECMWF ')
1883  CALL read_grib(hgrib,kluout,141,iret,zfield)
1884  CASE('ARPEGE','ALADIN','MOCAGE','HIRLAM')
1885  CALL read_grib(hgrib,kluout,66,iret,zfield)
1886  CASE DEFAULT
1887  CALL abor1_sfx('MODE_READ_GRIB:READ_GRIB_SNOW_VEG_AND_DEPTH: OPTION NOT SUPPORTED '//hinmodel)
1888 END SELECT
1889 !
1890 IF (iret /= 0 ) THEN
1891  CALL abor1_sfx('MODE_READ_GRIB: SNOW AND VEG DEPTH MISSING (READ_GRIB_SNOW_VEG_AND_DEPTH)')
1892 END IF
1893 !
1894  CALL read_grib_snow_den(hgrib,kluout,hinmodel,pmask,zrho)
1895 !
1896 IF (PRESENT(psnv)) THEN
1897  ALLOCATE(psnv(SIZE(zfield)))
1898  psnv(:)=zfield(:)
1899  IF (hinmodel=='ECMWF ') psnv(:) = psnv(:) * zrho(:)
1900  IF (SIZE(pmask)==SIZE(psnv)) &
1901  WHERE (pmask(:)/=1.) psnv(:) = xundef
1902 ENDIF
1903 !
1904 IF (PRESENT(psnvd)) THEN
1905  ALLOCATE(psnvd(SIZE(zfield)))
1906  psnvd(:)=zfield(:)
1907  IF (hinmodel/='ECMWF ') psnvd = psnvd / zrho(:)
1908  IF (SIZE(pmask)==SIZE(psnvd)) &
1909  WHERE (pmask(:)/=1.) psnvd(:) = xundef
1910 ENDIF
1911 !
1912 DEALLOCATE(zfield)
1913 DEALLOCATE(zrho)
1914 !
1915 !
1916 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SNOW_VEG_AND_DEPTH',1,zhook_handle)
1917 END SUBROUTINE read_grib_snow_veg_and_depth
1918 !-------------------------------------------------------------------
1919 !-------------------------------------------------------------------
1920 ! #######################
1921  SUBROUTINE read_grib_snow_alb(HGRIB,KLUOUT,HINMODEL,PMASK,PSNVA)
1922 ! #######################
1923 !!
1924 !! AUTHOR
1925 !! -------------
1926 !! C Ardilouze 07/2013 : possibility to read snow albedo (ERAI-land)
1927 !!
1928 !-------------------------------------------------------------------------------
1929 !
1930 USE modd_grid_grib, ONLY : nni
1931 USE modd_surf_par, ONLY : xundef
1932 USE modd_snow_par, ONLY : xansmin, xansmax
1933 !
1934 IMPLICIT NONE
1935 !
1936 !* dummy arguments
1937 ! ---------------
1938  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1939 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1940  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1941 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1942 REAL, DIMENSION(:), POINTER :: PSNVA ! field to initialize
1943 !
1944 !* local variables
1945 ! ---------------
1946 INTEGER(KIND=kindOfInt) :: IRET ! return code
1947 REAL, DIMENSION(:), POINTER :: ZFIELD => null() ! field to initialize
1948 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1949 !--------------------------------------------------------------------------------
1950 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SNOW_ALB',0,zhook_handle)
1951 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_SNOW_ALB: | Reading snow albedo'
1952 !
1953 ALLOCATE(psnva(nni))
1954 psnva(:) = 0.5 * ( xansmin + xansmax )
1955 IF (hinmodel == 'ECMWF') THEN
1956  CALL read_grib(hgrib,kluout,32,iret,zfield)
1957  IF (iret == 0 ) THEN
1958  DEALLOCATE(psnva)
1959  ALLOCATE(psnva(SIZE(zfield)))
1960  psnva(:)=zfield(:)
1961  DEALLOCATE(zfield)
1962  END IF
1963 END IF
1964 !
1965 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SNOW_ALB',1,zhook_handle)
1966 END SUBROUTINE read_grib_snow_alb
1967 !!-------------------------------------------------------------------
1968 ! #######################
1969  SUBROUTINE read_grib_snow_den(HGRIB,KLUOUT,HINMODEL,PMASK,PSNV)
1970 ! #######################
1971 !!
1972 !! AUTHOR
1973 !! -------------
1974 !! C Ardilouze 08/2013 : possibility to read snow density (ERAI-land)
1975 !!
1976 !-------------------------------------------------------------------------------
1977 !
1978 USE modd_grid_grib, ONLY : nni
1979 USE modd_surf_par, ONLY : xundef
1980 USE modd_snow_par, ONLY : xrhosmax
1981 !
1982 IMPLICIT NONE
1983 !
1984 !* dummy arguments
1985 ! ---------------
1986  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
1987 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
1988  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
1989 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
1990 REAL, DIMENSION(:), POINTER :: PSNV ! field to initialize
1991 !
1992 !* local variables
1993 ! ---------------
1994 INTEGER(KIND=kindOfInt) :: IRET ! return code
1995 REAL, DIMENSION(:), POINTER :: ZFIELD => null() ! field to initialize
1996 REAL(KIND=JPRB) :: ZHOOK_HANDLE
1997 !--------------------------------------------------------------------------------
1998 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SNOW_DEN',0,zhook_handle)
1999 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_SNOW_DEN: | Reading snow density'
2000 !
2001 ALLOCATE(psnv(nni))
2002 psnv(:) = xrhosmax
2003 IF (hinmodel == 'ECMWF') THEN
2004  CALL read_grib(hgrib,kluout,33,iret,zfield)
2005  IF (iret == 0 ) THEN
2006  DEALLOCATE(psnv)
2007  ALLOCATE(psnv(SIZE(zfield)))
2008  psnv(:)=zfield(:)
2009  DEALLOCATE(zfield)
2010  END IF
2011 END IF
2012 !
2013 IF (SIZE(pmask)==SIZE(psnv)) &
2014  WHERE (pmask(:)/=1.) psnv = xundef
2015 !
2016 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_SNOW_DEN',1,zhook_handle)
2017 END SUBROUTINE read_grib_snow_den
2018 !-------------------------------------------------------------------
2019 !-------------------------------------------------------------------
2020 ! #######################
2021  SUBROUTINE read_grib_t_teb(HGRIB,KLUOUT,HINMODEL,PTI,PMASK,PT,PD)
2022 ! #######################
2023 !
2024 USE modd_grid_grib, ONLY : nni
2025 USE modd_surf_par, ONLY : xundef
2026 !
2027 IMPLICIT NONE
2028 !
2029 !* dummy arguments
2030 ! ---------------
2031  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
2032 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
2033  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
2034 REAL, INTENT(IN) :: PTI ! internal temperature
2035 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
2036 REAL, DIMENSION(:,:), POINTER :: PT ! field to initialize
2037 REAL, DIMENSION(:,:), POINTER :: PD ! normalized grid
2038 !
2039 !* local variables
2040 ! ---------------
2041 REAL, DIMENSION(:), POINTER :: ZFIELD => null() ! field to initialize
2042 INTEGER :: JL ! layer loop counter
2043 REAL(KIND=JPRB) :: ZHOOK_HANDLE
2044 !--------------------------------------------------------------------------------
2045 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_T_TEB',0,zhook_handle)
2046 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_T_TEB: | Reading temperature for buildings'
2047 !
2048  CALL read_grib_ts(hgrib,kluout,hinmodel,pmask,zfield)
2049 !
2050 ALLOCATE(pt(SIZE(zfield),3))
2051 ALLOCATE(pd(SIZE(zfield),3))
2052 !
2053 pt(:,1) = zfield
2054 DEALLOCATE(zfield)
2055 pd(:,1) = 0.
2056 !
2057 pt(:,2) = pti
2058 pd(:,2) = 0.5 ! deep temperature depth assumed at half of wall or roof
2059 !
2060 pt(:,3) = pti
2061 pd(:,3) = 1. ! temperature at building interior
2062 !
2063 IF (SIZE(pmask)==SIZE(pt,1)) THEN
2064  DO jl=1,SIZE(pt,2)
2065  WHERE (pmask(:)/=1.) pt(:,jl) = xundef
2066  END DO
2067 ENDIF
2068 !
2069 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_T_TEB',1,zhook_handle)
2070 END SUBROUTINE read_grib_t_teb
2071 !-------------------------------------------------------------------
2072 ! #######################
2073  SUBROUTINE read_grib_tf_teb(HGRIB,KLUOUT,HINMODEL,PTI,PMASK,PTF,PD)
2074 ! #######################
2075 !
2076 USE modd_surf_par, ONLY : xundef
2077 !
2078 IMPLICIT NONE
2079 !
2080 !* dummy arguments
2081 ! ---------------
2082  CHARACTER(LEN=*), INTENT(IN) :: HGRIB ! Grib file name
2083 INTEGER, INTENT(IN) :: KLUOUT ! logical unit of output listing
2084  CHARACTER(LEN=6), INTENT(IN) :: HINMODEL ! Grib originating model
2085 REAL, INTENT(IN) :: PTI ! internal temperature
2086 REAL, DIMENSION(:), INTENT(IN) :: PMASK ! grib land mask
2087 REAL, DIMENSION(:,:), POINTER :: PTF ! field to initialize
2088 REAL, DIMENSION(:,:), POINTER :: PD ! thickness of each layer
2089 !
2090 !* local variables
2091 ! ---------------
2092 INTEGER(KIND=kindOfInt) :: IRET ! return code
2093 REAL, DIMENSION(:), POINTER :: ZFIELD => null() ! field to read
2094 INTEGER :: JL ! layer loop counter
2095 REAL(KIND=JPRB) :: ZHOOK_HANDLE
2096 !--------------------------------------------------------------------------------
2097 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TF_TEB',0,zhook_handle)
2098 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_TF_TEB: | Reading temperature for building floor'
2099 !
2100 ! 1. Deep soil temperature
2101 ! ---------------------
2102 !
2103 WRITE (kluout,'(A)') 'MODE_READ_GRIB:READ_GRIB_TF_TEB: | Reading deep soil temperature'
2104 !
2105  CALL read_grib_t2_land(hgrib,kluout,hinmodel,pmask,zfield)
2106 !
2107 ALLOCATE(ptf(SIZE(zfield),3))
2108 ALLOCATE(pd(SIZE(zfield),3))
2109 !
2110 ptf(:,2) = zfield(:)
2111 pd(:,2) = 0.5 ! deep temperature depth assumed at half of the floor
2112 !
2113 DEALLOCATE(zfield)
2114 !
2115 ! 2. level 1 is the internal building temperature
2116 ! -----------------------
2117 !
2118 ptf(:,1) = pti
2119 pd(:,1) = 0.
2120 !
2121 ! 3. Assumes uniform temperature profile below
2122 ! -----------------------------------------
2123 !
2124 ptf(:,3) = ptf(:,2)
2125 pd(:,3) = 1. ! deep temperature value
2126 !
2127 IF (lhook) CALL dr_hook('MODE_READ_GRIB:READ_GRIB_TF_TEB',1,zhook_handle)
2128 !------------------------------------------------------------------------------
2129 END SUBROUTINE read_grib_tf_teb
2130 !-------------------------------------------------------------------
2131 END MODULE mode_read_grib
mode_read_grib::fill_pfield
subroutine fill_pfield(KLUOUT, HROUT, KNLAYERDEEP, PDIN, PFIELDIN, PMASK, PFIELDOUT, PDOUT)
Definition: mode_read_grib.F90:648
mode_read_grib::get_grib_message
subroutine get_grib_message(KLUOUT, KLTYPE, KLEV1, KLEV2, KGRIB, KFOUND)
Definition: mode_read_grib.F90:71
trim
static const char * trim(const char *name, int *n)
Definition: drhook.c:2383
modd_grid_grib::cgrib_file
character(len=28) cgrib_file
Definition: modd_grid_grib.F90:41
mode_read_grib::read_grib_tswater
subroutine read_grib_tswater(HGRIB, KLUOUT, HINMODEL, PMASK, PTS)
Definition: mode_read_grib.F90:505
modd_surf_par
Definition: modd_surf_par.F90:6
mode_read_grib::read_grib_t_teb
subroutine read_grib_t_teb(HGRIB, KLUOUT, HINMODEL, PTI, PMASK, PT, PD)
Definition: mode_read_grib.F90:2022
mode_read_grib::read_grib_zs
subroutine read_grib_zs(HGRIB, KLUOUT, HINMODEL, PZS)
Definition: mode_read_grib.F90:303
yomhook::dr_hook
Definition: yomhook.F90:20
mode_read_grib::read_grib_tf_teb
subroutine read_grib_tf_teb(HGRIB, KLUOUT, HINMODEL, PTI, PMASK, PTF, PD)
Definition: mode_read_grib.F90:2074
mode_read_grib::read_grib_tg_meteo_france
subroutine read_grib_tg_meteo_france(HGRIB, KLUOUT, HINMODEL, PMASK, PTG, PDT)
Definition: mode_read_grib.F90:1241
modd_grid_grib::nni
integer nni
Definition: modd_grid_grib.F90:38
mode_read_grib::read_grib_zs_land
subroutine read_grib_zs_land(HGRIB, KLUOUT, HINMODEL, PMASK, PZSL)
Definition: mode_read_grib.F90:344
mode_read_grib::ecmwf_wgi
subroutine ecmwf_wgi(PTG, PWG, PWGI)
Definition: mode_read_grib.F90:944
modd_csts::xpi
real, save xpi
Definition: modd_csts.F90:43
mode_read_grib::read_grib_t2
subroutine read_grib_t2(HGRIB, KLUOUT, HINMODEL, PMASK, PT2)
Definition: mode_read_grib.F90:539
mode_read_grib::read_grib_wgi_ecmwf
subroutine read_grib_wgi_ecmwf(HGRIB, KLUOUT, HINMODEL, PMASK, PFIELD, PD)
Definition: mode_read_grib.F90:1158
mode_read_grib::read_grib_ts
subroutine read_grib_ts(HGRIB, KLUOUT, HINMODEL, PMASK, PTS)
Definition: mode_read_grib.F90:446
abor1_sfx
subroutine abor1_sfx(YTEXT)
Definition: abor1_sfx.F90:7
modd_surf_par::xundef
real, parameter xundef
Definition: modd_surf_par.F90:43
mode_read_grib::read_grib_tg_hirlam
subroutine read_grib_tg_hirlam(HGRIB, KLUOUT, HINMODEL, PMASK, PTG, PDT)
Definition: mode_read_grib.F90:1643
mode_read_grib::read_grib_t2_land
subroutine read_grib_t2_land(HGRIB, KLUOUT, HINMODEL, PMASK, ZFIELD)
Definition: mode_read_grib.F90:587
modd_csts::xg
real, save xg
Definition: modd_csts.F90:55
mode_read_grib::read_grib_snow_alb
subroutine read_grib_snow_alb(HGRIB, KLUOUT, HINMODEL, PMASK, PSNVA)
Definition: mode_read_grib.F90:1922
mode_read_grib::read_grib_t
subroutine read_grib_t(HGRIB, KLUOUT, HINMODEL, PT)
Definition: mode_read_grib.F90:392
parkind1::jprb
integer, parameter jprb
Definition: parkind1.F90:32
mode_read_grib::read_grib_sand_clay_meteo_france
subroutine read_grib_sand_clay_meteo_france(HGRIB, KLUOUT, HINMODEL, PSAND, PCLAY, GISBA)
Definition: mode_read_grib.F90:1295
mode_read_grib::read_grib
subroutine read_grib(HGRIB, KLUOUT, KPARAM, KRET, PFIELD, KLTYPE, KLEV1, KLEV2, KPARAM2)
Definition: mode_read_grib.F90:175
mode_read_grib::read_grib_zs_sea
subroutine read_grib_zs_sea(HGRIB, KLUOUT, HINMODEL, PMASK, PZSS)
Definition: mode_read_grib.F90:368
modd_grid_grib::nidx
integer(kind=kindofint) nidx
Definition: modd_grid_grib.F90:42
modd_grid_grib
Definition: modd_grid_grib.F90:6
mode_read_grib::read_grib_snow_veg_and_depth
subroutine read_grib_snow_veg_and_depth(HGRIB, KLUOUT, HINMODEL, PMASK, PSNV, PSNVD)
Definition: mode_read_grib.F90:1847
mode_read_grib::harmonize_grib_wg_wgi_ecmwf
subroutine harmonize_grib_wg_wgi_ecmwf(HGRIB, KLUOUT, HINMODEL, PMASK, PWG, PD, PWGI)
Definition: mode_read_grib.F90:988
mode_read_grib::read_grib_snow_den
subroutine read_grib_snow_den(HGRIB, KLUOUT, HINMODEL, PMASK, PSNV)
Definition: mode_read_grib.F90:1970
mode_read_grib::read_grib_sst
subroutine read_grib_sst(HGRIB, KLUOUT, HINMODEL, PMASK, PSST)
Definition: mode_read_grib.F90:471
sum
intent(out) overrides sub arrays one Sort by the least significant key first sum(iindex(1:n))
yomhook::lhook
logical lhook
Definition: yomhook.F90:15
test_iret
subroutine test_iret(KLUOUT, VAL1, VAL0, KRET)
Definition: mode_read_grib.F90:143
mode_read_grib::read_grib_land_mask
subroutine read_grib_land_mask(HGRIB, KLUOUT, HINMODEL, PMASK)
Definition: mode_read_grib.F90:258
mode_read_grib::read_grib_wg_ecmwf_1
subroutine read_grib_wg_ecmwf_1(HGRIB, KLUOUT, HINMODEL, PMASK, PWG, PD)
Definition: mode_read_grib.F90:806
mode_read_grib::read_grib_wg_meteo_france
subroutine read_grib_wg_meteo_france(HGRIB, KLUOUT, HINMODEL, PMASK, PFIELD, PD)
Definition: mode_read_grib.F90:1350
mode_read_grib::read_grib_tg_ecmwf
subroutine read_grib_tg_ecmwf(HGRIB, KLUOUT, HINMODEL, PMASK, PTG, PD)
Definition: mode_read_grib.F90:692
mode_read_grib::read_grib_wg_hirlam
subroutine read_grib_wg_hirlam(HGRIB, KLUOUT, HINMODEL, PMASK, PFIELD, PD)
Definition: mode_read_grib.F90:1718
mode_read_grib::read_grib_wgi_hirlam
subroutine read_grib_wgi_hirlam(HGRIB, KLUOUT, PFIELD, PD)
Definition: mode_read_grib.F90:1815
mode_read_grib::put_layer_depth
subroutine put_layer_depth(KLUOUT, KLEV, HROUT, KLTYPE, KLEV1, KLEV2, KNLAYERDEEP, PV4, PV, PD)
Definition: mode_read_grib.F90:613
mode_read_grib::make_grib_index
subroutine make_grib_index(HGRIB)
Definition: mode_read_grib.F90:20
modd_csts
Definition: modd_csts.F90:6
modd_csts::xtt
real, save xtt
Definition: modd_csts.F90:66
mode_read_grib::clear_grib_index
subroutine clear_grib_index
Definition: mode_read_grib.F90:47
mode_read_grib::read_grib_wgi_meteo_france
subroutine read_grib_wgi_meteo_france(HGRIB, KLUOUT, HINMODEL, PMASK, PFIELD, PD)
Definition: mode_read_grib.F90:1502
mode_read_grib
Definition: mode_read_grib.F90:6
mode_read_grib::read_grib_wg_ecmwf
subroutine read_grib_wg_ecmwf(HGRIB, KLUOUT, HINMODEL, PMASK, PFIELD, PD)
Definition: mode_read_grib.F90:1065
Generated on Tue Jan 16 2018 16:23:29 for SURFEX v8.1 by   doxygen 1.8.13