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