prep_teb_extern.F90

!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 PREP_TEB_EXTERN (DTCO, GCP, TOP, BOP, &
                            HPROGRAM,HSURF,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,KLUOUT,KPATCH,PFIELD)
!     #################################################################################
!
!
!!    MODIFICATIONS
!!    -------------
!
!
USE MODD_DATA_COVER_n, ONLY : DATA_COVER_t
USE MODD_GRID_CONF_PROJ_n, ONLY : GRID_CONF_PROJ_t
USE MODD_TEB_OPTION_n, ONLY : TEB_OPTIONS_t
USE MODD_BEM_OPTION_n, ONLY : BEM_OPTIONS_t
!
USE MODD_SURFEX_MPI, ONLY : NRANK, NPIO
USE MODD_TYPE_DATE_SURF
USE MODD_SURF_PAR, ONLY: NFILENAMELGTMAX
!
USE MODI_PREP_GRID_EXTERN
USE MODI_READ_SURF
USE MODI_GET_TEB_DEPTHS
USE MODI_INTERP_GRID
USE MODI_OPEN_AUX_IO_SURF
USE MODI_CLOSE_AUX_IO_SURF
USE MODI_TOWN_PRESENCE
USE MODI_READ_TEB_PATCH
USE MODI_MAKE_CHOICE_ARRAY
!
USE MODD_PREP,       ONLY : CINGRID_TYPE, CINTERP_TYPE
USE MODD_PREP_TEB,   ONLY : XGRID_ROAD, XGRID_WALL, XGRID_ROOF, &
                            XGRID_FLOOR, XGRID_MASS, &
                            XTI_BLD_DEF, XWS_ROOF_DEF, XWS_ROAD_DEF
USE MODD_SURF_PAR, ONLY: XUNDEF, LEN_HREC
!
USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
USE PARKIND1  ,ONLY : JPRB
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(DATA_COVER_t), INTENT(INOUT) :: DTCO
TYPE(GRID_CONF_PROJ_t),INTENT(INOUT) :: GCP
TYPE(TEB_OPTIONS_t), INTENT(INOUT) :: TOP
TYPE(BEM_OPTIONS_t), INTENT(INOUT) :: BOP
!
 CHARACTER(LEN=6),   INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=7),   INTENT(IN)  :: HSURF     ! type of field
 CHARACTER(LEN=NFILENAMELGTMAX),  INTENT(IN)  :: HFILE     ! name of file
 CHARACTER(LEN=6),                INTENT(IN)  :: HFILETYPE ! type of input file
 CHARACTER(LEN=NFILENAMELGTMAX),  INTENT(IN)  :: HFILEPGD     ! name of file
 CHARACTER(LEN=6),                INTENT(IN)  :: HFILEPGDTYPE ! type of input file
INTEGER,            INTENT(IN)  :: KLUOUT    ! logical unit of output listing
INTEGER,            INTENT(IN)  :: KPATCH
REAL,DIMENSION(:,:), POINTER    :: PFIELD    ! field to interpolate horizontally
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZFIELD         ! field read
REAL, DIMENSION(:,:), ALLOCATABLE :: ZDEPTH         ! depth of each layer
REAL :: ZDEPTH_TOT     ! total depth of surface
!
REAL, DIMENSION(:,:),   ALLOCATABLE :: ZD  ! intermediate array
!
REAL, DIMENSION(:), ALLOCATABLE :: ZMASK
!
 CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
INTEGER           :: IRESP          ! reading return code
INTEGER           :: ILAYER         ! number of layers
INTEGER           :: JLAYER, JI         ! loop counter
INTEGER           :: IVERSION_PGD, IVERSION_PREP       ! SURFEX version
INTEGER           :: IBUGFIX_PGD, IBUGFIX_PREP        ! SURFEX bug version
LOGICAL           :: GOLD_NAME      ! old name flag for temperatures
LOGICAL           :: GOLD_NAME2     ! old name flag for temperatures
 CHARACTER(LEN=4)  :: YWALL_OPT      ! option of walls
 CHARACTER(LEN=6)  :: YSURF          ! Surface type
 CHARACTER(LEN=3)  :: YBEM ! key of the building energy model DEF for DEFault (Masson et al. 2002) ,
                          ! BEM for Building Energy Model (Bueno et al. 2012)
!
INTEGER           :: INI            ! total 1D dimension
!
LOGICAL :: GDIM
LOGICAL                              :: GTEB      ! flag if TEB fields are present
LOGICAL                              :: GINTERP   ! flag if TEB fields are interpolated from a profile from TEB input
INTEGER                              :: IPATCH    ! number of soil temperature patches
INTEGER                              :: ITEB_PATCH! number of TEB patches in file
 CHARACTER(LEN=3)                     :: YPATCH    ! indentificator for TEB patch
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------------
!
!*      1.     Preparation of IO for reading in the file
!              -----------------------------------------
!
!* Note that all points are read, even those without physical meaning.
!  These points will not be used during the horizontal interpolation step.
!  Their value must be defined as XUNDEF.
!
IF (LHOOK) CALL DR_HOOK('PREP_TEB_EXTERN',0,ZHOOK_HANDLE)
!
!
 CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'FULL  ')
 CALL READ_SURF(HFILETYPE,'VERSION',IVERSION_PREP,IRESP,HDIR='-')
 CALL READ_SURF(HFILETYPE,'BUG',IBUGFIX_PREP,IRESP,HDIR='-')
 GDIM = (IVERSION_PREP>8 .OR. IVERSION_PREP==8 .AND. IBUGFIX_PREP>0)
 IF (GDIM) CALL READ_SURF(HFILETYPE,'SPLIT_PATCH',GDIM,IRESP)
 CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
 !
!* reading of version of the file being read
CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'FULL  ')
CALL READ_SURF(HFILEPGDTYPE,'VERSION',IVERSION_PGD,IRESP,HDIR='-')
CALL READ_SURF(HFILEPGDTYPE,'BUG',IBUGFIX_PGD,IRESP,HDIR='-')
!
!-------------------------------------------------------------------------------------
!
!*      2.     Reading of grid
!              ---------------
!
!* reads the grid
CALL PREP_GRID_EXTERN(GCP,HFILEPGDTYPE,KLUOUT,CINGRID_TYPE,CINTERP_TYPE,INI)
!
IF (NRANK/=NPIO) INI = 0
!
!* reads if TEB fields exist in the input file
 CALL TOWN_PRESENCE(HFILEPGDTYPE,GTEB,HDIR='-')
!
ALLOCATE(ZMASK(INI))
IF (IVERSION_PGD>=7) THEN
  IF (GTEB) THEN 
     YRECFM='FRAC_TOWN'
  ELSE
     YRECFM='FRAC_NATURE'
  END IF  
  CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZMASK,IRESP,HDIR='A')
ELSE
  ZMASK(:) = 1.
ENDIF
!
 CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
!
!---------------------------------------------------------------------------------------
!
!*     3.      Orography
!              ---------
!
IF (HSURF=='ZS     ') THEN
  !
  ALLOCATE(PFIELD(INI,1))
  YRECFM='ZS'
  CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'FULL  ')
  CALL READ_SURF(HFILEPGDTYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='E')
  CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
  !
  !---------------------------------------------------------------------------------------
ELSE
!---------------------------------------------------------------------------------------
!
!*     4.     TEB fields are read
!             -------------------
!
  IF (GTEB) THEN
!
    CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'TOWN  ')
    GOLD_NAME=(IVERSION_PGD<7 .OR. (IVERSION_PGD==7 .AND. IBUGFIX_PGD<3))
    GOLD_NAME2=(IVERSION_PGD<8 .OR. (IVERSION_PGD==8 .AND. IBUGFIX_PGD<2))
    IF (.NOT.GOLD_NAME.AND.GTEB) THEN
      YRECFM='BEM'
      CALL READ_SURF(HFILEPGDTYPE,YRECFM,YBEM,IRESP,HDIR='-')
    ELSE
      YBEM='DEF'
    ENDIF
    CALL READ_TEB_PATCH(HFILEPGD,HFILEPGDTYPE,IVERSION_PGD,IBUGFIX_PGD,ITEB_PATCH,HDIR='-')
    YPATCH='   '
    IF (ITEB_PATCH>1) THEN
      WRITE(YPATCH,FMT='(A,I1,A)') 'T',MIN(KPATCH,ITEB_PATCH),'_'
    END IF
    CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
!    
!---------------------------------------------------------------------------------------
    SELECT CASE(HSURF)
!---------------------------------------------------------------------------------------
!
!*     4.1    Profile of temperatures in roads, roofs or walls
!             ------------------------------------------------
!
    CASE('T_ROAD','T_BLD ','T_ROOF','T_WALLA','T_WALLB','T_FLOOR','T_MASS')

      GINTERP = .TRUE.
      
      YSURF=HSURF(1:6)
      !
      !* interpolation on the fine vertical grid
      IF (YSURF=='T_ROAD') THEN
        ALLOCATE(PFIELD(INI,SIZE(XGRID_ROAD)))
      ELSEIF (YSURF=='T_BLD ') THEN
        ALLOCATE(PFIELD(INI,SIZE(XGRID_ROAD)))
      ELSEIF (YSURF=='T_ROOF') THEN
        ALLOCATE(PFIELD(INI,SIZE(XGRID_ROOF)))
      ELSEIF (YSURF=='T_WALL') THEN
        ALLOCATE(PFIELD(INI,SIZE(XGRID_WALL)))
      ELSEIF (YSURF=='T_FLOO') THEN
        ALLOCATE(PFIELD(INI,SIZE(XGRID_FLOOR)))
      ELSEIF (YSURF=='T_MASS') THEN
        ALLOCATE(PFIELD(INI,SIZE(XGRID_FLOOR)))
      END IF

      CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'TOWN  ')
      !* reading of number of layers
      IF (YSURF=='T_ROAD' .OR. YSURF=='T_BLD') THEN
        IF (GOLD_NAME2) THEN
          YRECFM='ROAD_LAYER'
        ELSE
          YRECFM='TEB_SOIL_LAY'
        END IF
      END IF  
      IF (YSURF=='T_ROOF') YRECFM='ROOF_LAYER'
      IF (YSURF=='T_WALL') YRECFM='WALL_LAYER'
      IF (YSURF=='T_FLOO') THEN 
        IF (YBEM=='DEF') THEN
          GINTERP=.FALSE.
        ELSE
          YRECFM='FLOOR_LAYER'
        END IF
      END IF
      IF (YSURF=='T_MASS') THEN
         IF (YBEM=='DEF') THEN
          GINTERP=.FALSE.
         ELSE
            IF (GOLD_NAME2) THEN
               YRECFM='FLOOR_LAYER'
            ELSE
               YRECFM='MASS_LAYER'
            ENDIF
         END IF
      END IF    
      IF (GINTERP) THEN
        CALL READ_SURF(HFILEPGDTYPE,YRECFM,ILAYER,IRESP,HDIR='-')
      ELSE
        ILAYER=1
      END IF
      CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
      !
      !* reading of version of the file being read
      GOLD_NAME=(IVERSION_PREP<7 .OR. (IVERSION_PREP==7 .AND. IBUGFIX_PREP<3))
      !      
      CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'TOWN  ')
      !
      !* reading option for road orientation
      YWALL_OPT = 'UNIF'
      IF (YSURF =='T_WALL' .AND. .NOT. GOLD_NAME) THEN
        CALL READ_SURF(HFILETYPE,'WALL_OPT',YWALL_OPT,IRESP,HDIR='-')
      END IF
      !
      !* reading of the profile
      ALLOCATE(ZFIELD(INI,ILAYER))
      DO JLAYER=1,ILAYER
          !
          !
          !---------------------------------------------
          ! Wall temperature profiles, for two different facing walls case
          !---------------------------------------------
          IF (YSURF =='T_WALL' .AND. YWALL_OPT/='UNIF') THEN
            IF (GOLD_NAME) THEN
                WRITE(YRECFM,'(A6,I1.1)') HSURF(1:6),JLAYER   !Keeping "_" in variable name
            ELSE
          	WRITE(YRECFM,'(A1,A5,I1.1)') HSURF(1:1),HSURF(3:7),JLAYER
            END IF
          !
          !---------------------------------------------
          ! Floor or Mass temperature profiles, if no BEM present in input data
          !---------------------------------------------
          ELSEIF ((YSURF=='T_FLOO' .OR. YSURF=='T_MASS') .AND. YBEM=='DEF') THEN
              IF (GOLD_NAME2) THEN
                WRITE(YRECFM,'(A6)') 'TI_BLD'
              ELSE
                WRITE(YRECFM,'(A6)') 'TIBLD1'
              END IF
          !---------------------------------------------
          ! Floor or Mass temperature profiles, if BEM is present in input data
          !---------------------------------------------
          ELSE IF ((YSURF=='T_FLOO' .OR. YSURF=='T_MASS') .AND. YBEM=='BEM') THEN
            !
            ! Only the value for one compartment is read and then
            ! applied to all compartiments
            !
	    IF (GOLD_NAME) THEN 
              IF (YSURF=='T_FLOO') THEN 
                WRITE(YRECFM,'(A7,I1.1,A1,I1.1)') 'T_FLOO',JLAYER !A verifier eventuellement
              ELSE IF (YSURF=='T_MASS') THEN 
                WRITE(YRECFM,'(A6,I1.1,A1,I1.1)') 'T_MASS',JLAYER
              ENDIF
              
	    ELSE IF (GOLD_NAME2) THEN
              IF (YSURF=='T_FLOO') THEN 
                WRITE(YRECFM,'(A6,I1.1,A1,I1.1)') 'TFLOO',JLAYER
              ELSE IF (YSURF=='T_MASS') THEN 
                WRITE(YRECFM,'(A5,I1.1,A1,I1.1)') 'TMASS',JLAYER
              ENDIF
            ELSE
              IF (YSURF=='T_FLOO') THEN 
                WRITE(YRECFM,'(A5,I1.1,A1,I1.1)') 'TFLOO',JLAYER,'_',1
              ELSE IF (YSURF=='T_MASS') THEN 
                WRITE(YRECFM,'(A5,I1.1,A1,I1.1)') 'TMASS',JLAYER,'_',1
              ENDIF
            ENDIF
            !
          !---------------------------------------------
          ! Below building soil temperature profiles
          !---------------------------------------------
          ELSEIF (YSURF=='T_BLD') THEN
            IF (GOLD_NAME2) THEN   ! old TEB version without soil below buildings: set a mixing between road and building temp.
              IF (JLAYER.LT.2) THEN   ! layers near the building set to building internal temperature, to avoid too cold or warm top road layers
              	  IF (GOLD_NAME2) THEN
                    	WRITE(YRECFM,'(A6)') 'TI_BLD'
                  ELSE
                       	WRITE(YRECFM,'(A6)') 'TIBLD1'
                  END IF
                    	
              ELSE IF (JLAYER.LT.10) THEN
              	 
              	 IF (GOLD_NAME) THEN
              	   WRITE(YRECFM,'(A6,I1.1)') 'T_ROAD',JLAYER 
		 ELSE
                   WRITE(YRECFM,'(A5,I1.1)') 'TROAD',JLAYER 
               	 END IF
              ELSE  
            	 IF (GOLD_NAME) THEN
              	   WRITE(YRECFM,'(A6,I2.1)') 'T_ROAD',JLAYER 
		 ELSE
                   WRITE(YRECFM,'(A5,I2.1)') 'TROAD',JLAYER 
               	 END IF
              END IF
           ELSE               ! TEB version with soil temperature profile below buildings
              IF (JLAYER.LT.10) THEN
                WRITE(YRECFM,'(A5,I1.1)') 'TBLD',JLAYER 
              ELSE  
                WRITE(YRECFM,'(A5,I2.1)') 'TBLD',JLAYER  
              END IF
           END IF
          !---------------------------------------------
          ! Road temperature profiles
          !---------------------------------------------
          ELSE IF (YSURF =='T_ROAD') THEN
             IF (JLAYER.LT.10) THEN
           	 IF (GOLD_NAME) THEN
              	   WRITE(YRECFM,'(A6,I1.1)') 'T_ROAD',JLAYER 
		 ELSE
                   WRITE(YRECFM,'(A5,I1.1)') 'TROAD',JLAYER 
               	 END IF
              ELSE
            	 IF (GOLD_NAME) THEN
              	   WRITE(YRECFM,'(A6,I2.1)') 'T_ROAD',JLAYER 
		 ELSE
                   WRITE(YRECFM,'(A5,I2.1)') 'TROAD',JLAYER 
               	 END IF
             END IF
          ELSE
          !---------------------------------------------
          ! Others: Roof or uniform Wall temperature profiles
          !---------------------------------------------
		  IF (GOLD_NAME) THEN
		    WRITE(YRECFM,'(A6,I1.1)') YSURF(1:6),JLAYER            
		  ELSE
		    WRITE(YRECFM,'(A1,A4,I1.1)') YSURF(1:1),YSURF(3:6),JLAYER            
		  END IF
          
          !
          END IF
         !
        YRECFM=YPATCH//YRECFM
        YRECFM=ADJUSTL(YRECFM)
        CALL READ_SURF(HFILETYPE,YRECFM,ZFIELD(:,JLAYER),IRESP,HDIR='E')
        !
      END DO
      CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
      !
      DO JLAYER=1,SIZE(ZFIELD,2)
        WHERE (ZMASK(:)==0.) ZFIELD(:,JLAYER) = XUNDEF
      ENDDO
      !
    !-----------------------------------------------------------------------------------------------------
    IF (GINTERP) THEN ! input data of the corresponding variable is intrepolated from the input profile
    !-----------------------------------------------------------------------------------------------------
      ALLOCATE(ZD(INI,ILAYER))
      IF (YSURF=='T_ROAD') CALL GET_TEB_DEPTHS(DTCO,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,PD_ROAD=ZD,HDIR='E')
      IF (YSURF=='T_BLD ') CALL GET_TEB_DEPTHS(DTCO,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,PD_ROAD=ZD,HDIR='E')
      IF (YSURF=='T_ROOF') CALL GET_TEB_DEPTHS(DTCO,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,PD_ROOF=ZD,HDIR='E')
      IF (YSURF=='T_WALL') CALL GET_TEB_DEPTHS(DTCO,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,PD_WALL=ZD,HDIR='E')
      IF (YSURF=='T_MASS') CALL GET_TEB_DEPTHS(DTCO,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,PD_MASS=ZD,HDIR='E')
      IF (YSURF=='T_FLOO') CALL GET_TEB_DEPTHS(DTCO,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,PD_FLOOR=ZD,HDIR='E')
      !
      IF (NRANK==NPIO) THEN
        !
        !* recovers middle layer depth (from the surface)
        ALLOCATE(ZDEPTH    (INI,ILAYER))
        DO JI=1,INI
          !
          ZDEPTH    (JI,1)= ZD(JI,1)/2.
          ZDEPTH_TOT      = ZD(JI,1)
          DO JLAYER=2,ILAYER
            ZDEPTH    (JI,JLAYER) = ZDEPTH_TOT + ZD(JI,JLAYER)/2.
            ZDEPTH_TOT = ZDEPTH_TOT + ZD(JI,JLAYER)
          ENDDO
          !
          !* in case of wall or roof, normalizes by total wall or roof thickness
          IF (YSURF=='T_ROOF' .OR. YSURF=='T_WALL' .OR. YSURF == 'T_FLOO' .OR. YSURF == 'T_MASS') THEN
            DO JLAYER=1,ILAYER
              ZDEPTH(JI,JLAYER) = ZDEPTH(JI,JLAYER) / ZDEPTH_TOT
            END DO
          END IF
          !
        ENDDO
        !
        !* interpolation on the fine vertical grid
        IF (YSURF=='T_ROAD') THEN
          CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_ROAD,PFIELD)
        ELSEIF (YSURF=='T_BLD ') THEN
          CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_ROAD,PFIELD)
        ELSEIF (YSURF=='T_ROOF') THEN
          CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_ROOF,PFIELD)
        ELSEIF (YSURF=='T_WALL') THEN
          CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_WALL,PFIELD)
        ELSEIF (YSURF=='T_FLOO') THEN
          CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_FLOOR,PFIELD)
        ELSEIF (YSURF=='T_MASS') THEN
          CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_MASS,PFIELD)
        END IF
        DEALLOCATE(ZDEPTH)        
        !
      ENDIF
      !
      DEALLOCATE(ZD)
    !-----------------------------------------------------------------------------------------------------
    ELSE ! no interpolation because input data does not exist. A constant field is applied instead.
    !-----------------------------------------------------------------------------------------------------
      DO JLAYER=1,SIZE(PFIELD,2)
        PFIELD(:,JLAYER)=ZFIELD(:,1)
      END DO
    !-----------------------------------------------------------------------------------------------------
    END IF
    !-----------------------------------------------------------------------------------------------------
      !
      DEALLOCATE(ZFIELD)
!---------------------------------------------------------------------------------------
!
!*      4.2    Internal moisture and temperature
!              ---------------------------------
!
    CASE('QI_BLD ','TI_BLD ')
      ALLOCATE(PFIELD(INI,1))
      IF (HSURF=='TI_BLD' .OR. YBEM=='BEM') THEN
        IF (GOLD_NAME2) THEN
           YRECFM=HSURF
        ELSE
           IF (HSURF=='QI_BLD ') YRECFM='QIBLD1'
           IF (HSURF=='TI_BLD ') YRECFM='TIBLD1'
        END IF
        YRECFM=YPATCH//YRECFM
        YRECFM=ADJUSTL(YRECFM)
        CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'TOWN  ')
        CALL READ_SURF(HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='E')
        CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
        WHERE (ZMASK(:)==0.) PFIELD(:,1) = XUNDEF        
      ELSE
        IF (INI>0) PFIELD(:,1) = XUNDEF
      ENDIF
!
!---------------------------------------------------------------------------------------
!
!*      4.2    Other variables
!              ---------------
!
    CASE DEFAULT
      ALLOCATE(PFIELD(INI,1))
      CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'TOWN  ')
      YRECFM=HSURF
      GOLD_NAME=(IVERSION_PREP<7 .OR. (IVERSION_PREP==7 .AND. IBUGFIX_PREP<3))
      IF (HSURF=='T_CAN  ') THEN
        YRECFM='TCANYON'
        IF (GOLD_NAME) YRECFM='T_CANYON'
      ELSEIF (HSURF=='Q_CAN  ') THEN
        YRECFM='QCANYON'
        IF (GOLD_NAME) YRECFM='Q_CANYON'
      ELSEIF (HSURF=='T_WIN2 ' .OR. HSURF=='T_WIN1') THEN
        IF (YBEM=='BEM') THEN
          YRECFM=HSURF
        ELSE
          IF (GOLD_NAME2) THEN
            YRECFM='TI_BLD'
          ELSE
            YRECFM='TIBLD1'
          END IF
        ENDIF
      ELSEIF (HSURF=='VENTNIG') THEN
        !
        ! Only one compartiment is read
        ! 
        YRECFM='VENTNIG1'
        !
      ELSEIF (HSURF=='SHADVAC') THEN
        !
        ! Only one compartiment is read
        ! 
        YRECFM='SHADVAC1'
        !
      ELSEIF (HSURF=='TDEEP_T') THEN
        YRECFM='TDEEP_TEB'
        IF (GOLD_NAME2) YRECFM='TI_ROAD'
      ENDIF
      YRECFM=YPATCH//YRECFM
      YRECFM=ADJUSTL(YRECFM)
      IF (HSURF=='PSOLD ') THEN
         IF (GOLD_NAME2 .OR.  YBEM == 'DEF') THEN
            PFIELD(:,1) = XUNDEF
         ELSE
            CALL READ_SURF(HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='E')
         ENDIF
      ELSEIF (HSURF=='VENTNIG') THEN
         IF (GOLD_NAME2 .OR.  YBEM == 'DEF') THEN
            PFIELD(:,1) = 0.0
         ELSE
            CALL READ_SURF(HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='E')
         ENDIF
      ELSEIF (HSURF=='SHADVAC') THEN
         IF (GOLD_NAME2 .OR.  YBEM == 'DEF') THEN
            PFIELD(:,1) = 0.0
         ELSE
            CALL READ_SURF(HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='E')
         ENDIF
      ELSE
         CALL READ_SURF(HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='E')
      ENDIF
      CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
      WHERE (ZMASK(:)==0.) PFIELD(:,1) = XUNDEF      
!
!---------------------------------------------------------------------------------------
    END SELECT
!---------------------------------------------------------------------------------------
!
!*     5.     Subtitutes if TEB fields do not exist
!             -------------------------------------
!
  ELSE

    SELECT CASE(HSURF)

    !* temperature profiles
    CASE('T_ROAD','T_BLD ','T_ROOF','T_WALL','T_WIN1','T_CAN','TDEEP_T','T_WALLA','T_WALLB')
      YSURF=HSURF(1:6)
      !
      !* reading of the soil surface temperature
      CALL OPEN_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE,'NATURE')
      IPATCH = 0
      CALL READ_SURF(HFILEPGDTYPE,'PATCH_NUMBER',IPATCH,IRESP,HDIR='-')
      CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
      !
      ALLOCATE(ZFIELD(INI,IPATCH))
      !
      CALL OPEN_AUX_IO_SURF(HFILE,HFILETYPE,'NATURE')
      IF (YSURF=='T_FLOO' .OR. YSURF=='T_CAN ' .OR. YSURF=='TI_ROA') THEN
        CALL MAKE_CHOICE_ARRAY(HFILETYPE, IPATCH, GDIM, 'TG2', ZFIELD(:,:),HDIR='E')  
      ELSE          
        CALL MAKE_CHOICE_ARRAY(HFILETYPE, IPATCH, GDIM, 'TG1', ZFIELD(:,:),HDIR='E') 
      ENDIF
      CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
      DO JLAYER=1,SIZE(ZFIELD,2)
        WHERE (ZMASK(:)==0.) ZFIELD(:,JLAYER) = XUNDEF
      ENDDO      
      !* fills the whole temperature profile by this soil temperature
      IF (YSURF=='T_ROAD') ILAYER=SIZE(XGRID_ROAD)
      IF (YSURF=='T_BLD ') ILAYER=SIZE(XGRID_ROAD)
      IF (YSURF=='T_ROOF') ILAYER=SIZE(XGRID_ROOF)
      IF (YSURF=='T_WALL') ILAYER=SIZE(XGRID_WALL)
      IF (YSURF=='T_FLOO') ILAYER=SIZE(XGRID_FLOOR)
      IF (YSURF=='T_WIN1' .OR. YSURF=='T_CAN ' .OR. YSURF=='TI_ROA' .OR. YSURF=='TDEEP_') ILAYER=1
      ALLOCATE(PFIELD(INI,ILAYER))
      PFIELD(:,1) = ZFIELD(:,1)
      DO JLAYER=2,ILAYER
        PFIELD(:,JLAYER) = ZFIELD(:,1)
      END DO
      DEALLOCATE(ZFIELD)

    CASE('T_MASS','TI_BLD','T_WIN2','T_FLOOR')
      YSURF=HSURF(1:6)
      IF (YSURF=='T_FLOO') ILAYER = SIZE(XGRID_FLOOR)
      IF (YSURF=='T_MASS') ILAYER = SIZE(XGRID_MASS)
      IF (YSURF=='TI_BLD'.OR.YSURF=='T_WIN2') ILAYER=1
      ALLOCATE(PFIELD(INI, ILAYER))
      PFIELD(:,:) = XTI_BLD_DEF
 
    !* building moisture
    CASE('QI_BLD ')
      ALLOCATE(PFIELD(INI,1))
      PFIELD(:,1) = XUNDEF

    !* water reservoirs
    CASE('WS_ROOF','WS_ROAD')
      ALLOCATE(PFIELD(INI,1))
      IF (HSURF=='WS_ROOF') PFIELD = XWS_ROOF_DEF
      IF (HSURF=='WS_ROAD') PFIELD = XWS_ROAD_DEF
   !
   ! Robert: These values are hardcoded at the moment
   !
    CASE('PSOLD ')
      ALLOCATE(PFIELD(INI,1))
      PFIELD = XUNDEF
   !
   CASE('VENTNIG')    
      ALLOCATE(PFIELD(INI,1))
      PFIELD = 0.0
   !
   CASE('SHADVAC')    
      ALLOCATE(PFIELD(INI,1))
      PFIELD = 0.0
   !

   !* other fields
    CASE DEFAULT
      ALLOCATE(PFIELD(INI,1))
      PFIELD = 0.

    END SELECT

  END IF
!-------------------------------------------------------------------------------------
END IF
!-------------------------------------------------------------------------------------
!
DEALLOCATE(ZMASK)
!
!*      6.     End of IO
!              ---------
!
IF (LHOOK) CALL DR_HOOK('PREP_TEB_EXTERN',1,ZHOOK_HANDLE)
!
!---------------------------------------------------------------------------------------
!
END SUBROUTINE PREP_TEB_EXTERN