Juan 10/07/2015: phasage MNH-514
[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 !
105 ALLOCATE(ZDG   (KNI,KLAYER,KPATCH))
106 ALLOCATE(IWG_LAYER   (KNI,KPATCH))
107 IWG_LAYER(:,:) = NUNDEF
108 IHYDRO_LAYER = KLAYER
109 !
110 IF (GECOCLIMAP) THEN
111   !
112   !* reading of the cover to obtain the depth of inter-layers
113   !
114   CALL OLD_NAME(HPROGRAM,'COVER_LIST      ',YRECFM)
115   CALL READ_SURF(HPROGRAM,YRECFM,GCOVER(:),IRESP,HDIR='-')
116   !
117   ALLOCATE(ZCOVER(KNI,JPCOVER))
118   YRECFM='COVER'
119   CALL READ_SURF(HPROGRAM,YRECFM,ZCOVER(:,:),GCOVER(:),IRESP,HDIR='A')  
120   !
121   !* computes soil layers
122   !  
123   CALL CONVERT_COVER_ISBA(HISBA,NUNDEF,ZCOVER,'   ',HNAT,PSOILGRID=PSOILGRID,PDG=ZDG,KWG_LAYER=IWG_LAYER)
124   IF (HISBA=='DIF') IHYDRO_LAYER = MAXVAL(IWG_LAYER(:,:),IWG_LAYER(:,:)/=NUNDEF)
125   !
126   DEALLOCATE(ZCOVER)
127   !
128 ENDIF
129 !
130 YRECFM='VERSION'
131  CALL READ_SURF(HPROGRAM,YRECFM,IVERSION,IRESP)
132 !
133 YRECFM='BUG'
134  CALL READ_SURF(HPROGRAM,YRECFM,IBUGFIX,IRESP)
135 !
136 !-------------------------------------------------------------------
137 IF (HNAT=='NAT' .AND. (IVERSION>=7 .OR. .NOT.GECOCLIMAP)) THEN
138   !
139   !* directly read soil layers in the file for nature ISBA soil layers
140   !
141   GDATA_DG = .TRUE.
142   IF (IVERSION>=7) THEN
143     YRECFM='L_DG'
144     YCOMMENT=YRECFM
145     CALL READ_SURF(HPROGRAM,YRECFM,GDATA_DG,IRESP,HCOMMENT=YCOMMENT)
146   ENDIF
147   !
148   IF (GDATA_DG) THEN
149     !
150     ALLOCATE(ZWORK(KNI,KPATCH))
151     DO JLAYER=1,KLAYER
152       IF (JLAYER<10)  WRITE(YRECFM,FMT='(A4,I1.1)') 'D_DG',JLAYER
153       IF (JLAYER>=10) WRITE(YRECFM,FMT='(A4,I2.2)') 'D_DG',JLAYER
154       CALL READ_SURF_ISBA_PAR_n(HPROGRAM,YRECFM,KLUOUT,KNI,ZWORK,IRESP,IVERSION,HDIR='A')
155       DO JPATCH=1,KPATCH
156         ZDG(:,JLAYER,JPATCH) = ZWORK(:,JPATCH)
157       END DO
158     END DO
159     DEALLOCATE(ZWORK)
160     !
161   ENDIF
162   !
163     GDATA_GROUND_DEPTH=.FALSE.
164   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=2) THEN
165     !
166     YRECFM2='L_GROUND_DEPTH'
167     IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM2='L_GROUND_DPT'
168     YCOMMENT=YRECFM2
169     CALL READ_SURF(HPROGRAM,YRECFM2,GDATA_GROUND_DEPTH,IRESP,HCOMMENT=YCOMMENT)
170     !
171     IF (GDATA_GROUND_DEPTH) THEN
172       !
173       YRECFM2='D_GROUND_DETPH'
174       IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM2='D_GROUND_DPT'
175       ALLOCATE(ZGROUND_DEPTH(KNI,KPATCH))
176       CALL READ_SURF_ISBA_PAR_n(HPROGRAM,YRECFM2,KLUOUT,KNI,ZGROUND_DEPTH(:,:),IRESP,IVERSION,HDIR='A')
177       !
178       DO JPATCH=1,KPATCH
179         DO JJ=1,KNI
180           DO JLAYER=1,KLAYER
181             IF ( ZDG(JJ,JLAYER,JPATCH) <= ZGROUND_DEPTH(JJ,JPATCH) .AND. ZGROUND_DEPTH(JJ,JPATCH) < XUNDEF ) &
182                 IWG_LAYER(JJ,JPATCH) = JLAYER
183           ENDDO
184         ENDDO
185       ENDDO
186       DEALLOCATE(ZGROUND_DEPTH)
187       !
188       IF (HISBA=='DIF') IHYDRO_LAYER = MAXVAL(IWG_LAYER(:,:),IWG_LAYER(:,:)/=NUNDEF)
189       !
190     ENDIF
191     !
192   ENDIF
193   !
194 ELSE IF (HNAT=='GRD' .AND. .NOT.GECOCLIMAP) THEN
195   !
196   !* computes soil layers from vegetation fractions read in the file
197   !
198   CALL READ_SURF(HPROGRAM,'D_TYPE_HVEG',YHVEG,IRESP)
199   CALL READ_SURF(HPROGRAM,'D_TYPE_LVEG',YLVEG,IRESP)
200   CALL READ_SURF(HPROGRAM,'D_TYPE_NVEG',YNVEG,IRESP)
201   CALL READ_SURF(HPROGRAM,'D_FRAC_HVEG',ZHVEG,IRESP,HDIR='A')
202   CALL READ_SURF(HPROGRAM,'D_FRAC_LVEG',ZLVEG,IRESP,HDIR='A')
203   CALL READ_SURF(HPROGRAM,'D_FRAC_NVEG',ZNVEG,IRESP,HDIR='A')
204   ! Ground layers
205   CALL GARDEN_SOIL_DEPTH(YNVEG,YLVEG,YHVEG,ZNVEG,ZLVEG,ZHVEG,ZDG)
206   !
207 END IF
208 !
209 DEALLOCATE(IWG_LAYER)
210 !
211 IF (HFIELD=='WG    ' .OR. HFIELD=='WGI   ' .OR. HFIELD=='TWN_WG  ' .OR. HFIELD=='TWN_WGI ' .OR. &
212       HFIELD=='GD_WG  ' .OR. HFIELD=='GD_WGI ') THEN
213   KLAYER = IHYDRO_LAYER
214 ENDIF
215 !
216 !-------------------------------------------------------------------
217 !
218 !* In force-restore ISBA, adds a layer at bottom of surface layer and a layer
219 !  between root and deep layers.
220 !
221 IF (HISBA=='2-L' .OR. HISBA=='3-L') THEN
222   ILAYER = KLAYER + 1
223   IF (HISBA=='3-L') ILAYER = ILAYER + 1
224   ALLOCATE(ZD    (KNI,ILAYER,KPATCH))
225   DO JPATCH=1,KPATCH
226     ! for interpolations, middle of surface layer must be at least at 1cm
227     ZD(:,1,JPATCH) = MIN(3.*ZDG(:,1,JPATCH),MAX(ZDG(:,1,JPATCH),0.02))
228     ! new layer below surface layer. This layer will be at root depth layer humidity
229     ZD(:,2,JPATCH) = MIN(4.*ZDG(:,1,JPATCH),0.5*(ZDG(:,1,JPATCH)+ZDG(:,2,JPATCH)))
230     ! root layer
231     ZD(:,3,JPATCH) = ZDG(:,2,JPATCH)
232     IF (HISBA=='3-L') THEN
233       ! between root and deep layers. This layer will have deep soil humidity.
234       WHERE (ZDG(:,2,JPATCH)<ZDG(:,3,JPATCH))
235         ZD(:,4,JPATCH) = 0.75 * ZDG(:,2,JPATCH) + 0.25 * ZDG(:,3,JPATCH)
236       ELSEWHERE
237         ZD(:,4,JPATCH) = ZDG(:,3,JPATCH)
238       END WHERE
239       ! deep layer
240       ZD(:,5,JPATCH) = ZDG(:,3,JPATCH)
241     END IF
242   END DO
243 ELSE
244   ILAYER = KLAYER
245   ALLOCATE(ZD    (KNI,ILAYER,KPATCH))
246   ZD(:,:,:) = ZDG(:,1:KLAYER,:)
247 END IF
248 !
249 DEALLOCATE(ZDG)
250 !
251 !-------------------------------------------------------------------
252 !* recovers middle layer depth (from the surface)
253 ALLOCATE(ZDEPTH    (KNI,ILAYER,KPATCH))
254 ZDEPTH = XUNDEF
255 DO JPATCH=1,KPATCH
256   WHERE(ZD(:,1,JPATCH)/=XUNDEF) &
257     ZDEPTH    (:,1,JPATCH)=ZD(:,1,JPATCH)/2.  
258   DO JLAYER=2,ILAYER
259     WHERE(ZD(:,1,JPATCH)/=XUNDEF) &
260       ZDEPTH    (:,JLAYER,JPATCH) = (ZD(:,JLAYER-1,JPATCH) + ZD(:,JLAYER,JPATCH))/2.  
261   END DO
262 END DO
263 DEALLOCATE(ZD)
264 !-------------------------------------------------------------------
265 !
266 ALLOCATE(PDEPTH    (KNI,ILAYER,NVEGTYPE))
267  CALL PUT_ON_ALL_VEGTYPES(KNI,ILAYER,KPATCH,NVEGTYPE,ZDEPTH,PDEPTH)
268 DEALLOCATE(ZDEPTH)
269
270 IF (LHOOK) CALL DR_HOOK('MODE_READ_EXTERN:READ_EXTERN_DEPTH',1,ZHOOK_HANDLE)
271 !-------------------------------------------------------------------
272 !
273 END SUBROUTINE READ_EXTERN_DEPTH
274 !
275 !
276 !-------------------------------------------------------------------
277 !---------------------------------------------------------------------------------------
278 !
279 !     #######################
280       SUBROUTINE READ_EXTERN_ISBA(HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,&
281                                   KLUOUT,KNI,HFIELD,HNAME,PFIELD,PDEPTH)
282 !     #######################
283 !
284 USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE
285 USE MODD_SURF_PAR,       ONLY : XUNDEF
286 USE MODD_ISBA_PAR,    ONLY : XOPTIMGRID
287 !
288 USE MODI_OPEN_AUX_IO_SURF
289 USE MODI_CLOSE_AUX_IO_SURF
290 USE MODI_READ_SURF
291 USE MODE_SOIL
292 !
293 IMPLICIT NONE
294 !
295 !* dummy arguments
296 !  ---------------
297 !
298  CHARACTER(LEN=28),  INTENT(IN)  :: HFILE     ! name of file
299  CHARACTER(LEN=6),   INTENT(IN)  :: HFILETYPE ! type of input file
300  CHARACTER(LEN=28),  INTENT(IN)  :: HFILEPGD     ! name of file
301  CHARACTER(LEN=6),   INTENT(IN)  :: HFILEPGDTYPE ! type of input file
302 INTEGER,              INTENT(IN)    :: KLUOUT    ! logical unit of output listing
303 INTEGER,              INTENT(IN)    :: KNI       ! number of points
304  CHARACTER(LEN=7),     INTENT(IN)    :: HFIELD    ! field name
305  CHARACTER(LEN=*),     INTENT(IN)    :: HNAME     ! field name in the file
306 REAL, DIMENSION(:,:,:), POINTER       :: PFIELD    ! field to initialize
307 REAL, DIMENSION(:,:,:), POINTER       :: PDEPTH    ! middle depth of each layer
308 !
309 !
310 !* local variables
311 !  ---------------
312 !
313  CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
314  CHARACTER(LEN=4)  :: YLVL
315  CHARACTER(LEN=3)  :: YISBA          ! type of ISBA soil scheme
316  CHARACTER(LEN=3)  :: YNAT           ! type of surface (nature, garden)
317  CHARACTER(LEN=4)  :: YPEDOTF        ! type of pedo-transfert function
318 INTEGER           :: IRESP          ! reading return code
319 INTEGER           :: ILAYER         ! number of layers
320 INTEGER           :: JLAYER         ! loop counter
321 INTEGER           :: IPATCH         ! number of patch
322 INTEGER           :: JPATCH         ! loop counter
323 INTEGER           :: JVEGTYPE       ! loop counter
324 LOGICAL           :: GTEB           ! TEB field
325 !
326 REAL,  DIMENSION(:,:,:), ALLOCATABLE :: ZFIELD ! field read, one level, all patches
327 REAL,  DIMENSION(:,:),   ALLOCATABLE :: ZWORK  ! field read, one level, all patches
328 !
329 REAL,  DIMENSION(:,:,:), ALLOCATABLE :: ZVAR      ! profile of physical variable
330 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZCLAY     ! clay fraction
331 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZSAND     ! sand fraction
332 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZWWILT    ! wilting point
333 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZWFC      ! field capacity
334 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZWSAT     ! saturation
335 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZSOILGRID
336 REAL,  DIMENSION(:),   ALLOCATABLE   :: ZNAT      ! natural surface fraction 
337 !
338 INTEGER :: IVERSION   ! surface version
339 INTEGER :: IBUGFIX
340 !
341 REAL(KIND=JPRB) :: ZHOOK_HANDLE
342 !-------------------------------------------------------------------------------
343 IF (LHOOK) CALL DR_HOOK('MODE_READ_EXTERN:READ_EXTERN_ISBA',0,ZHOOK_HANDLE)
344 WRITE  (KLUOUT,*) ' | Reading ',HFIELD,' in externalized file'
345 !
346 GTEB = (HNAME(1:3)=='TWN' .OR. HNAME(1:3)=='GD_' .OR. HNAME(1:3)=='GR_' &
347         .OR. HNAME(4:6)=='GD_' .OR. HNAME(4:6)=='GR_')
348 !
349 !------------------------------------------------------------------------------
350 !
351 IF (GTEB) THEN
352   CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'TOWN  ')
353 ELSE
354   CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'NATURE')
355 ENDIF
356 !
357 YRECFM='VERSION'
358  CALL READ_SURF(HFILEPGDTYPE,YRECFM,IVERSION,IRESP)
359 !
360 YRECFM='BUG'
361  CALL READ_SURF(HFILEPGDTYPE,YRECFM,IBUGFIX,IRESP)
362 !
363 !* Read number of soil layers
364 !
365 YRECFM='GROUND_LAYER'
366 IF (GTEB) THEN 
367   YRECFM='TWN_LAYER'
368   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_LAYER'
369 ENDIF
370  CALL READ_SURF(HFILEPGDTYPE,YRECFM,ILAYER,IRESP)
371 !
372 !* number of tiles
373 !
374 IPATCH=1
375 IF (.NOT. GTEB) THEN
376   YRECFM='PATCH_NUMBER'
377   CALL READ_SURF(HFILEPGDTYPE,YRECFM,IPATCH,IRESP)
378 END IF
379 !
380 !* soil scheme
381 !
382 YRECFM='ISBA'
383 IF (GTEB) THEN 
384   YRECFM='TWN_ISBA'
385   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_ISBA'
386 ENDIF
387  CALL READ_SURF(HFILEPGDTYPE,YRECFM,YISBA,IRESP)
388 !
389 IF (IVERSION>=7) THEN
390   !
391   !* Pedo-transfert function
392   !
393   YRECFM='PEDOTF'
394   IF (GTEB) THEN 
395     YRECFM='TWN_PEDOTF'
396     IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_PEDOTF'
397   ENDIF
398   CALL READ_SURF(HFILEPGDTYPE,YRECFM,YPEDOTF,IRESP)
399   !
400 ELSE
401   YPEDOTF = 'CH78'
402 ENDIF
403 !
404 !Only Brook and Corey with Force-Restore scheme
405 IF(YISBA/='DIF')THEN
406   YPEDOTF='CH78'
407 ENDIF
408 !
409 !-------------------------------------------------------------------------------
410 !
411 ! *.  Read clay fraction
412 !     ------------------
413 !
414 ALLOCATE(ZCLAY(KNI))
415 YRECFM='CLAY'
416 IF (GTEB) THEN 
417   YRECFM='TWN_CLAY'
418   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_CLAY'
419 ENDIF
420  CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZCLAY(:),IRESP,HDIR='A')
421 !
422 !-------------------------------------------------------------------------------
423 !
424 ! *.  Read sand fraction
425 !     ------------------
426 !
427 ALLOCATE(ZSAND(KNI))
428 YRECFM='SAND'
429 IF (GTEB) THEN 
430   YRECFM='TWN_SAND'
431   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_SAND'
432 ENDIF
433  CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZSAND(:),IRESP,HDIR='A')
434 !
435 !-------------------------------------------------------------------------------
436 !
437 ! *.  Read soil grid
438 !     --------------
439 !
440 !* Reference grid for DIF
441 !
442 IF(YISBA=='DIF') THEN
443   ALLOCATE(ZSOILGRID(ILAYER))
444   ZSOILGRID=XUNDEF
445   IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=2) THEN
446     YRECFM='SOILGRID'
447     IF (GTEB) THEN 
448       YRECFM='TWN_SOILGRID'
449       IF (IVERSION>7 .OR. IVERSION==7 .AND. IBUGFIX>=3) YRECFM='GD_SOILGRID'
450     ENDIF
451     CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZSOILGRID,IRESP,HDIR='-')
452   ELSE
453     ZSOILGRID(1:ILAYER) = XOPTIMGRID(1:ILAYER)
454   ENDIF
455 ELSE
456   ALLOCATE(ZSOILGRID(0))
457 ENDIF
458 !
459 IF ((HFIELD=='TG    ') .AND. (YISBA=='2-L' .OR. YISBA=='3-L')) THEN
460   ALLOCATE(PDEPTH    (KNI,ILAYER,NVEGTYPE))
461   DO JVEGTYPE=1,NVEGTYPE
462     PDEPTH(:,1,JVEGTYPE) = 0.
463     PDEPTH(:,2,JVEGTYPE) = 0.2
464     IF (ILAYER==3) PDEPTH(:,3,JVEGTYPE) = 3.
465   END DO
466 ELSE
467   YNAT='NAT'
468   IF (GTEB) YNAT='GRD'
469   CALL READ_EXTERN_DEPTH(HFILEPGDTYPE,KLUOUT,YISBA,YNAT,HFIELD,KNI,ILAYER,IPATCH,&
470                          ZSOILGRID,PDEPTH,IVERSION)
471 END IF
472 !
473 DEALLOCATE(ZSOILGRID)
474 !
475 ! *.  Read fraction of nature
476 !     --------------
477 !
478 ALLOCATE(ZNAT(KNI))
479 IF (IVERSION>=7) THEN
480   CALL READ_SURF(HFILEPGDTYPE,'FRAC_NATURE',ZNAT,IRESP,HDIR='A')
481 ENDIF
482
483 !
484  CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
485 !
486 !* Allocate soil variable profile
487 !  ------------------------------
488 !
489 !
490 ALLOCATE(ZVAR(KNI,ILAYER,IPATCH))
491 ALLOCATE(ZWORK(KNI,IPATCH))
492 ZWORK(:,:) = XUNDEF
493 !
494 ! *.  Read soil variable profile
495 !     --------------------------
496 !
497 IF (GTEB) THEN
498   CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'TOWN  ')
499 ELSE
500   CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'NATURE')
501 ENDIF
502 !
503 DO JLAYER=1,ILAYER
504   WRITE(YLVL,'(I4)') JLAYER
505   YRECFM=TRIM(HNAME)//ADJUSTL(YLVL(:LEN_TRIM(YLVL)))
506   CALL READ_SURF(HFILETYPE,YRECFM,ZWORK(:,:),IRESP,HDIR='A')
507   DO JPATCH=1,IPATCH
508     ZVAR(:,JLAYER,JPATCH)=ZWORK(:,JPATCH)
509   END DO
510 END DO
511 !
512  CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
513 !
514 DEALLOCATE(ZWORK)
515 !
516 !
517 ! *.  Compute relative humidity from units kg/m^2 (SWI)
518 !     ------------------------------------------------
519 !
520 !* In case of force-restore ISBA, adds one layer at bottom of surface layer
521 IF ((HFIELD=='WG    ' .OR. HFIELD=='WGI   ') .AND. (YISBA=='2-L' .OR. YISBA=='3-L')) THEN
522   ALLOCATE(ZFIELD(KNI,ILAYER,IPATCH))
523   ZFIELD(:,:,:) = ZVAR(:,:,:)
524   DEALLOCATE(ZVAR)
525   !
526   ILAYER = ILAYER + 1
527   IF ( YISBA=='3-L' ) ILAYER = ILAYER + 1
528   ALLOCATE(ZVAR(KNI,ILAYER,IPATCH))
529   DO JPATCH=1,IPATCH
530     ZVAR(:,1,JPATCH)=ZFIELD(:,1,JPATCH)
531     ZVAR(:,2,JPATCH)=ZFIELD(:,2,JPATCH)  ! new layer at root layer humidity but below surface layer
532     ZVAR(:,3,JPATCH)=ZFIELD(:,2,JPATCH)
533     IF ( YISBA=='3-L' ) THEN
534       ZVAR(:,4,JPATCH)=ZFIELD(:,3,JPATCH)
535       ZVAR(:,5,JPATCH)=ZFIELD(:,3,JPATCH)
536     END IF
537   END DO
538   DEALLOCATE(ZFIELD)
539 END IF
540 !
541 ALLOCATE(ZFIELD(KNI,ILAYER,IPATCH))
542 ZFIELD = ZVAR
543 !
544 IF (HFIELD=='WG    ' .OR. HFIELD=='WGI   ') THEN
545   !
546   ! Compute ISBA model constants
547   !
548   ALLOCATE (ZWFC  (KNI))
549   ALLOCATE (ZWWILT(KNI))
550   ALLOCATE (ZWSAT (KNI))
551   !
552   ZWSAT (:) = WSAT_FUNC (ZCLAY(:),ZSAND(:),YPEDOTF)
553   ZWWILT(:) = WWILT_FUNC(ZCLAY(:),ZSAND(:),YPEDOTF)
554   ZWFC  (:) = WFC_FUNC  (ZCLAY(:),ZSAND(:),YPEDOTF)
555   !
556   DEALLOCATE (ZSAND)
557   DEALLOCATE (ZCLAY)
558
559   ZFIELD(:,:,:) = XUNDEF
560   !
561   IF (HFIELD=='WG    ') THEN
562     DO JPATCH=1,IPATCH
563       DO JLAYER=1,ILAYER
564         WHERE(ZNAT(:)>0.0)
565           ZVAR(:,JLAYER,JPATCH) = MAX(MIN(ZVAR(:,JLAYER,JPATCH),ZWSAT(:)),0.)
566           !
567           ZFIELD(:,JLAYER,JPATCH) = (ZVAR(:,JLAYER,JPATCH) - ZWWILT(:)) / (ZWFC(:) - ZWWILT(:))
568         END WHERE
569       END DO
570     END DO
571   ELSE IF (HFIELD=='WGI   ') THEN
572     DO JPATCH=1,IPATCH
573       DO JLAYER=1,ILAYER
574         WHERE(ZNAT(:)>0.0)
575           ZFIELD(:,JLAYER,JPATCH) = ZVAR(:,JLAYER,JPATCH) / ZWSAT(:)  
576         END WHERE 
577       END DO
578     END DO
579   END IF
580 !
581   DEALLOCATE (ZNAT)
582   DEALLOCATE (ZWSAT)
583   DEALLOCATE (ZWWILT)
584   DEALLOCATE (ZWFC)
585 !
586 !
587 END IF
588 !
589 DEALLOCATE(ZVAR)
590 !-------------------------------------------------------------------------------
591 !
592 ! *.  Set the field on all vegtypes
593 !     -----------------------------
594 !
595 ALLOCATE(PFIELD(KNI,ILAYER,NVEGTYPE))
596  CALL PUT_ON_ALL_VEGTYPES(KNI,ILAYER,IPATCH,NVEGTYPE,ZFIELD,PFIELD)
597 DEALLOCATE(ZFIELD)
598 IF (LHOOK) CALL DR_HOOK('MODE_READ_EXTERN:READ_EXTERN_ISBA',1,ZHOOK_HANDLE)
599 !
600 !------------------------------------------------------------------------------
601 !
602 END SUBROUTINE READ_EXTERN_ISBA
603 !
604 !------------------------------------------------------------------------------
605 !
606 END MODULE MODE_READ_EXTERN