Juan 24/11/2015: modif for PREPLL from M.Mogié
[MNH-git_open_source-lfs.git] / src / SURFEX / mode_read_extern.F90
1 !SURFEX_LIC Copyright 1994-2014 Meteo-France 
2 !SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C  licence
3 !SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !SURFEX_LIC for details. version 1.
5 ! Modifications :
6 ! P.Marguinaud : 11-09-2012 : shorten field name
7 ! G.Delautier : 24-06-2015 : bug for arome compressed files
8 !     #####################
9 MODULE MODE_READ_EXTERN
10 !     #####################
11 !-------------------------------------------------------------------
12 !
13 USE MODI_READ_LECOCLIMAP
14 !
15 USE MODI_PUT_ON_ALL_VEGTYPES
16 USE MODI_OLD_NAME
17 !
18 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
19 USE PARKIND1  ,ONLY : JPRB
20 !
21 CONTAINS
22 !
23 !---------------------------------------------------------------------------------------
24 !
25 !     #######################
26       SUBROUTINE READ_EXTERN_DEPTH(HPROGRAM,KLUOUT,HISBA,HNAT,HFIELD,KNI,KLAYER, &
27                                    KPATCH,PSOILGRID,PDEPTH,KVERSION  )
28 !     #######################
29 !
30 USE MODD_SURF_PAR,       ONLY : NUNDEF, XUNDEF
31 USE MODD_DATA_COVER_PAR, ONLY : JPCOVER, NVEGTYPE
32 !
33 USE MODI_READ_SURF_ISBA_PAR_n
34 USE MODI_READ_SURF
35 USE MODI_CONVERT_COVER_ISBA
36 USE MODI_GARDEN_SOIL_DEPTH
37
38 !
39 IMPLICIT NONE
40 !
41 !* dummy arguments
42 !  ---------------
43 !
44  CHARACTER(LEN=6),     INTENT(IN)    :: HPROGRAM  ! type of input file
45 INTEGER,              INTENT(IN)    :: KLUOUT    ! logical unit of output listing
46  CHARACTER(LEN=3),     INTENT(IN)    :: HISBA     ! type of ISBA soil scheme
47  CHARACTER(LEN=3),     INTENT(IN)    :: HNAT      ! type of surface (nature, gardens)
48  CHARACTER(LEN=7),     INTENT(IN)    :: HFIELD    ! field name
49 INTEGER,              INTENT(IN)    :: KNI       ! number of points
50 INTEGER,           INTENT(INOUT)    :: KLAYER    ! number of layers
51 INTEGER,              INTENT(IN)    :: KPATCH    ! number of patch
52 INTEGER,              INTENT(IN)    :: KVERSION  ! surface version
53 REAL, DIMENSION(:),   INTENT(IN)    :: PSOILGRID
54 REAL, DIMENSION(:,:,:), POINTER     :: PDEPTH    ! middle depth of each layer
55 !
56 !* local variables
57 !  ---------------
58 !
59  CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
60  CHARACTER(LEN=16) :: YRECFM2
61  CHARACTER(LEN=100):: YCOMMENT       ! Comment string
62 INTEGER           :: IRESP          ! reading return code
63 INTEGER           :: ILAYER         ! number of soil layers
64 INTEGER           :: JLAYER         ! loop counter
65 INTEGER           :: JPATCH         ! loop counter
66 INTEGER           :: JJ
67 INTEGER           :: IVERSION
68 INTEGER           :: IBUGFIX
69 !
70 LOGICAL, DIMENSION(JPCOVER)          :: GCOVER ! flag to read the covers
71 REAL,    DIMENSION(:,:), ALLOCATABLE :: ZCOVER ! cover fractions
72 REAL,    DIMENSION(:,:), ALLOCATABLE :: ZGROUND_DEPTH ! cover fractions
73 INTEGER, DIMENSION(:,:), ALLOCATABLE :: IWG_LAYER
74 REAL,  DIMENSION(:,:,:), ALLOCATABLE :: ZD     ! depth of each inter-layer
75 REAL,  DIMENSION(:,:,:), ALLOCATABLE :: ZDG    ! depth of each inter-layer
76 REAL,  DIMENSION(:,:,:), ALLOCATABLE :: ZDEPTH ! middle of each layer for each patch
77 REAL,  DIMENSION(:,:),   ALLOCATABLE :: ZWORK  ! work array
78 REAL,  DIMENSION(KNI)                :: ZHVEG  ! high vegetation fraction
79 REAL,  DIMENSION(KNI)                :: ZLVEG  ! low  vegetation fraction
80 REAL,  DIMENSION(KNI)                :: ZNVEG  ! no   vegetation fraction
81  CHARACTER(LEN=4)                     :: YHVEG  ! type of high vegetation
82  CHARACTER(LEN=4)                     :: YLVEG  ! type of low  vegetation
83  CHARACTER(LEN=4)                     :: YNVEG  ! type of no   vegetation
84 LOGICAL                              :: GECOCLIMAP ! T if ecoclimap is used
85 LOGICAL                              :: GPAR_GARDEN! T if garden data are used
86 LOGICAL                              :: GDATA_DG
87 LOGICAL                              :: GDATA_GROUND_DEPTH
88 INTEGER                              :: IHYDRO_LAYER
89 REAL(KIND=JPRB) :: ZHOOK_HANDLE
90 !
91 !
92 !------------------------------------------------------------------------------
93 !
94 IF (LHOOK) CALL DR_HOOK('MODE_READ_EXTERN:READ_EXTERN_DEPTH',0,ZHOOK_HANDLE)
95 !
96 IF (HNAT=='NAT') THEN
97   CALL READ_LECOCLIMAP(HPROGRAM,GECOCLIMAP)
98 ELSE
99   CALL READ_SURF(HPROGRAM,'PAR_GARDEN',GPAR_GARDEN,IRESP)
100   GECOCLIMAP = .NOT. GPAR_GARDEN
101 END IF
102 !
103 !
104 YRECFM='VERSION'
105  CALL READ_SURF(HPROGRAM,YRECFM,IVERSION,IRESP)
106 !
107 YRECFM='BUG'
108  CALL READ_SURF(HPROGRAM,YRECFM,IBUGFIX,IRESP)
109 !
110 !------------------------------------------------------------------------------
111 !
112 ALLOCATE(ZDG   (KNI,KLAYER,KPATCH))
113 ALLOCATE(IWG_LAYER   (KNI,KPATCH))
114 IWG_LAYER(:,:) = NUNDEF
115 IHYDRO_LAYER = KLAYER
116 !
117 IF (GECOCLIMAP) THEN
118
119  IF (IVERSION<7 .OR. IVERSION==7 .AND. IBUGFIX<=3) THEN
120   !
121   !* reading of the cover to obtain the depth of inter-layers
122   !
123   CALL OLD_NAME(HPROGRAM,'COVER_LIST      ',YRECFM)
124   CALL READ_SURF(HPROGRAM,YRECFM,GCOVER(:),IRESP,HDIR='-')
125   !
126   ALLOCATE(ZCOVER(KNI,JPCOVER))
127   YRECFM='COVER'
128   CALL READ_SURF(HPROGRAM,YRECFM,ZCOVER(:,:),GCOVER(:),IRESP,HDIR='A')  
129   !
130   !* computes soil layers
131   !  
132   CALL CONVERT_COVER_ISBA(HISBA,NUNDEF,ZCOVER,'   ',HNAT,PSOILGRID=PSOILGRID,PDG=ZDG,KWG_LAYER=IWG_LAYER)
133   !
134   DEALLOCATE(ZCOVER)
135  ELSE
136 print*, '-----------------------------------------------'
137 print*, '-----------------------------------------------'
138 print*, '-----------------------------------------------'
139 print*, '-----------------------------------------------'
140 print*, 'MODE_READ_EXTERN : ==> ON NE LIT PAS LES COVERS'
141 print*, '-----------------------------------------------'
142 print*, '-----------------------------------------------'
143 print*, '-----------------------------------------------'
144 print*, '-----------------------------------------------'
145 #ifdef MNH_PARALLEL
146 DO JPATCH=1,SIZE(ZDG,3)
147   DO JLAYER=1,SIZE(ZDG,2)
148     IF (JLAYER<10) THEN
149       IF (HNAT=='NAT') THEN
150         WRITE(YRECFM,FMT='(A6,I1,I4.4)') 'ECO_DG',JLAYER,JPATCH
151       ELSE
152         WRITE(YRECFM,FMT='(A9,I1,I4.4)') 'GD_ECO_DG',JLAYER,JPATCH
153       END IF
154     ELSE
155       IF (HNAT=='NAT') THEN
156         WRITE(YRECFM,FMT='(A6,I2,I4.4)') 'ECO_DG',JLAYER,JPATCH
157       ELSE
158         WRITE(YRECFM,FMT='(A9,I2,I4.4)') 'GD_ECO_DG',JLAYER,JPATCH
159       END IF
160     ENDIF
161     CALL READ_SURF(HPROGRAM,YRECFM,ZDG(:,JLAYER,JPATCH),IRESP,HDIR='A')
162   END DO
163 END DO
164 #else
165   DO JLAYER=1,SIZE(ZDG,2)
166     IF (JLAYER<10) THEN
167       IF (HNAT=='NAT') THEN
168         WRITE(YRECFM,FMT='(A6,I1)') 'ECO_DG',JLAYER
169       ELSE
170         WRITE(YRECFM,FMT='(A9,I1)') 'GD_ECO_DG',JLAYER
171       END IF
172     ELSE
173       IF (HNAT=='NAT') THEN
174         WRITE(YRECFM,FMT='(A6,I2)') 'ECO_DG',JLAYER
175       ELSE
176         WRITE(YRECFM,FMT='(A9,I2)') 'GD_ECO_DG',JLAYER
177       END IF
178     ENDIF
179     CALL READ_SURF(HPROGRAM,YRECFM,ZDG(:,JLAYER,:),IRESP,HDIR='A')
180   END DO
181 #endif
182   IF (HISBA=='DIF') THEN
183     YRECFM='ECO_WG_L'
184     IF (HNAT=='GRD') YRECFM='GD_ECO_WG_L'
185     ALLOCATE(ZWORK(KNI,KPATCH))
186     CALL READ_SURF(HPROGRAM,YRECFM,ZWORK(:,:),IRESP,HDIR='A')  
187     WHERE (ZWORK==XUNDEF) ZWORK=NUNDEF
188     IWG_LAYER=NINT(ZWORK)
189     DEALLOCATE(ZWORK)
190   END IF
191  END IF
192   !
193   IF (HISBA=='DIF') IHYDRO_LAYER = MAXVAL(IWG_LAYER(:,:),IWG_LAYER(:,:)/=NUNDEF)
194 ENDIF
195
196 !-------------------------------------------------------------------
197 IF (HNAT=='NAT' .AND. (IVERSION>=7 .OR. .NOT.GECOCLIMAP)) THEN
198   !
199   !* directly read soil layers in the file for nature ISBA soil layers
200   !
201   GDATA_DG = .TRUE.
202   IF (IVERSION>=7) THEN
203     YRECFM='L_DG'
204     YCOMMENT=YRECFM
205     CALL READ_SURF(HPROGRAM,YRECFM,GDATA_DG,IRESP,HCOMMENT=YCOMMENT)
206   ENDIF
207   !
208   IF (GDATA_DG) THEN
209     !
210     ALLOCATE(ZWORK(KNI,KPATCH))
211     DO JLAYER=1,KLAYER
212       IF (JLAYER<10)  WRITE(YRECFM,FMT='(A4,I1.1)') 'D_DG',JLAYER
213       IF (JLAYER>=10) WRITE(YRECFM,FMT='(A4,I2.2)') 'D_DG',JLAYER
214       CALL READ_SURF_ISBA_PAR_n(HPROGRAM,YRECFM,KLUOUT,KNI,ZWORK,IRESP,IVERSION,HDIR='A')
215       DO JPATCH=1,KPATCH
216         ZDG(:,JLAYER,JPATCH) = ZWORK(:,JPATCH)
217       END DO
218     END DO
219     DEALLOCATE(ZWORK)
220     !
221   ENDIF
222   !
223     GDATA_GROUND_DEPTH=.FALSE.
224   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=2) THEN
225     !
226     YRECFM2='L_GROUND_DEPTH'
227     IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM2='L_GROUND_DPT'
228     YCOMMENT=YRECFM2
229     CALL READ_SURF(HPROGRAM,YRECFM2,GDATA_GROUND_DEPTH,IRESP,HCOMMENT=YCOMMENT)
230     !
231     IF (GDATA_GROUND_DEPTH) THEN
232       !
233       YRECFM2='D_GROUND_DETPH'
234       IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM2='D_GROUND_DPT'
235       ALLOCATE(ZGROUND_DEPTH(KNI,KPATCH))
236       CALL READ_SURF_ISBA_PAR_n(HPROGRAM,YRECFM2,KLUOUT,KNI,ZGROUND_DEPTH(:,:),IRESP,IVERSION,HDIR='A')
237       !
238       DO JPATCH=1,KPATCH
239         DO JJ=1,KNI
240           DO JLAYER=1,KLAYER
241             IF ( ZDG(JJ,JLAYER,JPATCH) <= ZGROUND_DEPTH(JJ,JPATCH) .AND. ZGROUND_DEPTH(JJ,JPATCH) < XUNDEF ) &
242                 IWG_LAYER(JJ,JPATCH) = JLAYER
243           ENDDO
244         ENDDO
245       ENDDO
246       DEALLOCATE(ZGROUND_DEPTH)
247       !
248       IF (HISBA=='DIF') IHYDRO_LAYER = MAXVAL(IWG_LAYER(:,:),IWG_LAYER(:,:)/=NUNDEF)
249       !
250     ENDIF
251     !
252   ENDIF
253   !
254 ELSE IF (HNAT=='GRD' .AND. .NOT.GECOCLIMAP) THEN
255   !
256   !* computes soil layers from vegetation fractions read in the file
257   !
258   CALL READ_SURF(HPROGRAM,'D_TYPE_HVEG',YHVEG,IRESP)
259   CALL READ_SURF(HPROGRAM,'D_TYPE_LVEG',YLVEG,IRESP)
260   CALL READ_SURF(HPROGRAM,'D_TYPE_NVEG',YNVEG,IRESP)
261   CALL READ_SURF(HPROGRAM,'D_FRAC_HVEG',ZHVEG,IRESP,HDIR='A')
262   CALL READ_SURF(HPROGRAM,'D_FRAC_LVEG',ZLVEG,IRESP,HDIR='A')
263   CALL READ_SURF(HPROGRAM,'D_FRAC_NVEG',ZNVEG,IRESP,HDIR='A')
264   ! Ground layers
265   CALL GARDEN_SOIL_DEPTH(YNVEG,YLVEG,YHVEG,ZNVEG,ZLVEG,ZHVEG,ZDG)
266   !
267 END IF
268 !
269 DEALLOCATE(IWG_LAYER)
270 !
271 IF (HFIELD=='WG    ' .OR. HFIELD=='WGI   ' .OR. HFIELD=='TWN_WG  ' .OR. HFIELD=='TWN_WGI ' .OR. &
272       HFIELD=='GD_WG  ' .OR. HFIELD=='GD_WGI ') THEN
273   KLAYER = IHYDRO_LAYER
274 ENDIF
275 !
276 !-------------------------------------------------------------------
277 !
278 !* In force-restore ISBA, adds a layer at bottom of surface layer and a layer
279 !  between root and deep layers.
280 !
281 IF (HISBA=='2-L' .OR. HISBA=='3-L') THEN
282   ILAYER = KLAYER + 1
283   IF (HISBA=='3-L') ILAYER = ILAYER + 1
284   ALLOCATE(ZD    (KNI,ILAYER,KPATCH))
285   DO JPATCH=1,KPATCH
286     ! for interpolations, middle of surface layer must be at least at 1cm
287     ZD(:,1,JPATCH) = MIN(3.*ZDG(:,1,JPATCH),MAX(ZDG(:,1,JPATCH),0.02))
288     ! new layer below surface layer. This layer will be at root depth layer humidity
289     ZD(:,2,JPATCH) = MIN(4.*ZDG(:,1,JPATCH),0.5*(ZDG(:,1,JPATCH)+ZDG(:,2,JPATCH)))
290     ! root layer
291     ZD(:,3,JPATCH) = ZDG(:,2,JPATCH)
292     IF (HISBA=='3-L') THEN
293       ! between root and deep layers. This layer will have deep soil humidity.
294       WHERE (ZDG(:,2,JPATCH)<ZDG(:,3,JPATCH))
295         ZD(:,4,JPATCH) = 0.75 * ZDG(:,2,JPATCH) + 0.25 * ZDG(:,3,JPATCH)
296       ELSEWHERE
297         ZD(:,4,JPATCH) = ZDG(:,3,JPATCH)
298       END WHERE
299       ! deep layer
300       ZD(:,5,JPATCH) = ZDG(:,3,JPATCH)
301     END IF
302   END DO
303 ELSE
304   ILAYER = KLAYER
305   ALLOCATE(ZD    (KNI,ILAYER,KPATCH))
306   ZD(:,:,:) = ZDG(:,1:KLAYER,:)
307 END IF
308 !
309 DEALLOCATE(ZDG)
310 !
311 !-------------------------------------------------------------------
312 !* recovers middle layer depth (from the surface)
313 ALLOCATE(ZDEPTH    (KNI,ILAYER,KPATCH))
314 ZDEPTH = XUNDEF
315 DO JPATCH=1,KPATCH
316   WHERE(ZD(:,1,JPATCH)/=XUNDEF) &
317     ZDEPTH    (:,1,JPATCH)=ZD(:,1,JPATCH)/2.  
318   DO JLAYER=2,ILAYER
319     WHERE(ZD(:,1,JPATCH)/=XUNDEF) &
320       ZDEPTH    (:,JLAYER,JPATCH) = (ZD(:,JLAYER-1,JPATCH) + ZD(:,JLAYER,JPATCH))/2.  
321   END DO
322 END DO
323 DEALLOCATE(ZD)
324 !-------------------------------------------------------------------
325 !
326 ALLOCATE(PDEPTH    (KNI,ILAYER,NVEGTYPE))
327  CALL PUT_ON_ALL_VEGTYPES(KNI,ILAYER,KPATCH,NVEGTYPE,ZDEPTH,PDEPTH)
328 DEALLOCATE(ZDEPTH)
329
330 IF (LHOOK) CALL DR_HOOK('MODE_READ_EXTERN:READ_EXTERN_DEPTH',1,ZHOOK_HANDLE)
331 !-------------------------------------------------------------------
332 !
333 END SUBROUTINE READ_EXTERN_DEPTH
334 !
335 !
336 !-------------------------------------------------------------------
337 !---------------------------------------------------------------------------------------
338 !
339 !     #######################
340       SUBROUTINE READ_EXTERN_ISBA(HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,&
341                                   KLUOUT,KNI,HFIELD,HNAME,PFIELD,PDEPTH)
342 !     #######################
343 !
344 USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE
345 USE MODD_SURF_PAR,       ONLY : XUNDEF
346 USE MODD_ISBA_PAR,    ONLY : XOPTIMGRID
347 !
348 USE MODI_OPEN_AUX_IO_SURF
349 USE MODI_CLOSE_AUX_IO_SURF
350 USE MODI_READ_SURF
351 USE MODE_SOIL
352 !
353 IMPLICIT NONE
354 !
355 !* dummy arguments
356 !  ---------------
357 !
358  CHARACTER(LEN=28),  INTENT(IN)  :: HFILE     ! name of file
359  CHARACTER(LEN=6),   INTENT(IN)  :: HFILETYPE ! type of input file
360  CHARACTER(LEN=28),  INTENT(IN)  :: HFILEPGD     ! name of file
361  CHARACTER(LEN=6),   INTENT(IN)  :: HFILEPGDTYPE ! type of input file
362 INTEGER,              INTENT(IN)    :: KLUOUT    ! logical unit of output listing
363 INTEGER,              INTENT(IN)    :: KNI       ! number of points
364  CHARACTER(LEN=7),     INTENT(IN)    :: HFIELD    ! field name
365  CHARACTER(LEN=*),     INTENT(IN)    :: HNAME     ! field name in the file
366 REAL, DIMENSION(:,:,:), POINTER       :: PFIELD    ! field to initialize
367 REAL, DIMENSION(:,:,:), POINTER       :: PDEPTH    ! middle depth of each layer
368 !
369 !
370 !* local variables
371 !  ---------------
372 !
373  CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
374  CHARACTER(LEN=4)  :: YLVL
375 #ifdef MNH_PARALLEL
376  CHARACTER(LEN=8)  :: YPATCH
377 #endif
378  CHARACTER(LEN=3)  :: YISBA          ! type of ISBA soil scheme
379  CHARACTER(LEN=3)  :: YNAT           ! type of surface (nature, garden)
380  CHARACTER(LEN=4)  :: YPEDOTF        ! type of pedo-transfert function
381 INTEGER           :: IRESP          ! reading return code
382 INTEGER           :: ILAYER         ! number of layers
383 INTEGER           :: JLAYER         ! loop counter
384 INTEGER           :: IPATCH         ! number of patch
385 INTEGER           :: JPATCH         ! loop counter
386 INTEGER           :: JVEGTYPE       ! loop counter
387 LOGICAL           :: GTEB           ! TEB field
388 !
389 REAL,  DIMENSION(:,:,:), ALLOCATABLE :: ZFIELD ! field read, one level, all patches
390 REAL,  DIMENSION(:,:),   ALLOCATABLE :: ZWORK  ! field read, one level, all patches
391 !
392 REAL,  DIMENSION(:,:,:), ALLOCATABLE :: ZVAR      ! profile of physical variable
393 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZCLAY     ! clay fraction
394 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZSAND     ! sand fraction
395 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZWWILT    ! wilting point
396 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZWFC      ! field capacity
397 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZWSAT     ! saturation
398 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZSOILGRID
399 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZNAT      ! natural surface fraction 
400 !
401 INTEGER :: IVERSION   ! surface version
402 INTEGER :: IBUGFIX
403 !
404 REAL(KIND=JPRB) :: ZHOOK_HANDLE
405 !-------------------------------------------------------------------------------
406 IF (LHOOK) CALL DR_HOOK('MODE_READ_EXTERN:READ_EXTERN_ISBA',0,ZHOOK_HANDLE)
407 WRITE  (KLUOUT,*) ' | Reading ',HFIELD,' in externalized file'
408 !
409 GTEB = (HNAME(1:3)=='TWN' .OR. HNAME(1:3)=='GD_' .OR. HNAME(1:3)=='GR_' &
410         .OR. HNAME(4:6)=='GD_' .OR. HNAME(4:6)=='GR_')
411 !
412 !------------------------------------------------------------------------------
413 !
414 IF (GTEB) THEN
415   CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'TOWN  ')
416 ELSE
417   CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'NATURE')
418 ENDIF
419 !
420 YRECFM='VERSION'
421  CALL READ_SURF(HFILEPGDTYPE,YRECFM,IVERSION,IRESP)
422 !
423 YRECFM='BUG'
424  CALL READ_SURF(HFILEPGDTYPE,YRECFM,IBUGFIX,IRESP)
425 !
426 !* Read number of soil layers
427 !
428 YRECFM='GROUND_LAYER'
429 IF (GTEB) THEN 
430   YRECFM='TWN_LAYER'
431   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_LAYER'
432 ENDIF
433  CALL READ_SURF(HFILEPGDTYPE,YRECFM,ILAYER,IRESP)
434 !
435 !* number of tiles
436 !
437 IPATCH=1
438 IF (.NOT. GTEB) THEN
439   YRECFM='PATCH_NUMBER'
440   CALL READ_SURF(HFILEPGDTYPE,YRECFM,IPATCH,IRESP)
441 END IF
442 !
443 !* soil scheme
444 !
445 YRECFM='ISBA'
446 IF (GTEB) THEN 
447   YRECFM='TWN_ISBA'
448   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_ISBA'
449 ENDIF
450  CALL READ_SURF(HFILEPGDTYPE,YRECFM,YISBA,IRESP)
451 !
452 IF (IVERSION>=7) THEN
453   !
454   !* Pedo-transfert function
455   !
456   YRECFM='PEDOTF'
457   IF (GTEB) THEN 
458     YRECFM='TWN_PEDOTF'
459     IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_PEDOTF'
460   ENDIF
461   CALL READ_SURF(HFILEPGDTYPE,YRECFM,YPEDOTF,IRESP)
462   !
463 ELSE
464   YPEDOTF = 'CH78'
465 ENDIF
466 !
467 !Only Brook and Corey with Force-Restore scheme
468 IF(YISBA/='DIF')THEN
469   YPEDOTF='CH78'
470 ENDIF
471 !
472 !-------------------------------------------------------------------------------
473 !
474 ! *.  Read clay fraction
475 !     ------------------
476 !
477 ALLOCATE(ZCLAY(KNI))
478 YRECFM='CLAY'
479 IF (GTEB) THEN 
480   YRECFM='TWN_CLAY'
481   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_CLAY'
482 ENDIF
483  CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZCLAY(:),IRESP,HDIR='A')
484 !
485 !-------------------------------------------------------------------------------
486 !
487 ! *.  Read sand fraction
488 !     ------------------
489 !
490 ALLOCATE(ZSAND(KNI))
491 YRECFM='SAND'
492 IF (GTEB) THEN 
493   YRECFM='TWN_SAND'
494   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_SAND'
495 ENDIF
496  CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZSAND(:),IRESP,HDIR='A')
497 !
498 !-------------------------------------------------------------------------------
499 !
500 ! *.  Read soil grid
501 !     --------------
502 !
503 !* Reference grid for DIF
504 !
505 IF(YISBA=='DIF') THEN
506   ALLOCATE(ZSOILGRID(ILAYER))
507   ZSOILGRID=XUNDEF
508   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=2) THEN
509     YRECFM='SOILGRID'
510     IF (GTEB) THEN 
511       YRECFM='TWN_SOILGRID'
512       IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_SOILGRID'
513     ENDIF
514     CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZSOILGRID,IRESP,HDIR='-')
515   ELSE
516     ZSOILGRID(1:ILAYER) = XOPTIMGRID(1:ILAYER)
517   ENDIF
518 ELSE
519   ALLOCATE(ZSOILGRID(0))
520 ENDIF
521 !
522 IF ((HFIELD=='TG    ') .AND. (YISBA=='2-L' .OR. YISBA=='3-L')) THEN
523   ALLOCATE(PDEPTH    (KNI,ILAYER,NVEGTYPE))
524   DO JVEGTYPE=1,NVEGTYPE
525     PDEPTH(:,1,JVEGTYPE) = 0.
526     PDEPTH(:,2,JVEGTYPE) = 0.2
527     IF (ILAYER==3) PDEPTH(:,3,JVEGTYPE) = 3.
528   END DO
529 ELSE
530   YNAT='NAT'
531   IF (GTEB) YNAT='GRD'
532   CALL READ_EXTERN_DEPTH(HFILEPGDTYPE,KLUOUT,YISBA,YNAT,HFIELD,KNI,ILAYER,IPATCH,&
533                          ZSOILGRID,PDEPTH,IVERSION)
534 END IF
535 !
536 DEALLOCATE(ZSOILGRID)
537 !
538 ! *.  Read fraction of nature
539 !     --------------
540 !
541 ALLOCATE(ZNAT(KNI))
542 IF (IVERSION>=7) THEN
543   CALL READ_SURF(HFILEPGDTYPE,'FRAC_NATURE',ZNAT,IRESP,HDIR='A')
544 ENDIF
545
546 !
547  CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
548 !
549 !* Allocate soil variable profile
550 !  ------------------------------
551 !
552 !
553 ALLOCATE(ZVAR(KNI,ILAYER,IPATCH))
554 ALLOCATE(ZWORK(KNI,IPATCH))
555 ZWORK(:,:) = XUNDEF
556 !
557 ! *.  Read soil variable profile
558 !     --------------------------
559 !
560 IF (GTEB) THEN
561   CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'TOWN  ')
562 ELSE
563   CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'NATURE')
564 ENDIF
565 !
566 DO JLAYER=1,ILAYER
567   WRITE(YLVL,'(I4)') JLAYER
568 #ifdef MNH_PARALLEL
569   DO JPATCH=1,IPATCH
570     IF (JLAYER >= 10) WRITE(YPATCH,'(I2,I4.4)') JLAYER,JPATCH
571     IF (JLAYER < 10)  WRITE(YPATCH,FMT='(I1,I4.4)') JLAYER,JPATCH
572     YRECFM=TRIM(HNAME)//ADJUSTL(YPATCH(:LEN_TRIM(YPATCH)))
573     CALL READ_SURF(HFILETYPE,YRECFM,ZWORK(:,JPATCH),IRESP,HDIR='A')
574     ZVAR(:,JLAYER,JPATCH)=ZWORK(:,JPATCH)
575   END DO
576 #else
577   YRECFM=TRIM(HNAME)//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
578   CALL READ_SURF(HFILETYPE,YRECFM,ZWORK(:,:),IRESP,HDIR='A')
579   DO JPATCH=1,IPATCH
580     ZVAR(:,JLAYER,JPATCH)=ZWORK(:,JPATCH)
581   END DO
582 #endif
583 END DO
584 !
585  CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
586 !
587 DEALLOCATE(ZWORK)
588 !
589 !
590 ! *.  Compute relative humidity from units kg/m^2 (SWI)
591 !     ------------------------------------------------
592 !
593 !* In case of force-restore ISBA, adds one layer at bottom of surface layer
594 IF ((HFIELD=='WG    ' .OR. HFIELD=='WGI   ') .AND. (YISBA=='2-L' .OR. YISBA=='3-L')) THEN
595   ALLOCATE(ZFIELD(KNI,ILAYER,IPATCH))
596   ZFIELD(:,:,:) = ZVAR(:,:,:)
597   DEALLOCATE(ZVAR)
598   !
599   ILAYER = ILAYER + 1
600   IF ( YISBA=='3-L' ) ILAYER = ILAYER + 1
601   ALLOCATE(ZVAR(KNI,ILAYER,IPATCH))
602   DO JPATCH=1,IPATCH
603     ZVAR(:,1,JPATCH)=ZFIELD(:,1,JPATCH)
604     ZVAR(:,2,JPATCH)=ZFIELD(:,2,JPATCH)  ! new layer at root layer humidity but below surface layer
605     ZVAR(:,3,JPATCH)=ZFIELD(:,2,JPATCH)
606     IF ( YISBA=='3-L' ) THEN
607       ZVAR(:,4,JPATCH)=ZFIELD(:,3,JPATCH)
608       ZVAR(:,5,JPATCH)=ZFIELD(:,3,JPATCH)
609     END IF
610   END DO
611   DEALLOCATE(ZFIELD)
612 END IF
613 !
614 ALLOCATE(ZFIELD(KNI,ILAYER,IPATCH))
615 ZFIELD = ZVAR
616 !
617 IF (HFIELD=='WG    ' .OR. HFIELD=='WGI   ') THEN
618   !
619   ! Compute ISBA model constants
620   !
621   ALLOCATE (ZWFC  (KNI))
622   ALLOCATE (ZWWILT(KNI))
623   ALLOCATE (ZWSAT (KNI))
624   !
625   ZWSAT (:) = WSAT_FUNC (ZCLAY(:),ZSAND(:),YPEDOTF)
626   ZWWILT(:) = WWILT_FUNC(ZCLAY(:),ZSAND(:),YPEDOTF)
627   ZWFC  (:) = WFC_FUNC  (ZCLAY(:),ZSAND(:),YPEDOTF)
628   !
629   DEALLOCATE (ZSAND)
630   DEALLOCATE (ZCLAY)
631
632   ZFIELD(:,:,:) = XUNDEF
633   !
634   IF (HFIELD=='WG    ') THEN
635     DO JPATCH=1,IPATCH
636       DO JLAYER=1,ILAYER
637         WHERE(ZNAT(:)>0.0)
638           ZVAR(:,JLAYER,JPATCH) = MAX(MIN(ZVAR(:,JLAYER,JPATCH),ZWSAT(:)),0.)
639           !
640           ZFIELD(:,JLAYER,JPATCH) = (ZVAR(:,JLAYER,JPATCH) - ZWWILT(:)) / (ZWFC(:) - ZWWILT(:))
641         END WHERE
642       END DO
643     END DO
644   ELSE IF (HFIELD=='WGI   ') THEN
645     DO JPATCH=1,IPATCH
646       DO JLAYER=1,ILAYER
647         WHERE(ZNAT(:)>0.0)
648           ZFIELD(:,JLAYER,JPATCH) = ZVAR(:,JLAYER,JPATCH) / ZWSAT(:)  
649         END WHERE 
650       END DO
651     END DO
652   END IF
653 !
654   DEALLOCATE (ZNAT)
655   DEALLOCATE (ZWSAT)
656   DEALLOCATE (ZWWILT)
657   DEALLOCATE (ZWFC)
658 !
659 !
660 END IF
661 !
662 DEALLOCATE(ZVAR)
663 !-------------------------------------------------------------------------------
664 !
665 ! *.  Set the field on all vegtypes
666 !     -----------------------------
667 !
668 ALLOCATE(PFIELD(KNI,ILAYER,NVEGTYPE))
669  CALL PUT_ON_ALL_VEGTYPES(KNI,ILAYER,IPATCH,NVEGTYPE,ZFIELD,PFIELD)
670 DEALLOCATE(ZFIELD)
671 IF (LHOOK) CALL DR_HOOK('MODE_READ_EXTERN:READ_EXTERN_ISBA',1,ZHOOK_HANDLE)
672 !
673 !------------------------------------------------------------------------------
674 !
675 END SUBROUTINE READ_EXTERN_ISBA
676 !
677 !------------------------------------------------------------------------------
678 !
679 END MODULE MODE_READ_EXTERN