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