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