Juan 13/01/2014: add header SURFEX_LIC to all SURFEX files
[MNH-git_open_source-lfs.git] / src / SURFEX / pgd_isba.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       SUBROUTINE PGD_ISBA(HPROGRAM,OECOCLIMAP)
7 !     ##############################################################
8 !
9 !!**** *PGD_ISBA* monitor for averaging and interpolations of ISBA physiographic fields
10 !!
11 !!    PURPOSE
12 !!    -------
13 !!
14 !!    METHOD
15 !!    ------
16 !!   
17 !
18 !!    EXTERNAL
19 !!    --------
20 !!
21 !!    IMPLICIT ARGUMENTS
22 !!    ------------------
23 !!
24 !!    REFERENCE
25 !!    ---------
26 !!
27 !!    AUTHOR
28 !!    ------
29 !!
30 !!    V. Masson        Meteo-France
31 !!
32 !!    MODIFICATION
33 !!    ------------
34 !!
35 !!    Original    10/12/97
36 !!    P. Le Moigne  12/2004 : add type of photosynthesis and correct computation
37 !!                            of ground layers number in diffusion case
38 !!    P. Le Moigne  09/2005 : AGS modifs of L. Jarlan
39 !!    B. Decharme      2008 :  XWDRAIN
40 !!    E. Martin     12/2008 : files of data for runoffb and wdrain
41 !!    B. Decharme   06/2009 : files of data for topographic index
42 !!    A.L. Gibelin  04/2009 : dimension NBIOMASS for ISBA-A-gs
43 !!
44 !----------------------------------------------------------------------------
45 !
46 !*    0.     DECLARATION
47 !            -----------
48 !
49 USE MODD_SURF_PAR,       ONLY : XUNDEF, NUNDEF
50 USE MODD_PGD_GRID,       ONLY : NL
51 USE MODD_PGDWORK,        ONLY : CATYPE
52 USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE, JPCOVER
53 USE MODD_ISBA_n,         ONLY : NPATCH, NGROUND_LAYER, NNBIOMASS, CISBA, &
54                                 CPEDOTF, XCOVER, LCOVER, XZS,            &
55                                 XZ0EFFJPDIR, CPHOTO, LTR_ML, XRM_PATCH,  &
56                                 XCLAY, XSAND, XSOC, LSOCP, LNOF,         &
57                                 XRUNOFFB, XWDRAIN, LECOCLIMAP,           &
58                                 XSOILGRID, LPERM, XPERM, XPH, XFERT 
59 USE MODD_ISBA_GRID_n,    ONLY : CGRID, XGRID_PAR, XLAT, XLON, XMESH_SIZE
60 !
61 USE MODD_ISBA_PAR,       ONLY : NOPTIMLAYER, XOPTIMGRID
62 !
63 USE MODI_GET_LUOUT
64 USE MODI_READ_NAM_PGD_ISBA
65 USE MODI_PGD_FIELD
66 USE MODI_TEST_NAM_VAR_SURF
67 !
68 USE MODI_GET_AOS_n
69 USE MODI_GET_SSO_n
70 USE MODI_GET_SURF_SIZE_n
71 USE MODI_PACK_PGD_ISBA
72 USE MODI_PACK_PGD
73 USE MODI_WRITE_COVER_TEX_ISBA
74 USE MODI_WRITE_COVER_TEX_ISBA_PAR
75 USE MODI_PGD_TOPO_INDEX
76 USE MODI_PGD_ISBA_PAR
77 USE MODI_PGD_TOPD
78 !
79 USE MODI_READ_SURF
80 USE MODI_INIT_IO_SURF_n
81 USE MODI_END_IO_SURF_n
82 #ifdef ASC
83 USE MODD_IO_SURF_ASC, ONLY : CFILEIN
84 #endif
85 #ifdef FA
86 USE MODD_IO_SURF_FA,  ONLY : CFILEIN_FA
87 #endif
88 #ifdef LFI
89 USE MODD_IO_SURF_LFI, ONLY : CFILEIN_LFI
90 #endif
91 !
92 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
93 USE PARKIND1  ,ONLY : JPRB
94 !
95 USE MODI_ABOR1_SFX
96 !
97 IMPLICIT NONE
98 !
99 !*    0.1    Declaration of arguments
100 !            ------------------------
101 !
102  CHARACTER(LEN=6), INTENT(IN)  :: HPROGRAM   ! program calling surf. schemes
103 LOGICAL,          INTENT(IN)  :: OECOCLIMAP ! T if parameters are computed with ecoclimap
104 !                                           ! F if all parameters must be specified
105 !
106 !
107 !*    0.2    Declaration of local variables
108 !            ------------------------------
109 !
110 INTEGER                           :: ILUOUT    ! output listing logical unit
111 INTEGER                           :: JLAYER    ! loop counter
112 INTEGER                           :: ILU       ! number of points
113 REAL, DIMENSION(NL)               :: ZAOSIP    ! A/S i+ on all surface points
114 REAL, DIMENSION(NL)               :: ZAOSIM    ! A/S i- on all surface points
115 REAL, DIMENSION(NL)               :: ZAOSJP    ! A/S j+ on all surface points
116 REAL, DIMENSION(NL)               :: ZAOSJM    ! A/S j- on all surface points
117 REAL, DIMENSION(NL)               :: ZHO2IP    ! h/2 i+ on all surface points
118 REAL, DIMENSION(NL)               :: ZHO2IM    ! h/2 i- on all surface points
119 REAL, DIMENSION(NL)               :: ZHO2JP    ! h/2 j+ on all surface points
120 REAL, DIMENSION(NL)               :: ZHO2JM    ! h/2 j- on all surface points
121 REAL, DIMENSION(NL)               :: ZSSO_SLOPE! subgrid slope on all surface points
122 INTEGER                           :: IRESP     ! error code
123 !
124 !*    0.3    Declaration of namelists
125 !            ------------------------
126 !
127 !
128 INTEGER                  :: IPATCH           ! number of patches
129 INTEGER                  :: IGROUND_LAYER    ! number of soil layers
130  CHARACTER(LEN=3)         :: YISBA            ! ISBA option
131  CHARACTER(LEN=4)         :: YPEDOTF          ! Pedo transfert function for DIF
132  CHARACTER(LEN=3)         :: YPHOTO           ! photosynthesis option
133 LOGICAL                  :: GTR_ML           ! new radiative transfert
134 REAL                     :: ZRM_PATCH        ! threshold to remove little fractions of patches
135  CHARACTER(LEN=28)        :: YSAND            ! file name for sand fraction
136  CHARACTER(LEN=28)        :: YCLAY            ! file name for clay fraction
137  CHARACTER(LEN=28)        :: YSOC_TOP         ! file name for organic carbon top soil
138  CHARACTER(LEN=28)        :: YSOC_SUB         ! file name for organic carbon sub soil
139  CHARACTER(LEN=28)        :: YCTI             ! file name for topographic index
140  CHARACTER(LEN=28)        :: YRUNOFFB         ! file name for runoffb parameter
141  CHARACTER(LEN=28)        :: YWDRAIN          ! file name for wdrain parameter
142  CHARACTER(LEN=28)        :: YPERM            ! file name for permafrost distribution
143  CHARACTER(LEN=6)         :: YSANDFILETYPE    ! sand data file type
144  CHARACTER(LEN=6)         :: YCLAYFILETYPE    ! clay data file type
145  CHARACTER(LEN=6)         :: YSOCFILETYPE     ! organic carbon data file type
146  CHARACTER(LEN=6)         :: YCTIFILETYPE     ! topographic index data file type
147  CHARACTER(LEN=6)         :: YRUNOFFBFILETYPE ! subgrid runoff data file type
148  CHARACTER(LEN=6)         :: YWDRAINFILETYPE  ! subgrid drainage data file type
149  CHARACTER(LEN=6)         :: YPERMFILETYPE    ! permafrost distribution data file type
150 REAL                     :: XUNIF_SAND       ! uniform value of sand fraction  (-)
151 REAL                     :: XUNIF_CLAY       ! uniform value of clay fraction  (-)
152 REAL                     :: XUNIF_SOC_TOP    ! uniform value of organic carbon top soil (kg/m2)
153 REAL                     :: XUNIF_SOC_SUB    ! uniform value of organic carbon sub soil (kg/m2)
154 REAL                     :: XUNIF_RUNOFFB    ! uniform value of subgrid runoff coefficient
155 REAL                     :: XUNIF_WDRAIN     ! uniform subgrid drainage parameter
156 REAL                     :: XUNIF_PERM       ! uniform permafrost distribution
157 LOGICAL                  :: LIMP_SAND        ! Imposed maps of Sand
158 LOGICAL                  :: LIMP_CLAY        ! Imposed maps of Clay
159 LOGICAL                  :: LIMP_SOC         ! Imposed maps of organic carbon
160 LOGICAL                  :: LIMP_CTI         ! Imposed maps of topographic index statistics
161 LOGICAL                  :: LIMP_PERM        ! Imposed maps of permafrost distribution
162 REAL, DIMENSION(150)     :: ZSOILGRID        ! Soil grid reference for DIF
163  CHARACTER(LEN=28)        :: YPH           ! file name for pH
164  CHARACTER(LEN=28)        :: YFERT         ! file name for fertilisation rate
165  CHARACTER(LEN=6)         :: YPHFILETYPE   ! pH data file type
166  CHARACTER(LEN=6)         :: YFERTFILETYPE ! fertilisation data file type
167 REAL                     :: XUNIF_PH      ! uniform value of pH
168 REAL                     :: XUNIF_FERT    ! uniform value of fertilisation rate
169 !
170 REAL(KIND=JPRB) :: ZHOOK_HANDLE
171 !
172 !-------------------------------------------------------------------------------
173 !
174 IF (LHOOK) CALL DR_HOOK('PGD_ISBA',0,ZHOOK_HANDLE)
175  CALL GET_LUOUT(HPROGRAM,ILUOUT)
176 !
177 !-------------------------------------------------------------------------------
178 !
179 !*    2.      Reading of namelist
180 !             -------------------
181 !
182  CALL READ_NAM_PGD_ISBA(HPROGRAM, IPATCH, IGROUND_LAYER,                          &
183                        YISBA,  YPEDOTF, YPHOTO, GTR_ML, ZRM_PATCH,               &
184                        YCLAY, YCLAYFILETYPE, XUNIF_CLAY, LIMP_CLAY,              &
185                        YSAND, YSANDFILETYPE, XUNIF_SAND, LIMP_SAND,              &
186                        YSOC_TOP, YSOC_SUB, YSOCFILETYPE, XUNIF_SOC_TOP,          &
187                        XUNIF_SOC_SUB, LIMP_SOC, YCTI, YCTIFILETYPE, LIMP_CTI,    &
188                        YPERM, YPERMFILETYPE, XUNIF_PERM, LIMP_PERM,              &                       
189                        YRUNOFFB, YRUNOFFBFILETYPE, XUNIF_RUNOFFB,                &
190                        YWDRAIN,  YWDRAINFILETYPE , XUNIF_WDRAIN, ZSOILGRID,      &
191                        YPH, YPHFILETYPE, XUNIF_PH, YFERT, YFERTFILETYPE,         &
192                        XUNIF_FERT                          )  
193 !
194 NPATCH        = IPATCH
195 NGROUND_LAYER = IGROUND_LAYER
196 CISBA         = YISBA
197 CPEDOTF       = YPEDOTF
198 CPHOTO        = YPHOTO
199 LTR_ML        = GTR_ML
200 XRM_PATCH     = MAX(MIN(ZRM_PATCH,1.),0.)
201 !
202 !-------------------------------------------------------------------------------
203 !
204 !*    3.      Coherence of options
205 !             --------------------
206 !
207  CALL TEST_NAM_VAR_SURF(ILUOUT,'CISBA',CISBA,'2-L','3-L','DIF')
208  CALL TEST_NAM_VAR_SURF(ILUOUT,'CPEDOTF',CPEDOTF,'CH78','CO84')
209  CALL TEST_NAM_VAR_SURF(ILUOUT,'CPHOTO',CPHOTO,'NON','AGS','LAI','AST','LST','NIT','NCB')
210 !
211 SELECT CASE (CISBA)
212 !
213   CASE ('2-L')
214 !          
215     NGROUND_LAYER = 2
216     CPEDOTF       ='CH78'   
217     WRITE(ILUOUT,*) '*****************************************'
218     WRITE(ILUOUT,*) '* With option CISBA = ',CISBA,'         *'
219     WRITE(ILUOUT,*) '* the number of soil layers is set to 2 *'
220     WRITE(ILUOUT,*) '* Pedo transfert function = CH78        *'    
221     WRITE(ILUOUT,*) '*****************************************'
222 !    
223   CASE ('3-L')
224 !          
225     NGROUND_LAYER = 3
226     CPEDOTF       ='CH78'    
227     WRITE(ILUOUT,*) '*****************************************'
228     WRITE(ILUOUT,*) '* With option CISBA = ',CISBA,'         *'
229     WRITE(ILUOUT,*) '* the number of soil layers is set to 3 *'
230     WRITE(ILUOUT,*) '* Pedo transfert function = CH78        *'    
231     WRITE(ILUOUT,*) '*****************************************'
232 !    
233   CASE ('DIF')
234 !          
235     IF(NGROUND_LAYER==NUNDEF)THEN
236       IF(OECOCLIMAP)THEN
237         NGROUND_LAYER=NOPTIMLAYER
238       ELSE
239         WRITE(ILUOUT,*) '****************************************'
240         WRITE(ILUOUT,*) '* Number of ground layer not specified *'
241         WRITE(ILUOUT,*) '****************************************'
242         CALL ABOR1_SFX('PGD_ISBA: NGROUND_LAYER MUST BE DONE IN NAM_ISBA')
243       ENDIF
244     ENDIF
245
246     ALLOCATE(XSOILGRID(NGROUND_LAYER))
247     XSOILGRID(:)=XUNDEF
248     XSOILGRID(:)=ZSOILGRID(1:NGROUND_LAYER) 
249     IF (ALL(ZSOILGRID(:)==XUNDEF)) THEN
250       IF(OECOCLIMAP) XSOILGRID(1:NGROUND_LAYER)=XOPTIMGRID(1:NGROUND_LAYER)
251     ELSEIF (COUNT(XSOILGRID/=XUNDEF)/=NGROUND_LAYER) THEN
252       WRITE(ILUOUT,*) '********************************************************'
253       WRITE(ILUOUT,*) '* Soil grid reference values /= number of ground layer *'
254       WRITE(ILUOUT,*) '********************************************************'
255       CALL ABOR1_SFX('PGD_ISBA: XSOILGRID must be coherent with NGROUND_LAYER in NAM_ISBA') 
256     ELSEIF (XSOILGRID(1).GT.0.01) THEN
257       CALL ABOR1_SFX('PGD_ISBA: First layer of XSOILGRID must be lower than 1cm')
258     ENDIF
259 !
260     WRITE(ILUOUT,*) '*****************************************'
261     WRITE(ILUOUT,*) '* Option CISBA            = ',CISBA
262     WRITE(ILUOUT,*) '* Pedo transfert function = ',CPEDOTF    
263     WRITE(ILUOUT,*) '* Number of soil layers   = ',NGROUND_LAYER
264     IF(OECOCLIMAP)THEN
265       WRITE(ILUOUT,*) '* Soil layers grid (m)    = ',XSOILGRID(1:NGROUND_LAYER)
266     ENDIF
267     WRITE(ILUOUT,*) '*****************************************'
268 !    
269 END SELECT
270 !
271 SELECT CASE (CPHOTO)
272   CASE ('AGS','LAI','AST','LST')
273     NNBIOMASS = 1
274   CASE ('NIT')
275     NNBIOMASS = 3
276   CASE ('NCB')
277     NNBIOMASS = 6
278 END SELECT
279 WRITE(ILUOUT,*) '*****************************************'
280 WRITE(ILUOUT,*) '* With option CPHOTO = ',CPHOTO,'               *'
281 WRITE(ILUOUT,*) '* the number of biomass pools is set to ', NNBIOMASS
282 WRITE(ILUOUT,*) '*****************************************'
283 !
284 IF (NPATCH<1 .OR. NPATCH>NVEGTYPE) THEN
285   WRITE(ILUOUT,*) '*****************************************'
286   WRITE(ILUOUT,*) '* Number of patch must be between 1 and ', NVEGTYPE
287   WRITE(ILUOUT,*) '* You have chosen NPATCH = ', NPATCH
288   WRITE(ILUOUT,*) '*****************************************'
289   CALL ABOR1_SFX('PGD_ISBA: NPATCH MUST BE BETWEEN 1 AND NVEGTYPE')
290 END IF
291 !
292 IF ( CPHOTO/='NON' .AND. NPATCH/=12 ) THEN
293   WRITE(ILUOUT,*) '*****************************************'
294   WRITE(ILUOUT,*) '* With option CPHOTO = ', CPHOTO
295   WRITE(ILUOUT,*) '* Number of patch must be equal to 12 '
296   WRITE(ILUOUT,*) '* But you have chosen NPATCH = ', NPATCH
297   WRITE(ILUOUT,*) '*****************************************'
298   CALL ABOR1_SFX('PGD_ISBA: CPHOTO='//CPHOTO//' REQUIRES NPATCH=12')
299 END IF
300 !
301 IF ( CPHOTO=='NON' .AND. LTR_ML ) THEN
302   WRITE(ILUOUT,*) '*****************************************'
303   WRITE(ILUOUT,*) '* With option CPHOTO == NON '
304   WRITE(ILUOUT,*) '* New radiative transfert TR_ML  '
305   WRITE(ILUOUT,*) '* cant be used '
306   WRITE(ILUOUT,*) '*****************************************'
307   CALL ABOR1_SFX('PGD_ISBA: WITH CPHOTO= NON LTR_ML MUST BE FALSE')
308 END IF
309 !
310 !-------------------------------------------------------------------------------
311 !
312 !*    4.      Number of points and packing of general fields
313 !             ----------------------------------------------
314 !
315  CALL GET_SURF_SIZE_n('NATURE',ILU)
316 !
317 ALLOCATE(LCOVER     (JPCOVER))
318 ALLOCATE(XCOVER     (ILU,JPCOVER))
319 ALLOCATE(XZS        (ILU))
320 ALLOCATE(XLAT       (ILU))
321 ALLOCATE(XLON       (ILU))
322 ALLOCATE(XMESH_SIZE (ILU))
323 ALLOCATE(XZ0EFFJPDIR(ILU))
324 !
325  CALL PACK_PGD(HPROGRAM, 'NATURE',                    &
326                 CGRID,  XGRID_PAR,                     &
327                 LCOVER, XCOVER, XZS,                   &
328                 XLAT, XLON, XMESH_SIZE, XZ0EFFJPDIR    )  
329 !
330 !-------------------------------------------------------------------------------
331 !
332 !*    5.      Packing of ISBA specific fields
333 !             -------------------------------
334 !
335  CALL GET_AOS_n(HPROGRAM,NL,ZAOSIP,ZAOSIM,ZAOSJP,ZAOSJM,ZHO2IP,ZHO2IM,ZHO2JP,ZHO2JM)
336  CALL GET_SSO_n(HPROGRAM,NL,ZSSO_SLOPE)
337 !
338  CALL PACK_PGD_ISBA(HPROGRAM,                                    &
339                      ZAOSIP, ZAOSIM, ZAOSJP, ZAOSJM,              &
340                      ZHO2IP, ZHO2IM, ZHO2JP, ZHO2JM,              &
341                      ZSSO_SLOPE                                   )  
342 !
343 !-------------------------------------------------------------------------------
344 !
345 !*    6.      Topographic index for TOPMODEL
346 !             ------------------------------
347 !
348  CALL PGD_TOPO_INDEX(HPROGRAM,ILU,YCTI,YCTIFILETYPE,LIMP_CTI)
349 !
350 !-------------------------------------------------------------------------------
351 !
352 !*    7.      Sand fraction
353 !             -------------
354 !
355 CATYPE='ARI'
356 !
357 ALLOCATE(XSAND(ILU,NGROUND_LAYER))
358 !
359 IF(LIMP_SAND)THEN
360 !
361   IF(YSANDFILETYPE=='NETCDF')THEN
362      CALL ABOR1_SFX('Use another format than netcdf for sand input file with LIMP_SAND')
363   ELSE
364 #ifdef ASC
365      CFILEIN     = ADJUSTL(ADJUSTR(YSAND)//'.txt')
366 #endif
367 #ifdef FA
368      CFILEIN_FA  = ADJUSTL(ADJUSTR(YSAND)//'.fa')
369 #endif
370 #ifdef LFI
371      CFILEIN_LFI = ADJUSTL(YSAND)
372 #endif
373      CALL INIT_IO_SURF_n(YSANDFILETYPE,'NATURE','ISBA  ','READ ')
374   ENDIF     
375 !   
376   CALL READ_SURF(YSANDFILETYPE,'SAND',XSAND(:,1),IRESP) 
377 !
378   CALL END_IO_SURF_n(YSANDFILETYPE)
379 !
380 ELSE
381    CALL PGD_FIELD(HPROGRAM,'sand fraction','NAT',YSAND,YSANDFILETYPE,XUNIF_SAND,XSAND(:,1))
382 ENDIF
383 !
384 DO JLAYER=1,NGROUND_LAYER
385   XSAND(:,JLAYER) = XSAND(:,1)
386 END DO
387 !-------------------------------------------------------------------------------
388 !
389 !*    8.      Clay fraction
390 !             -------------
391 !
392 ALLOCATE(XCLAY(ILU,NGROUND_LAYER))
393 !
394 IF(LIMP_CLAY)THEN
395 !
396   IF(YCLAYFILETYPE=='NETCDF')THEN
397      CALL ABOR1_SFX('Use another format than netcdf for clay input file with LIMP_CLAY')
398   ELSE
399 #ifdef ASC
400      CFILEIN     = ADJUSTL(ADJUSTR(YSAND)//'.txt')
401 #endif
402 #ifdef FA
403      CFILEIN_FA  = ADJUSTL(ADJUSTR(YSAND)//'.fa')
404 #endif
405 #ifdef LFI
406      CFILEIN_LFI = ADJUSTL(YSAND)
407 #endif
408      CALL INIT_IO_SURF_n(YCLAYFILETYPE,'NATURE','ISBA  ','READ ')
409   ENDIF     
410 !   
411   CALL READ_SURF(YCLAYFILETYPE,'CLAY',XCLAY(:,1),IRESP) 
412 !
413   CALL END_IO_SURF_n(YCLAYFILETYPE)
414 !
415 ELSE
416   CALL PGD_FIELD(HPROGRAM,'clay fraction','NAT',YCLAY,YCLAYFILETYPE,XUNIF_CLAY,XCLAY(:,1))
417 ENDIF
418 !
419 DO JLAYER=1,NGROUND_LAYER
420   XCLAY(:,JLAYER) = XCLAY(:,1)
421 END DO
422 !
423 !-------------------------------------------------------------------------------
424 !
425 !*    9.      organic carbon profile
426 !             ----------------------
427 !
428 IF(LEN_TRIM(YSOCFILETYPE)/=0.OR.(XUNIF_SOC_TOP/=XUNDEF.AND.XUNIF_SOC_SUB/=XUNDEF))THEN
429 !
430   ALLOCATE(XSOC(ILU,NGROUND_LAYER))
431 !
432   LSOCP=.TRUE.
433 !
434   IF((LEN_TRIM(YSOC_TOP)==0.AND.LEN_TRIM(YSOC_SUB)/=0).OR.(LEN_TRIM(YSOC_TOP)/=0.AND.LEN_TRIM(YSOC_SUB)==0))THEN
435     WRITE(ILUOUT,*) ' '
436     WRITE(ILUOUT,*) '***********************************************************'
437     WRITE(ILUOUT,*) '* Error in soil organic carbon preparation                *'
438     WRITE(ILUOUT,*) '* If used, sub and top soil input file must be given      *'
439     WRITE(ILUOUT,*) '***********************************************************'
440     WRITE(ILUOUT,*) ' '
441     CALL ABOR1_SFX('PGD_ISBA: TOP AND SUB SOC INPUT FILE REQUIRED')        
442   ENDIF
443 !
444   IF(LIMP_SOC)THEN
445 !
446 !   Topsoil
447 !
448     IF(YSOCFILETYPE=='NETCDF')THEN
449        CALL ABOR1_SFX('Use another format than netcdf for organic carbon input file with LIMP_SOC')
450     ELSE
451 #ifdef ASC
452        CFILEIN     = ADJUSTL(ADJUSTR(YSOC_TOP)//'.txt')
453 #endif
454 #ifdef FA
455        CFILEIN_FA  = ADJUSTL(ADJUSTR(YSOC_TOP)//'.fa')
456 #endif
457 #ifdef LFI
458        CFILEIN_LFI = ADJUSTL(YSOC_TOP)
459 #endif
460        CALL INIT_IO_SURF_n(YSOCFILETYPE,'NATURE','ISBA  ','READ ')
461     ENDIF     
462 !   
463     CALL READ_SURF(YSOCFILETYPE,'SOC_TOP',XSOC(:,1),IRESP) 
464 !
465     CALL END_IO_SURF_n(YSOCFILETYPE)
466 !
467 !   Subsoil
468 !
469     IF(YSOCFILETYPE=='NETCDF')THEN
470        CALL ABOR1_SFX('Use another format than netcdf for organic carbon input file with LIMP_SOC')
471     ELSE
472 #ifdef ASC
473        CFILEIN     = ADJUSTL(ADJUSTR(YSOC_SUB)//'.txt')
474 #endif
475 #ifdef FA
476        CFILEIN_FA  = ADJUSTL(ADJUSTR(YSOC_SUB)//'.fa')
477 #endif
478 #ifdef LFI
479        CFILEIN_LFI = ADJUSTL(YSOC_SUB)
480 #endif
481        CALL INIT_IO_SURF_n(YSOCFILETYPE,'NATURE','ISBA  ','READ ')
482     ENDIF     
483 !   
484     CALL READ_SURF(YSOCFILETYPE,'SOC_SUB',XSOC(:,2),IRESP) 
485 !
486     CALL END_IO_SURF_n(YSOCFILETYPE)
487 !
488   ELSE
489     CALL PGD_FIELD(HPROGRAM,'organic carbon','NAT',YSOC_TOP,YSOCFILETYPE,XUNIF_SOC_TOP,XSOC(:,1))
490     CALL PGD_FIELD(HPROGRAM,'organic carbon','NAT',YSOC_SUB,YSOCFILETYPE,XUNIF_SOC_SUB,XSOC(:,2))
491   ENDIF
492 !
493   DO JLAYER=2,NGROUND_LAYER
494     XSOC(:,JLAYER) = XSOC(:,2)
495   END DO
496 !
497 ELSE
498 !
499   LSOCP=.FALSE.
500   ALLOCATE(XSOC(0,0))
501 !
502 ENDIF
503 !
504 !*    10.     Permafrost distribution
505 !             -----------------------
506 !
507 IF(LEN_TRIM(YPERM)/=0.OR.XUNIF_PERM/=XUNDEF)THEN
508 !
509   ALLOCATE(XPERM(ILU))
510 !
511   LPERM=.TRUE.
512 !
513   IF(LIMP_PERM)THEN
514 !
515     IF(YPERMFILETYPE=='NETCDF')THEN
516        CALL ABOR1_SFX('Use another format than netcdf for permafrost input file with LIMP_PERM')
517     ELSE
518 #ifdef ASC
519        CFILEIN     = ADJUSTL(ADJUSTR(YPERM)//'.txt')
520 #endif
521 #ifdef FA
522        CFILEIN_FA  = ADJUSTL(ADJUSTR(YPERM)//'.fa')
523 #endif
524 #ifdef LFI
525        CFILEIN_LFI = ADJUSTL(YPERM)
526 #endif
527        CALL INIT_IO_SURF_n(YPERMFILETYPE,'NATURE','ISBA  ','READ ')
528     ENDIF     
529 !   
530     CALL READ_SURF(YPERMFILETYPE,'PERM',XPERM(:),IRESP) 
531 !
532     CALL END_IO_SURF_n(YPERMFILETYPE)
533   ELSE
534     CALL PGD_FIELD(HPROGRAM,'permafrost','NAT',YPERM,YPERMFILETYPE,XUNIF_PERM,XPERM(:))
535   ENDIF
536 !
537 ELSE
538 !
539   LPERM=.FALSE.  
540   ALLOCATE(XPERM(0))
541 !
542 ENDIF
543 !
544 !-------------------------------------------------------------------------------
545 !
546 !*    11.  pH and fertlisation data
547 !             --------------------------
548 !
549 IF((LEN_TRIM(YPHFILETYPE)/=0.OR.XUNIF_PH/=XUNDEF) .AND. (LEN_TRIM(YFERTFILETYPE)/=0.OR.XUNIF_FERT/=XUNDEF)) THEN
550   !
551   ALLOCATE(XPH(ILU))
552   ALLOCATE(XFERT(ILU))
553   !
554   LNOF = .TRUE.
555   !
556   CALL PGD_FIELD(HPROGRAM,'pH value','NAT',YPH,YPHFILETYPE,XUNIF_PH,XPH(:))
557   CALL PGD_FIELD(HPROGRAM,'fertilisation','NAT',YFERT,YFERTFILETYPE,XUNIF_FERT,XFERT(:))
558   !
559 ENDIF
560 !
561 !-------------------------------------------------------------------------------
562 !
563 !*    12.      Subgrid runoff 
564 !             --------------
565 !
566 ALLOCATE(XRUNOFFB(ILU))
567  CALL PGD_FIELD                                                                              &
568        (HPROGRAM,'subgrid runoff','NAT',YRUNOFFB,YRUNOFFBFILETYPE,XUNIF_RUNOFFB,XRUNOFFB(:))  
569 !
570 !-------------------------------------------------------------------------------
571 !
572 !*    13.     Drainage coefficient
573 !             --------------------
574 !
575 ALLOCATE(XWDRAIN(ILU))
576  CALL PGD_FIELD                                                                              &
577        (HPROGRAM,'subgrid drainage','NAT',YWDRAIN,YWDRAINFILETYPE,XUNIF_WDRAIN,XWDRAIN(:))  
578 !
579 !-------------------------------------------------------------------------------
580 !
581 !*   14.      ISBA specific fields
582 !             --------------------
583 !
584 LECOCLIMAP = OECOCLIMAP
585 !
586  CALL PGD_ISBA_PAR(HPROGRAM)
587 !
588 !-------------------------------------------------------------------------------
589 !
590  CALL PGD_TOPD(HPROGRAM)
591 !
592 !-------------------------------------------------------------------------------
593 !
594 !*   15.     Prints of cover parameters in a tex file
595 !            ----------------------------------------
596 !
597 IF (OECOCLIMAP) THEN
598   CALL WRITE_COVER_TEX_ISBA    (NPATCH,NGROUND_LAYER,CISBA)
599   CALL WRITE_COVER_TEX_ISBA_PAR(NPATCH,NGROUND_LAYER,CISBA,CPHOTO,XSOILGRID)
600 END IF
601 IF (LHOOK) CALL DR_HOOK('PGD_ISBA',1,ZHOOK_HANDLE)
602 !
603 !-------------------------------------------------------------------------------
604 !
605 END SUBROUTINE PGD_ISBA