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