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