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