init_veg_pgdn.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!#############################################################
SUBROUTINE init_veg_pgd_n (ISSK, DTI, IO, S, K, KK, PK, PEK, AGK, KI, &
                           HPROGRAM, HSURF, KLUOUT, KSIZE, KMONTH,    &
                           ODEEPSOIL, OPHYSDOMC, PTDEEP_CLI, PGAMMAT_CLI,     &
                           OAGRIP, PTHRESHOLD, HINIT, PCO2, PRHOA     )  
!#############################################################
!
!!****  *INIT_VEG_PGD_n_n* - routine to initialize ISBA
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      23/07/13     (Decharme) Surface / Water table depth coupling
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_surfex_mpi, ONLY : nrank
!
USE modd_sso_n, ONLY : sso_t
!
USE modd_data_isba_n, ONLY : data_isba_t
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_k_t, isba_p_t, isba_pe_t
USE modd_agri_n, ONLY : agri_t
!
USE modd_surf_atm,       ONLY : lcpl_arp
USE modd_data_cover_par, ONLY : nvegtype
USE modd_surf_par,       ONLY : xundef, nundef
USE modd_csts,           ONLY : xcpd, xlvtt, xlstt
USE modd_snow_par,       ONLY : xemissn
USE modd_isba_par,       ONLY : xtau_ice
!
USE modd_sgh_par,        ONLY : xice_deph_max
!
USE mode_cotwo,          ONLY : gauleg
!
USE modi_surf_patch
USE modi_get_1d_mask
USE modi_co2_init_n
USE modi_subscale_z0eff
!
USE mode_soil
!
USE modi_heatcapz
USE modi_thrmcondz
USE modi_abor1_sfx
USE modi_dif_layer
USE modi_dry_wet_soil_albedos
USE modi_pack_same_rank
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              --²-----------------------
!
!
TYPE(sso_t), INTENT(INOUT) :: ISSK
TYPE(data_isba_t), INTENT(INOUT) :: DTI
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(isba_k_t), INTENT(INOUT) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
TYPE(agri_t), INTENT(INOUT) :: AGK
!
INTEGER, INTENT(IN) :: KI
!
 CHARACTER(LEN=6), INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=6), INTENT(IN)  :: HSURF     ! Type of surface
INTEGER, INTENT(IN)  :: KLUOUT
!
INTEGER, INTENT(IN)  :: KSIZE
!
INTEGER, INTENT(IN)  :: KMONTH
!
LOGICAL, INTENT(IN) :: ODEEPSOIL
LOGICAL, INTENT(IN) :: OPHYSDOMC
REAL, DIMENSION(:), INTENT(IN) :: PTDEEP_CLI
REAL, DIMENSION(:), INTENT(IN) :: PGAMMAT_CLI
!
LOGICAL, INTENT(IN) :: OAGRIP
REAL, DIMENSION(:), INTENT(IN) :: PTHRESHOLD
!
 CHARACTER(LEN=3), INTENT(IN) :: HINIT
 !
REAL, DIMENSION(:), INTENT(IN) :: PCO2
REAL, DIMENSION(:), INTENT(IN) :: PRHOA
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
INTEGER :: JPATCH  ! loop counter on tiles
INTEGER :: JILU,JP, JMAXLOC    ! loop increment
INTEGER :: JL  ! loop counter on layers
!
INTEGER :: IABC
!
REAL, DIMENSION(SIZE(PCO2))       :: ZCO2  ! CO2 concentration  (kg/kg)
!
INTEGER, DIMENSION(:), ALLOCATABLE :: IR_NATURE_P
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!               Initialisation for IO
!
IF (lhook) CALL dr_hook('INIT_VEG_PGD_n',0,zhook_handle)
!
!------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------
!
!        PART 1: fields that are needed unpacked and packed: defined unpacked
!        -------------------------------------------------------------------
!
!*          Soil hydraulic characteristics:
!           -------------------------------
!
IF (.NOT.ASSOCIATED(k%XMPOTSAT)) THEN
  !
  ALLOCATE(k%XMPOTSAT (ki,io%NGROUND_LAYER))
  ALLOCATE(k%XBCOEF   (ki,io%NGROUND_LAYER))
  ALLOCATE(k%XWWILT   (ki,io%NGROUND_LAYER)) ! wilting point
  ALLOCATE(k%XWFC     (ki,io%NGROUND_LAYER)) ! field capacity
  ALLOCATE(k%XWSAT    (ki,io%NGROUND_LAYER)) ! saturation
  !
  DO jl=1,io%NGROUND_LAYER
    IF (dti%LDATA_BCOEF) THEN
      k%XBCOEF  (:,jl) = dti%XPAR_BCOEF  (:,jl)
    ELSE
      k%XBCOEF  (:,jl) = bcoef_func(k%XCLAY(:,jl),k%XSAND(:,jl),io%CPEDOTF)
    ENDIF
    IF (dti%LDATA_MPOTSAT) THEN
      k%XMPOTSAT(:,jl) = dti%XPAR_MPOTSAT(:,jl)
    ELSE
      k%XMPOTSAT(:,jl) = matpotsat_func(k%XCLAY(:,jl),k%XSAND(:,jl),io%CPEDOTF)
    ENDIF
    IF (dti%LDATA_WSAT) THEN
      k%XWSAT   (:,jl) = dti%XPAR_WSAT   (:,jl)
    ELSE
      k%XWSAT   (:,jl) = wsat_func(k%XCLAY(:,jl),k%XSAND(:,jl),io%CPEDOTF)
    ENDIF
    IF (dti%LDATA_WWILT) THEN
      k%XWWILT   (:,jl) = dti%XPAR_WWILT  (:,jl)
    ELSE
      k%XWWILT  (:,jl) = wwilt_func(k%XCLAY(:,jl),k%XSAND(:,jl),io%CPEDOTF)
    ENDIF
  END DO
  IF (dti%LDATA_BCOEF  ) DEALLOCATE(dti%XPAR_BCOEF)
  IF (dti%LDATA_MPOTSAT) DEALLOCATE(dti%XPAR_MPOTSAT)
  IF (dti%LDATA_WSAT   ) DEALLOCATE(dti%XPAR_WSAT)
  IF (dti%LDATA_WWILT  ) DEALLOCATE(dti%XPAR_WWILT)
  !
  IF (dti%LDATA_WFC) THEN
    k%XWFC(:,:) = dti%XPAR_WFC(:,:)
    DEALLOCATE(dti%XPAR_WFC)
  ELSEIF (io%CISBA=='2-L' .OR. io%CISBA=='3-L') THEN
    !  field capacity at hydraulic conductivity = 0.1mm/day
    k%XWFC(:,:) = wfc_func(k%XCLAY(:,:),k%XSAND(:,:),io%CPEDOTF)
  ELSE IF (io%CISBA=='DIF') THEN
    !  field capacity at water potential = 0.33bar        
    k%XWFC(:,:) = w33_func(k%XCLAY(:,:),k%XSAND(:,:),io%CPEDOTF)
  END IF
 !
  IF (io%CISBA=='2-L' .OR. io%CISBA=='3-L') THEN
    ALLOCATE(k%XCGSAT (ki))
    ALLOCATE(k%XC4B   (ki))
    ALLOCATE(k%XACOEF (ki))
    ALLOCATE(k%XPCOEF (ki))
    k%XCGSAT(:)  = cgsat_func(k%XCLAY(:,1),k%XSAND(:,1))
    k%XC4B  (:)  = c4b_func(k%XCLAY(:,1))
    k%XACOEF(:)  = acoef_func(k%XCLAY(:,1))
    k%XPCOEF(:)  = pcoef_func(k%XCLAY(:,1))
  ELSE IF (io%CISBA=='DIF') THEN
    ALLOCATE(k%XCGSAT (0))
    ALLOCATE(k%XC4B   (0))
    ALLOCATE(k%XACOEF (0))
    ALLOCATE(k%XPCOEF (0))
  ENDIF
  !
  IF(io%CRUNOFF=='SGH')THEN
    !
    ALLOCATE(k%XWD0   (ki,io%NGROUND_LAYER))
    ALLOCATE(k%XKANISO(ki,io%NGROUND_LAYER))
    !
    IF(io%CISBA=='DIF')THEN
      k%XWD0(:,:) = wfc_func(k%XCLAY(:,:),k%XSAND(:,:),io%CPEDOTF)
    ELSE
      k%XWD0(:,:) = k%XWWILT(:,:)
    ENDIF
    k%XKANISO(:,:) = aniso_func(k%XCLAY(:,:))
    !
  ELSE
    !
    ALLOCATE(k%XWD0   (0,0))
    ALLOCATE(k%XKANISO(0,0))
    !
  ENDIF
  !
  IF (io%CSCOND=='PL98'.OR.io%CISBA=='DIF') THEN
    ALLOCATE(k%XHCAPSOIL(ki,io%NGROUND_LAYER))
    ALLOCATE(k%XCONDDRY (ki,io%NGROUND_LAYER))
    ALLOCATE(k%XCONDSLD (ki,io%NGROUND_LAYER))
    ! 
    CALL heatcapz(k%XSAND,k%XHCAPSOIL)
    CALL thrmcondz(k%XSAND,k%XWSAT,k%XCONDDRY,k%XCONDSLD)
  ELSE
    ALLOCATE(k%XHCAPSOIL(0,0))
    ALLOCATE(k%XCONDDRY (0,0))
    ALLOCATE(k%XCONDSLD (0,0))
  END IF
  !
ENDIF
!
!CSCOND used in soil.F90 and soildif.F90
!
IF (io%CSCOND=='NP89'.AND.io%CISBA=='DIF') THEN
   WRITE(kluout,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
   WRITE(kluout,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
   WRITE(kluout,*)'IF CISBA=DIF, CSCOND=NP89 is not available'
   WRITE(kluout,*)'because not physic. CSCOND is put to PL98 '
   WRITE(kluout,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
   WRITE(kluout,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
ENDIF
!
!CPSURF used in drag.F90
!CPL_ARP used in drag.F90 and e_budget.F90
IF(io%CCPSURF=='DRY'.AND.lcpl_arp) THEN
  CALL abor1_sfx('CCPSURF=DRY must not be used with LCPL_ARP')
ENDIF
!
!------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------
!
!        PART 2: fields that are needed only packed: defined packed directly
!        -------------------------------------------------------------------
!
!        PART 2: A: fields that don't depend on patches: KK, AGK, ISSK
!        -------------------------------------------------------------
!
!*       2.A.1.    Miscellaneous fields for ISBA:
!        ----------------------------------------
!
!* default value for:
! lateral water flux, deep soil temperature climatology and its relaxation time-scale
!
!these arrays are used only packed: we define them directly packed
ALLOCATE(kk%XTDEEP (ksize))
ALLOCATE(kk%XGAMMAT(ksize))
kk%XTDEEP (:) = xundef
kk%XGAMMAT(:) = xundef
!
IF (odeepsoil) THEN
  DO jilu = 1, ksize
    kk%XTDEEP (jilu) = ptdeep_cli(kmonth)
    kk%XGAMMAT(jilu) = 1. / pgammat_cli(kmonth)
  END DO
  !
  WRITE(kluout,*)' LDEEPSOIL = ',odeepsoil,' LPHYSDOMC = ',ophysdomc
  WRITE(kluout,*)' XTDEEP    = ',minval(kk%XTDEEP(:)) ,maxval(kk%XTDEEP(:))
  WRITE(kluout,*)' XGAMMAT   = ',minval(kk%XGAMMAT(:)),maxval(kk%XGAMMAT(:))
ENDIF
!
!
!*         2.A.2. Initialize hydrology
!          ---------------------------
!
IF (io%CISBA == 'DIF') THEN
  !
  ALLOCATE(kk%XFWTD(ksize))
  ALLOCATE(kk%XWTD (ksize))
  kk%XFWTD(:) = 0.0
  kk%XWTD (:) = xundef
  !
ELSE
  !
  ALLOCATE(kk%XFWTD(0))
  ALLOCATE(kk%XWTD (0))
  !   
ENDIF
!
!
!*         Physiographic Radiative fields:  
!               ------------------------------
!
!
!*        2.A.3. dry and wet bare soil albedos
!         ------------------------------------
!
ALLOCATE(kk%XALBNIR_DRY  (ksize))
ALLOCATE(kk%XALBVIS_DRY  (ksize))
ALLOCATE(kk%XALBUV_DRY   (ksize))
ALLOCATE(kk%XALBNIR_WET  (ksize))
ALLOCATE(kk%XALBVIS_WET  (ksize))
ALLOCATE(kk%XALBUV_WET   (ksize))
!
 CALL dry_wet_soil_albedos(kk )
!

!
!*       2.A.4. Irrigation
!        -----------------
!
IF (oagrip) THEN
  !
  ALLOCATE(agk%NIRRINUM     (ksize))
  ALLOCATE(agk%LIRRIDAY     (ksize))
  ALLOCATE(agk%LIRRIGATE    (ksize))
  ALLOCATE(agk%XTHRESHOLDSPT(ksize))
  !
  agk%NIRRINUM (:) = 1
  agk%LIRRIDAY (:) = .false.                          
  agk%LIRRIGATE(:) = .false.                          
  !
  DO jilu = 1, ksize
    agk%XTHRESHOLDSPT(jilu) = pthreshold(agk%NIRRINUM(jilu))
  END DO
ELSE
  ALLOCATE(agk%NIRRINUM     (0))
  ALLOCATE(agk%LIRRIDAY     (0))
  ALLOCATE(agk%LIRRIGATE    (0))
  ALLOCATE(agk%XTHRESHOLDSPT(0))
ENDIF
!
!*       2.A.5. Orographic roughness length
!        ----------------------------------
!
ALLOCATE(issk%XZ0EFFIP(ksize))
ALLOCATE(issk%XZ0EFFIM(ksize))
ALLOCATE(issk%XZ0EFFJP(ksize))
ALLOCATE(issk%XZ0EFFJM(ksize))
!
issk%XZ0EFFIP(:) = xundef
issk%XZ0EFFIM(:) = xundef
issk%XZ0EFFJP(:) = xundef
issk%XZ0EFFJM(:) = xundef
!
IF (SIZE(issk%XAOSIP)>0) CALL subscale_z0eff(issk,pek%XZ0,.false.)
!
!-----------------------------------------------------------------------
!
!        PART 2: B: fields that depend on patches: PK, PEK
!        -------------------------------------------------
!
!
!*       2.B.1. Additional fields for ISBA-AGS:
!        --------------------------------------                 
!
IF(io%CPHOTO /= 'NON' .AND. hinit == 'ALL') THEN
  !
  IF (.NOT.ASSOCIATED(s%XABC)) THEN
    IF (io%LTR_ML) THEN
      iabc = 10
    ELSE
      iabc = 3
    ENDIF
    ALLOCATE(s%XABC(iabc))
    ALLOCATE(s%XPOI(iabc))
    s%XABC(:) = 0.
    s%XPOI(:) = 0.       
    CALL gauleg(0.0,1.0,s%XABC,s%XPOI,iabc)
  ENDIF
  !
  zco2(:) = pco2(:) / prhoa(:)
  ALLOCATE(pk%XANMAX        (ksize))
  ALLOCATE(pk%XFZERO        (ksize))
  ALLOCATE(pk%XEPSO         (ksize))
  ALLOCATE(pk%XGAMM         (ksize))
  ALLOCATE(pk%XQDGAMM       (ksize))
  ALLOCATE(pk%XQDGMES       (ksize))
  ALLOCATE(pk%XT1GMES       (ksize))
  ALLOCATE(pk%XT2GMES       (ksize))
  ALLOCATE(pk%XAMAX         (ksize))
  ALLOCATE(pk%XQDAMAX       (ksize))
  ALLOCATE(pk%XT1AMAX       (ksize))
  ALLOCATE(pk%XT2AMAX       (ksize))
  ALLOCATE(pk%XAH           (ksize))
  ALLOCATE(pk%XBH           (ksize))
  ALLOCATE(pk%XTAU_WOOD     (ksize))
  ALLOCATE(pk%XINCREASE     (ksize,io%NNBIOMASS))
  ALLOCATE(pk%XTURNOVER     (ksize,io%NNBIOMASS))
  CALL co2_init_n(io, s, pk, pek, ksize, zco2  )
  !
ELSEIF(io%CPHOTO == 'NON' .AND. io%LTR_ML) THEN ! Case for MEB
   !
   IF (.NOT.ASSOCIATED(s%XABC)) THEN
     iabc = 10
     ALLOCATE (s%XABC(iabc))
     ALLOCATE (s%XPOI(iabc)) ! Working
     s%XABC(:) = 0.
     s%XPOI(:) = 0.
     CALL gauleg(0.0,1.0,s%XABC,s%XPOI,iabc)
     DEALLOCATE (s%XPOI)
     ALLOCATE   (s%XPOI(0))
   ENDIF
   !
ELSE
  !
  IF (.NOT.ASSOCIATED(s%XABC)) THEN
    ALLOCATE(s%XABC(0))
    ALLOCATE(s%XPOI(0))
  ENDIF
  !
  ALLOCATE(pk%XANMAX        (0))
  ALLOCATE(pk%XFZERO        (0))
  ALLOCATE(pk%XEPSO         (0))
  ALLOCATE(pk%XGAMM         (0))
  ALLOCATE(pk%XQDGAMM       (0))
  ALLOCATE(pk%XQDGMES       (0))
  ALLOCATE(pk%XT1GMES       (0))
  ALLOCATE(pk%XT2GMES       (0))
  ALLOCATE(pk%XAMAX         (0))
  ALLOCATE(pk%XQDAMAX       (0))
  ALLOCATE(pk%XT1AMAX       (0))
  ALLOCATE(pk%XT2AMAX       (0))
  ALLOCATE(pk%XAH           (0))
  ALLOCATE(pk%XBH           (0))
  ALLOCATE(pk%XTAU_WOOD     (0))
  ALLOCATE(pk%XINCREASE     (0,0))
  ALLOCATE(pk%XTURNOVER     (0,0))
  !
END IF
!
!
!*          2.B.2. Soil hydraulic characteristics (rest) :
!           --------------------------------------------
!
!
ALLOCATE(pk%XCONDSAT (ksize,io%NGROUND_LAYER))
ALLOCATE(pk%XTAUICE  (ksize))
!
IF (dti%LDATA_CONDSAT) THEN
  CALL pack_same_rank(pk%NR_P,dti%XPAR_CONDSAT(:,:),pk%XCONDSAT(:,:))
ELSE
  DO jl=1,io%NGROUND_LAYER
    pk%XCONDSAT(:,jl) = hydcondsat_func(kk%XCLAY(:,jl),kk%XSAND(:,jl),io%CPEDOTF) 
  END DO
ENDIF
pk%XTAUICE(:) = xtau_ice
!
IF (io%CISBA=='2-L' .OR. io%CISBA=='3-L') THEN
  !
  ALLOCATE(pk%XC1SAT (ksize))
  ALLOCATE(pk%XC2REF (ksize))
  ALLOCATE(pk%XC3    (ksize,2))  
  ALLOCATE(pk%XC4REF (ksize))      
  pk%XC1SAT(:) = c1sat_func(kk%XCLAY(:,1))
  pk%XC2REF(:) = c2ref_func(kk%XCLAY(:,1))
  pk%XC3 (:,1) = c3_func(kk%XCLAY(:,1))
  pk%XC3 (:,2) = c3_func(kk%XCLAY(:,2))
  !
  pk%XC4REF(:) = c4ref_func(kk%XCLAY(:,1),kk%XSAND(:,1),pk%XDG(:,2), &
                                        pk%XDG(:,io%NGROUND_LAYER)  )
  !
ELSE IF (io%CISBA=='DIF') THEN
  !
  ALLOCATE(pk%XC1SAT (0))
  ALLOCATE(pk%XC2REF (0))
  ALLOCATE(pk%XC3    (0,0))  
  ALLOCATE(pk%XC4REF (0))  
  !
END IF
!
ALLOCATE(pk%XCPS (ksize))
ALLOCATE(pk%XLVTT(ksize))
ALLOCATE(pk%XLSTT(ksize))
pk%XCPS (:) = xcpd
pk%XLVTT(:) = xlvtt
pk%XLSTT(:) = xlstt
!
!
!*       2.B.3.  Initialize hydrology
!        ----------------------------
!
ALLOCATE(pk%XRUNOFFD (ksize))
pk%XRUNOFFD(:)=xundef
!
IF (io%CISBA == 'DIF') THEN
!  
  ALLOCATE(pk%XDZG       (ksize,io%NGROUND_LAYER))
  ALLOCATE(pk%XDZDIF     (ksize,io%NGROUND_LAYER))
  ALLOCATE(pk%XSOILWGHT  (ksize,io%NGROUND_LAYER))
  CALL dif_layer(ksize, io, pk )
!
ELSEIF (count(io%LMEB_PATCH(:))/=0) THEN
!
  ALLOCATE(pk%XDZG       (ksize,io%NGROUND_LAYER))
  CALL dif_layer(ksize, io, pk, omeb_3l = .true.)
!
ELSE
!   
  ALLOCATE(pk%XDZG       (0,0))
  ALLOCATE(pk%XDZDIF     (0,0))
  ALLOCATE(pk%XSOILWGHT  (0,0))
  !
  WHERE(pk%XPATCH(:)>0.0)
    pk%XRUNOFFD(:) = pk%XDG(:,2)
  ENDWHERE
!   
ENDIF
!
!Horton (also used by the flooding sheme)
! 
ALLOCATE(pk%XKSAT_ICE(ksize))
!
IF(io%CISBA/='DIF')THEN
  pk%XD_ICE   (:) = min(pk%XDG(:,2),pk%XD_ICE(:))
  pk%XD_ICE   (:) = max(xice_deph_max,pk%XD_ICE(:))
  pk%XKSAT_ICE(:) = pk%XCONDSAT(:,1)
ELSE
  pk%XD_ICE   (:) = 0.0
  pk%XKSAT_ICE(:) = 0.0
ENDIF
!
!-------------------------------------------------------------------------------
!
!*        Physiographic Radiative fields:  
!         ------------------------------
!

!
!*       2.B.4.  Nitrogen version for isbaAgs
!        ------------------------------------                        
!
IF (io%CPHOTO=='NIT' .OR. io%CPHOTO=='NCB') THEN
  ALLOCATE(pk%XBSLAI_NITRO (ksize ))
  WHERE ((pek%XCE_NITRO(:) * pek%XCNA_NITRO(:) + pek%XCF_NITRO (:)) /= 0. )
      pk%XBSLAI_NITRO(:) = 1. / (pek%XCE_NITRO (:)*pek%XCNA_NITRO(:)+pek%XCF_NITRO (:))
  ELSEWHERE
      pk%XBSLAI_NITRO(:) = xundef
  ENDWHERE
ELSE
  ALLOCATE(pk%XBSLAI_NITRO (0))
ENDIF
!
IF (lhook) CALL dr_hook('INIT_VEG_PGD_n',1,zhook_handle)
!
END SUBROUTINE init_veg_pgd_n