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