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