read_isban.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 read_isba_n (DTCO, IO, S, NP, NPE, PCLAY, U, HPROGRAM)
!     ##################################
!
!!****  *READ_ISBA_n* - routine to initialise ISBA variables
!!                         
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2003
!!
!!      READ_SURF for general reading : 08/2003 (S.Malardel)
!!      B. Decharme  2008    : Floodplains
!!      B. Decharme  01/2009 : Optional Arpege deep soil temperature read
!!      A.L. Gibelin   03/09 : modifications for CENTURY model 
!!      A.L. Gibelin    04/2009 : BIOMASS and RESP_BIOMASS arrays 
!!      A.L. Gibelin    06/2009 : Soil carbon variables for CNT option
!!      B. Decharme  09/2012 : suppress NWG_LAYER (parallelization problems)
!!      T. Aspelien  08/2013 : Read diagnostics for assimilation
!!      P. Samuelsson   10/2014 : MEB
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_np_t, isba_npe_t, isba_p_t, isba_pe_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_co2v_par,       ONLY : xanfminit, xcondctmin
!                          
USE modd_assim,          ONLY : lassim,cassim_isba,xat2m_isba,xahu2m_isba,&
                                xazon10m_isba,xamer10m_isba,nific,nvar, &
                                cobs,nobstype,cvar,lprt,xtprt,nivar,cbio, &
                                xaddinfl,nens,xsigma,nie
!                                
USE modd_surf_par,       ONLY : xundef, nundef
USE modd_snow_par,       ONLY : xz0sn
USE modd_isba_par,       ONLY : xwgmin
!
USE mode_read_surf_layers
!
USE modi_read_surf
USE modi_make_choice_array
USE modi_pack_same_rank
!
USE modi_read_gr_snow
USE modi_abor1_sfx
USE modi_io_buff
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_io_buff_clean
USE modi_get_type_dim_n
USE mode_random
USE mode_ekf
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENT(INOUT) :: NPE
REAL, DIMENSION(:,:), INTENT(IN) :: PCLAY
TYPE(surf_atm_t), INTENT(INOUT) :: U
!
 CHARACTER(LEN=6),  INTENT(IN)  :: HPROGRAM ! calling program
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
TYPE(isba_p_t), POINTER :: PK
TYPE(isba_pe_t), POINTER :: PEK
!
INTEGER           :: ILU          ! 1D physical dimension
!
INTEGER           :: IRESP          ! Error code after redding
!
 CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
 CHARACTER(LEN=12) :: YCBIO          ! Name of biomass variable
!
 CHARACTER(LEN=4)  :: YLVL
!
REAL, DIMENSION(:,:,:),ALLOCATABLE :: ZLAI
REAL, DIMENSION(:,:,:),ALLOCATABLE  :: ZWORK3D    ! 2D array to write data in file
REAL, DIMENSION(:,:),ALLOCATABLE  :: ZWORK      ! 2D array to write data in file
REAL, DIMENSION(:), ALLOCATABLE :: ZCOFSWI
!
REAL,DIMENSION(IO%NPATCH) :: ZVLAIMIN
REAL :: ZCOEF
!
INTEGER :: IWORK   ! Work integer
!
INTEGER :: JP, JL, JNB, JNLITTER, JNSOILCARB, JNLITTLEVS  ! loop counter on layers
INTEGER :: JVAR, JI
!
INTEGER           :: IVERSION       ! surface version
INTEGER           :: IBUGFIX
INTEGER           :: IIVAR
INTEGER           :: IOBS
INTEGER           :: IBSUP
INTEGER           :: ISIZE_LMEB_PATCH, IMASK
!
LOGICAL :: GKNOWN, GDIM
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!
!* 1D physical dimension
!
IF (lhook) CALL dr_hook('READ_ISBA_N',0,zhook_handle)
yrecfm='SIZE_NATURE'
 CALL get_type_dim_n(dtco, u, 'NATURE',ilu)
!
yrecfm='VERSION'
 CALL read_surf(hprogram,yrecfm,iversion,iresp)
!
yrecfm='BUG'
 CALL read_surf(hprogram,yrecfm,ibugfix,iresp)
!
gdim = (iversion>8 .OR. iversion==8 .AND. ibugfix>0)
IF (gdim) CALL read_surf(hprogram,'SPLIT_PATCH',gdim,iresp)
!
yrecfm='BUG'
 CALL read_surf(hprogram,yrecfm,ibugfix,iresp)

!*       2.     Prognostic fields:
!               -----------------
!
ALLOCATE(zwork(ilu,io%NPATCH))
!* soil temperatures
!
IF(io%LTEMP_ARP)THEN
  iwork=io%NTEMPLAYER_ARP
ELSEIF(io%CISBA=='DIF')THEN
  iwork=io%NGROUND_LAYER
ELSE
  iwork=2 !Only 2 temperature layer in ISBA-FR
ENDIF
!
IF ( trim(cassim_isba)=="ENKF") THEN
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    ALLOCATE(pk%XRED_NOISE(pk%NSIZE_P,nvar))
    pk%XRED_NOISE(:,:) = 0.
  ENDDO
ELSE
  DO jp = 1,io%NPATCH
    ALLOCATE(np%AL(jp)%XRED_NOISE(0,0))
  ENDDO  
ENDIF
!
IF  ( trim(cassim_isba)=="ENKF" .OR. (trim(cassim_isba)=="EKF" .AND. lprt) ) THEN
  ALLOCATE(zcofswi(ilu))
  CALL cofswi(pclay(:,1),zcofswi)
ELSE
  ALLOCATE(zcofswi(0))
ENDIF
!
DO jp = 1,io%NPATCH
  pek => npe%AL(jp)
  pk => np%AL(jp)
  !
  ALLOCATE(pek%XTG(pk%NSIZE_P,iwork))
  !
  ALLOCATE(pek%XWG (pk%NSIZE_P,io%NGROUND_LAYER))
  ALLOCATE(pek%XWGI(pk%NSIZE_P,io%NGROUND_LAYER))
  !
  pek%XTG (:,:) = xundef
  pek%XWG (:,:) = xundef
  pek%XWGI(:,:) = xundef
  !
  ALLOCATE(pek%XWR(pk%NSIZE_P))
  !
ENDDO
!
ALLOCATE(zwork3d(ilu,iwork,io%NPATCH))
 CALL read_surf_layers(hprogram,'TG',gdim,zwork3d,iresp)
DO jl=1,iwork
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork3d(:,jl,jp),npe%AL(jp)%XTG(:,jl))
  ENDDO
ENDDO
DEALLOCATE(zwork3d)
!
! Perturb value if requested
IF ( trim(cassim_isba)=="EKF" .AND. lprt ) THEN
  !
  DO jl=1,iwork
  ! read in control variable
    IF ( (trim(cvar(nivar))=="TG1" .AND. jl==1) .OR. (trim(cvar(nivar))=="TG2" .AND. jl==2) ) THEN
      DO jp = 1,io%NPATCH
        pek => npe%AL(jp)
        WHERE ( pek%XTG(:,jl)/=xundef )
          pek%XTG(:,jl) = pek%XTG(:,jl) + xtprt(nivar)*pek%XTG(:,jl)
        ENDWHERE
      ENDDO
    ENDIF
  END DO
  !
ELSEIF ( trim(cassim_isba)=="ENKF" .AND. nie<nens+1 ) THEN
  !
  CALL make_ens_enkf(iwork,ilu,"TG ",zcofswi,np)
  !
ENDIF
!
!
!* soil liquid and ice water contents
!
ALLOCATE(zwork3d(ilu,io%NGROUND_LAYER,io%NPATCH))
 CALL read_surf_layers(hprogram,'WG',gdim,zwork3d,iresp)
DO jl=1,io%NGROUND_LAYER
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork3d(:,jl,jp),npe%AL(jp)%XWG(:,jl))
  ENDDO
ENDDO
DEALLOCATE(zwork3d)
!
! Perturb value if requested
IF ( trim(cassim_isba)=="EKF" .AND. lprt ) THEN
   !
   DO jl=1,io%NGROUND_LAYER
     ! read in control variable
     IF ( (trim(cvar(nivar))=="WG1" .AND. jl==1) .OR. & 
          (trim(cvar(nivar))=="WG2" .AND. jl==2) .OR. &
          (trim(cvar(nivar))=="WG3" .AND. jl==3) .OR. &
          (trim(cvar(nivar))=="WG4" .AND. jl==4) .OR. &
          (trim(cvar(nivar))=="WG5" .AND. jl==5) .OR. &
          (trim(cvar(nivar))=="WG6" .AND. jl==6) .OR. &
          (trim(cvar(nivar))=="WG7" .AND. jl==7) .OR. &
          (trim(cvar(nivar))=="WG8" .AND. jl==8) ) THEN
            
       DO jp = 1,io%NPATCH
         pek => npe%AL(jp)
         pk => np%AL(jp)
         DO ji = 1,pk%NSIZE_P
           imask = pk%NR_P(ji)
           IF (pek%XWG(ji,jl)/=xundef ) THEN
             pek%XWG(ji,jl) = pek%XWG(ji,jl) + xtprt(nivar) * zcofswi(imask) 
           ENDIF
         ENDDO
       END DO
       !
     ENDIF
     !
   END DO
   !
ELSEIF ( trim(cassim_isba)=="ENKF" .AND. nie<nens+1 ) THEN
  !
  CALL make_ens_enkf(iwork,ilu,"WG ",zcofswi,np)
  !
ENDIF
!
IF(io%CISBA=='DIF')THEN
  iwork=io%NGROUND_LAYER
ELSE
  iwork=2 !Only 2 soil ice layer in ISBA-FR
ENDIF
!
ALLOCATE(zwork3d(ilu,iwork,io%NPATCH))
 CALL read_surf_layers(hprogram,'WGI',gdim,zwork3d,iresp)
DO jl=1,iwork
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork3d(:,jl,jp),npe%AL(jp)%XWGI(:,jl))
  ENDDO
ENDDO
DEALLOCATE(zwork3d)
!
!* water intercepted on leaves
!
!
yrecfm = 'WR'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XWR(:))
ENDDO
!
!* Leaf Area Index
!
IF (io%CPHOTO=='NIT' .OR. io%CPHOTO=='NCB') THEN

  yrecfm = 'LAI'
  CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)  
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XLAI(:))
  ENDDO

  IF ( trim(cassim_isba)=="EKF" .AND. lprt ) THEN
    !
    ! read in control variable
    IF ( trim(cvar(nivar))=="LAI" ) THEN
      DO jp = 1,io%NPATCH
        pek => npe%AL(jp)
        WHERE ( pek%XLAI(:)/=xundef ) 
           pek%XLAI(:) =  pek%XLAI(:) + xtprt(nivar)* pek%XLAI(:)
        ENDWHERE
      ENDDO
    ENDIF
    !
  ELSEIF ( trim(cassim_isba)=="ENKF" .AND. nie<nens+1 ) THEN
    !
    IF (io%NPATCH==12) THEN
      zvlaimin = (/0.3,0.3,0.3,0.3,1.0,1.0,0.3,0.3,0.3,0.3,0.3,0.3/)
    ELSE
      zvlaimin = (/0.3/)
    ENDIF
    !
    ALLOCATE(zlai(ilu,1,io%NPATCH))
    zlai(:,:,:) = 0.
    DO jp = 1,io%NPATCH
      zlai(1:np%AL(jp)%NSIZE_P,1,jp) = npe%AL(jp)%XLAI(:)
    ENDDO
    CALL make_ens_enkf(1,ilu,"LAI",zcofswi,np,zlai)
    DO jp = 1,io%NPATCH
      npe%AL(jp)%XLAI(:) = max(zvlaimin(jp),zlai(1:np%AL(jp)%NSIZE_P,1,jp))
    ENDDO
    DEALLOCATE(zlai)
    !    
  ENDIF  
END IF
!
!* snow mantel
!
DO jp = 1,io%NPATCH
  IF (jp>1) THEN
    npe%AL(jp)%TSNOW%SCHEME = npe%AL(1)%TSNOW%SCHEME
    npe%AL(jp)%TSNOW%NLAYER = npe%AL(1)%TSNOW%NLAYER
  ENDIF
  CALL read_gr_snow(hprogram,'VEG','     ',ilu,np%AL(jp)%NSIZE_P,np%AL(jp)%NR_P,jp,&
                    npe%AL(jp)%TSNOW,knpatch=io%NPATCH  )
  CALL io_buff_clean
ENDDO
!
IF(io%LGLACIER)THEN
  DO jp = 1,io%NPATCH
    ALLOCATE(npe%AL(jp)%XICE_STO(np%AL(jp)%NSIZE_P))
  ENDDO
  IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=2) THEN
    yrecfm = 'ICE_STO'
    CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
    DO jp = 1,io%NPATCH
      CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XICE_STO(:))
    ENDDO
  ELSE
    DO jp = 1,io%NPATCH
      npe%AL(jp)%XICE_STO(:) = 0.0
    ENDDO
  ENDIF
ELSE  
  DO jp = 1,io%NPATCH
    ALLOCATE(npe%AL(jp)%XICE_STO(0))
  ENDDO
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       3.  MEB Prognostic or Semi-prognostic variables
!            -------------------------------------------
!
isize_lmeb_patch=count(io%LMEB_PATCH(:))
!
IF (isize_lmeb_patch>0) THEN
!
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    pek => npe%AL(jp)
    ALLOCATE(pek%XWRL (pk%NSIZE_P))
    ALLOCATE(pek%XWRLI(pk%NSIZE_P))
    ALLOCATE(pek%XWRVN(pk%NSIZE_P))
    ALLOCATE(pek%XTV  (pk%NSIZE_P))
    ALLOCATE(pek%XTL  (pk%NSIZE_P))
    ALLOCATE(pek%XTC  (pk%NSIZE_P))
    ALLOCATE(pek%XQC  (pk%NSIZE_P))
  ENDDO

!* water intercepted on litter
 
 yrecfm = 'WRL'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XWRL(:))
ENDDO

 yrecfm = 'WRLI'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XWRLI(:))
ENDDO 
!
!* snow intercepted on vegetation canopy leaves
!
  yrecfm = 'WRVN'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XWRVN(:))
ENDDO   
!
!* vegetation canopy temperature
!
  yrecfm = 'TV'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XTV(:))
ENDDO    
!
!* litter temperature
!
  yrecfm = 'TL'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XTL(:))
ENDDO  
!
!* vegetation canopy air temperature
!
  yrecfm = 'TC'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XTC(:))
ENDDO    
!
!* vegetation canopy air specific humidity
!
  yrecfm = 'QC'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XQC(:))
ENDDO     
!
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       4.  Semi-prognostic variables
!            -------------------------
!
DO jp = 1,io%NPATCH
  pk => np%AL(jp)
  pek => npe%AL(jp)

  ALLOCATE(pek%XRESA(pk%NSIZE_P))
  ALLOCATE(pek%XLE  (pk%NSIZE_P))  
  pek%XRESA(:) = 100.
  pek%XLE(:) = 0.
  
  IF (io%CPHOTO/='NON') THEN
    ALLOCATE(pek%XANFM  (pk%NSIZE_P))
    ALLOCATE(pek%XAN    (pk%NSIZE_P))
    ALLOCATE(pek%XANDAY (pk%NSIZE_P))
    pek%XANFM (:) = xanfminit  
    pek%XAN   (:) = 0.
    pek%XANDAY(:) = 0.      
    !
    ALLOCATE(pek%XBIOMASS         (pk%NSIZE_P,io%NNBIOMASS))
    ALLOCATE(pek%XRESP_BIOMASS    (pk%NSIZE_P,io%NNBIOMASS))
    pek%XBIOMASS(:,:) = 0.
    pek%XRESP_BIOMASS(:,:) = 0.    
  END IF

ENDDO
!
!
!* aerodynamical resistance
!
yrecfm = 'RESA'
 CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
DO jp = 1,io%NPATCH
  CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XRESA(:))
ENDDO
!
!* patch averaged radiative temperature (K)
!
ALLOCATE(s%XTSRAD_NAT(ilu))
IF (iversion<6) THEN
  s%XTSRAD_NAT(:)=0.
  DO jp=1,io%NPATCH
    DO ji = 1,np%AL(jp)%NSIZE_P
      imask = np%AL(jp)%NR_P(ji)
      s%XTSRAD_NAT(imask) = s%XTSRAD_NAT(imask)+npe%AL(jp)%XTG(ji,1)
    ENDDO
  ENDDO
  s%XTSRAD_NAT(:)=s%XTSRAD_NAT(:)/io%NPATCH
ELSE
  yrecfm='TSRAD_NAT'
  CALL read_surf(hprogram,yrecfm,s%XTSRAD_NAT(:),iresp)
ENDIF
!
DO jp = 1,io%NPATCH
  npe%AL(jp)%XLE(:) = xundef
ENDDO
!
!*       5. ISBA-AGS variables
!
IF (io%CPHOTO/='NON') THEN
  yrecfm = 'AN'
  CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XAN(:))
  ENDDO
  !
  yrecfm = 'ANDAY'
  CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XANDAY(:))
  ENDDO  
  !
  yrecfm = 'ANFM'
  CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XANFM(:))
  ENDDO  
  !
  yrecfm = 'LE_AGS'
  CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
  DO jp = 1,io%NPATCH
    CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XLE(:))
  ENDDO  
END IF
!
IF (io%CPHOTO=='NIT'.OR.io%CPHOTO=='NCB') THEN
  !
  ALLOCATE(zwork3d(ilu,io%NNBIOMASS,io%NPATCH))    
  IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) THEN
    yrecfm='BIOMA'
  ELSE
    yrecfm='BIOMASS'
  ENDIF
  CALL read_surf_layers(hprogram,yrecfm,gdim,zwork3d,iresp)
  DO jnb=1,io%NNBIOMASS
    WRITE(ylvl,'(I1)') jnb
    IF ( trim(cassim_isba)=="EKF" .AND. lprt ) THEN
      ycbio = yrecfm(:len_trim(yrecfm))//adjustl(ylvl(:len_trim(ylvl)))
      ! read in control variable
      IF ( trim(cvar(nivar)) == "LAI" .AND. trim(cbio)==trim(ycbio) ) THEN
        WHERE ( zwork3d(:,jnb,:)/=xundef ) 
          zwork3d(:,jnb,:) = zwork3d(:,jnb,:) * ( 1. + xtprt(nivar) )
        ENDWHERE
      ENDIF
    ELSEIF ( trim(cassim_isba)=="ENKF" .AND. nie<nens+1 .AND. .NOT.lassim ) THEN
      !
      IF ( trim(cbio)==trim(yrecfm) ) THEN
        DO jvar = 1,nvar
          IF (trim(cvar(jvar)) == "LAI") THEN
            DO ji = 1,ilu
              DO jp = 1,io%NPATCH
                zwork3d(ji,jnb,jp) = zwork3d(ji,jnb,jp) + xaddinfl(jvar)*random_normal()
              ENDDO
            ENDDO
            EXIT
          ENDIF
        ENDDO
      ENDIF
      !      
    ENDIF     
    DO jp = 1,io%NPATCH
      CALL pack_same_rank(np%AL(jp)%NR_P,zwork3d(:,jnb,jp),npe%AL(jp)%XBIOMASS(:,jnb))
    ENDDO 
  END DO
  DEALLOCATE(zwork3d)
!
  iwork=0
  IF(io%CPHOTO=='NCB'.OR.iversion<8)iwork=2
!
  DO jnb=2,io%NNBIOMASS-iwork
    WRITE(ylvl,'(I1)') jnb
    IF (iversion>7 .OR. (iversion==7 .AND. ibugfix>=3)) THEN
      yrecfm='RESPI'//adjustl(ylvl(:len_trim(ylvl)))
    ELSE
      yrecfm='RESP_BIOM'//adjustl(ylvl(:len_trim(ylvl)))
    ENDIF    
    CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
    IF ( trim(cassim_isba)=="EKF" .AND. lprt ) THEN
      ! read in control variable
      IF ( trim(cvar(nivar)) == "LAI" .AND. trim(cbio)==trim(yrecfm) ) THEN
        WHERE ( zwork(:,:)/=xundef ) 
          zwork(:,:) = zwork(:,:) + xtprt(nivar)*zwork(:,:)
        ENDWHERE
    ELSEIF ( trim(cassim_isba)=="ENKF" .AND. nie<nens+1 .AND. .NOT.lassim ) THEN
      !
      IF ( trim(cbio)==trim(yrecfm) ) THEN
        DO jvar = 1,nvar
          IF (trim(cvar(jvar)) == "LAI") THEN
            DO ji = 1,ilu
              DO jp = 1,io%NPATCH
                zwork(ji,jp) = zwork(ji,jp) + xaddinfl(jvar)*random_normal()
              ENDDO
            ENDDO
            EXIT
          ENDIF
        ENDDO
      ENDIF
      !  
      ENDIF
    ENDIF  
    DO jp = 1,io%NPATCH
      CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XRESP_BIOMASS(:,jnb))
    ENDDO
  END DO
  !
ENDIF
!
DEALLOCATE(zcofswi)
!
!*       6. Soil carbon
!
!
IF (io%CRESPSL=='CNT') THEN
  !
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    pek => npe%AL(jp)
    ALLOCATE(pek%XLITTER      (pk%NSIZE_P,io%NNLITTER,io%NNLITTLEVS))
    ALLOCATE(pek%XSOILCARB    (pk%NSIZE_P,io%NNSOILCARB))
    ALLOCATE(pek%XLIGNIN_STRUC(pk%NSIZE_P,io%NNLITTLEVS))
    !
    pek%XLITTER(:,:,:) = 0.
    pek%XSOILCARB(:,:) = 0. 
    pek%XLIGNIN_STRUC(:,:) = 0.  
    !
  ENDDO
  !
  DO jnlitter=1,io%NNLITTER
    DO jnlittlevs=1,io%NNLITTLEVS
      WRITE(ylvl,'(I1,A1,I1)') jnlitter,'_',jnlittlevs
      yrecfm='LITTER'//adjustl(ylvl(:len_trim(ylvl)))
      CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
      DO jp = 1,io%NPATCH
        CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XLITTER(:,jnlitter,jnlittlevs))
      ENDDO       
    END DO
  END DO

  DO jnsoilcarb=1,io%NNSOILCARB
    WRITE(ylvl,'(I4)') jnsoilcarb
    yrecfm='SOILCARB'//adjustl(ylvl(:len_trim(ylvl)))
    CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
    DO jp = 1,io%NPATCH
      CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XSOILCARB(:,jnsoilcarb))
    ENDDO      
  END DO
!
  DO jnlittlevs=1,io%NNLITTLEVS
    WRITE(ylvl,'(I4)') jnlittlevs
    yrecfm='LIGN_STR'//adjustl(ylvl(:len_trim(ylvl)))
    CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
    DO jp = 1,io%NPATCH
      CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),npe%AL(jp)%XSOILCARB(:,jnlittlevs))
    ENDDO     
  END DO
!
ENDIF

IF ( lassim ) THEN
  IF ( trim(cassim_isba) == "OI" ) THEN
    IF ( io%NPATCH /= 1 ) CALL abor1_sfx ('Reading of diagnostical values for'&
                       & //'assimilation at the moment only works for one patch for OI')          
    ! Diagnostic fields for assimilation
    IF ( .NOT. ALLOCATED(xat2m_isba)) ALLOCATE(xat2m_isba(ilu,1))
    xat2m_isba=xundef
    yrecfm='T2M'
    CALL io_buff(yrecfm,'R',gknown)
    CALL read_surf(hprogram,yrecfm,xat2m_isba(:,1),iresp)

    IF ( .NOT. ALLOCATED(xahu2m_isba)) ALLOCATE(xahu2m_isba(ilu,1))
    xahu2m_isba=xundef
    yrecfm='HU2M'
    CALL io_buff(yrecfm,'R',gknown)
    CALL read_surf(hprogram,yrecfm,xahu2m_isba(:,1),iresp)

    IF ( .NOT. ALLOCATED(xazon10m_isba)) ALLOCATE(xazon10m_isba(ilu,1))
    xazon10m_isba=xundef
    yrecfm='ZON10M'
    CALL io_buff(yrecfm,'R',gknown)
    CALL read_surf(hprogram,yrecfm,xazon10m_isba(:,1),iresp)

    IF ( .NOT. ALLOCATED(xamer10m_isba)) ALLOCATE(xamer10m_isba(ilu,1))
    xamer10m_isba=xundef
    yrecfm='MER10M'
    CALL io_buff(yrecfm,'R',gknown)
    CALL read_surf(hprogram,yrecfm,xamer10m_isba(:,1),iresp)
  ELSEIF ( nific/=nvar+2 ) THEN
    ! Diagnostic fields for EKF assimilation ("observations")
    DO iobs = 1,nobstype
     SELECT CASE (trim(cobs(iobs)))
       CASE("T2M")
         IF ( .NOT. ALLOCATED(xat2m_isba)) ALLOCATE(xat2m_isba(ilu,1))
         xat2m_isba=xundef
         yrecfm='T2M'
         CALL io_buff(yrecfm,'R',gknown)
         CALL read_surf(hprogram,yrecfm,xat2m_isba(:,1),iresp)
       CASE("HU2M")
         IF ( .NOT. ALLOCATED(xahu2m_isba)) ALLOCATE(xahu2m_isba(ilu,1))
         xahu2m_isba=xundef
         yrecfm='HU2M'
         CALL io_buff(yrecfm,'R',gknown)
         CALL read_surf(hprogram,yrecfm,xahu2m_isba(:,1),iresp)
       CASE("WG1")
         ! This is already read above
       CASE("WG2")
         ! This is already read above
       CASE("LAI")
         ! This is already read above   
       CASE("SWE")
         ! This is handled independently 
       CASE DEFAULT
         CALL abor1_sfx("Mapping of "//trim(cobs(iobs))//" is not defined in READ_ISBA_n!")
     END SELECT
    ENDDO
  ENDIF
ENDIF
!
DEALLOCATE(zwork)
!
DO jp = 1,io%NPATCH
  pk => np%AL(jp)
  pek => npe%AL(jp)
  ALLOCATE(pek%XSNOWFREE_ALB     (pk%NSIZE_P))
  ALLOCATE(pek%XSNOWFREE_ALB_VEG (pk%NSIZE_P))
  ALLOCATE(pek%XSNOWFREE_ALB_SOIL(pk%NSIZE_P))
ENDDO
!
IF (lhook) CALL dr_hook('READ_ISBA_N',1,zhook_handle)
!
CONTAINS
!
SUBROUTINE make_ens_enkf(KWORK,KLU,HREC,PCOFSWI,NP,PVAR)
!
USE modd_assim, ONLY : lens_gen, xaddtimecorr, xaddinfl, xassim_winh
!
USE modi_add_noise
USE mode_random
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KWORK
INTEGER, INTENT(IN) :: KLU
 CHARACTER(LEN=3), INTENT(IN) :: HREC
REAL, DIMENSION(:), INTENT(IN) :: PCOFSWI
TYPE(isba_np_t), INTENT(INOUT) :: NP
REAL, DIMENSION(:,:,:), INTENT(INOUT), OPTIONAL :: PVAR
!
 CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
 CHARACTER(LEN=4) :: YLVL
 CHARACTER(LEN=3) :: YVAR
REAL, DIMENSION(KLU) :: ZVAR
REAL :: ZWHITE_NOISE, ZVAR0
INTEGER :: JL, JI, JP, IVAR
LOGICAL :: GPASS
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('READ_ISBA_N:MAKE_ENS_ENKF',0,zhook_handle)
!
!
DO jl=1,kwork
  !
  IF (kwork>1) THEN
    WRITE(ylvl,'(I4)') jl
    yrecfm = trim(hrec)//adjustl(ylvl(:len_trim(ylvl)))
  ELSE
    yrecfm = trim(hrec)
  ENDIF
  !
  ivar = 0
  DO jvar = 1,nvar
    gpass = ( trim(cvar(jvar))==trim(yrecfm) )
    IF (gpass) THEN
      ivar = jvar
      EXIT
    ENDIF
  ENDDO
  !
  IF ( gpass ) THEN
    !
    IF (xaddinfl(ivar)>0.) THEN
      !
      IF (lassim .OR. (.NOT.lens_gen .AND. xaddtimecorr(ivar)>0.)) THEN
        !
        WRITE(yvar,'(I3)') ivar
        yrecfm='RD_NS'//adjustl(yvar(:len_trim(yvar)))
        CALL make_choice_array(hprogram, io%NPATCH, gdim, yrecfm, zwork)
        DO jp = 1,io%NPATCH
          CALL pack_same_rank(np%AL(jp)%NR_P,zwork(:,jp),np%AL(jp)%XRED_NOISE(:,ivar))
        ENDDO           
        !
        IF (.NOT.lassim) THEN
          !
          DO jp = 1,io%NPATCH
            pk => np%AL(jp)
            DO ji = 1,np%AL(jp)%NSIZE_P
              imask = pk%NR_P(ji)
              zwhite_noise = xaddinfl(ivar)*pcofswi(imask)*random_normal()
              CALL add_noise(xaddtimecorr(ivar),xassim_winh,zwhite_noise,pk%XRED_NOISE(ji,ivar))
            ENDDO
          ENDDO
          !
          zcoef = xassim_winh/24.
          !
        ENDIF
        !
      ELSE
        !
        DO jp = 1,io%NPATCH
          pk => np%AL(jp)
          DO ji = 1,np%AL(jp)%NSIZE_P
            imask = pk%NR_P(ji)
            np%AL(jp)%XRED_NOISE(ji,ivar) = xaddinfl(ivar)*pcofswi(imask)*random_normal()
          ENDDO
        ENDDO
        !
        zcoef = 1. 
        !
      ENDIF
      !
      IF (.NOT.lassim) THEN
        !
        DO jp = 1,io%NPATCH
          !
          zvar(:) = 0.
          IF (trim(hrec)=='TG') THEN
            zvar(1:np%AL(jp)%NSIZE_P) = npe%AL(jp)%XTG(:,jl)
          ELSEIF (trim(hrec)=='WG') THEN
            zvar(1:np%AL(jp)%NSIZE_P) = npe%AL(jp)%XWG(:,jl)
          ELSEIF (trim(hrec)=='LAI' .AND. PRESENT(pvar)) THEN
            zvar(1:np%AL(jp)%NSIZE_P) = pvar(1:np%AL(jp)%NSIZE_P,jl,jp)
          ELSE
            CALL abor1_sfx("READ_ISBAn: HREC "//hrec//" not permitted")
          ENDIF
          !
          DO ji = 1,np%AL(jp)%NSIZE_P
            IF ( zvar(ji)/=xundef ) THEN
              !
              zvar0 = zvar(ji)
              !
              zvar(ji) = zvar(ji) + zcoef * np%AL(jp)%XRED_NOISE(ji,ivar)
              !
              IF (zvar(ji) < 0.) THEN
                IF (lens_gen) THEN
                  zvar(ji) = abs(zvar(ji))
                ELSE
                  zvar(ji) = zvar0
                ENDIF
              ENDIF
            ENDIF
          ENDDO
          !
          IF (trim(hrec)=='TG') THEN
            npe%AL(jp)%XTG(:,jl) = zvar(1:np%AL(jp)%NSIZE_P)
          ELSEIF (trim(hrec)=='WG') THEN
            npe%AL(jp)%XWG(:,jl) = zvar(1:np%AL(jp)%NSIZE_P)
          ELSEIF (trim(hrec)=='LAI') THEN
            pvar(1:np%AL(jp)%NSIZE_P,jl,jp) = zvar(1:np%AL(jp)%NSIZE_P)
          ENDIF
          !
        ENDDO
        !
      ENDIF
      !
    ENDIF
    !
  ENDIF
  !
ENDDO
!
IF (lhook) CALL dr_hook('READ_ISBA_N:MAKE_ENS_ENKF',1,zhook_handle)
!
END SUBROUTINE make_ens_enkf
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE read_isba_n