!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.
! #########
!!
!! MODIFICATIONS
!! -------------
!! M.Moge 01/2016 using READ_SURF_FIELD2D/3D for 2D/3D surfex fields reads
!
SUBROUTINE PREP_TEB_EXTERN (DTCO,GCP, &
HPROGRAM,HSURF,HFILE,HFILETYPE,HFILEPGD,HFILEPGDTYPE,KLUOUT,KPATCH,PFIELD)
! #################################################################################
!
!
!
!
USE MODD_DATA_COVER_n, ONLY : DATA_COVER_t
USE MODD_GRID_CONF_PROJ, ONLY : GRID_CONF_PROJ_t
!
USE MODD_TYPE_DATE_SURF
!
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_READ_SURF_FIELD2D
!
USE MODD_PREP, ONLY : CINGRID_TYPE, CINTERP_TYPE
USE MODD_PREP_TEB, ONLY : XGRID_ROAD, XGRID_WALL, XGRID_ROOF, &
XGRID_FLOOR, XWS_ROOF, XWS_ROAD, &
XTI_BLD_DEF, XWS_ROOF_DEF, XWS_ROAD_DEF, XHUI_BLD_DEF
USE MODD_DATA_COVER_PAR, ONLY : JPCOVER
USE MODD_SURF_PAR, ONLY: XUNDEF
!
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
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=7), INTENT(IN) :: HSURF ! type of field
CHARACTER(LEN=28), INTENT(IN) :: HFILE ! name of file
CHARACTER(LEN=6), INTENT(IN) :: HFILETYPE ! type of input file
CHARACTER(LEN=28), 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, DIMENSION(:), ALLOCATABLE :: ZDEPTH_TOT ! total depth of surface
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZD ! intermediate array
!
REAL, DIMENSION(:), ALLOCATABLE :: ZMASK
!
CHARACTER(LEN=12) :: YRECFM ! Name of the article to be read
INTEGER :: IRESP ! reading return code
INTEGER :: ILAYER ! number of layers
INTEGER :: JLAYER ! loop counter
INTEGER :: IVERSION ! SURFEX version
INTEGER :: IBUGFIX ! SURFEX bug version
LOGICAL :: GOLD_NAME ! 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 :: GTEB ! flag if TEB fields are present
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)
!
!
!* reading of version of the file being read
CALL OPEN_AUX_IO_SURF(&
HFILEPGD,HFILEPGDTYPE,'FULL ')
CALL READ_SURF(&
HFILEPGDTYPE,'VERSION',IVERSION,IRESP)
CALL READ_SURF(&
HFILEPGDTYPE,'BUG',IBUGFIX,IRESP)
CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
GOLD_NAME=(IVERSION<7 .OR. (IVERSION==7 .AND. IBUGFIX<3))
!
!-------------------------------------------------------------------------------------
!
!* 2. Reading of grid
! ---------------
!
CALL OPEN_AUX_IO_SURF(&
HFILEPGD,HFILEPGDTYPE,'FULL ')
!* reads the grid
CALL PREP_GRID_EXTERN(GCP,&
HFILEPGDTYPE,KLUOUT,CINGRID_TYPE,CINTERP_TYPE,INI)
!
YRECFM='VERSION'
CALL READ_SURF(HFILEPGDTYPE,YRECFM,IVERSION,IRESP)
!
ALLOCATE(ZMASK(INI))
IF (IVERSION>=7) THEN
YRECFM='FRAC_TOWN'
CALL READ_SURF(HFILEPGDTYPE,YRECFM,ZMASK,IRESP,HDIR='A')
ELSE
ZMASK(:) = 1.
ENDIF
!
!* reads if TEB fields exist in the input file
CALL TOWN_PRESENCE(&
HFILEPGDTYPE,GTEB)
CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
!
IF (.NOT.GOLD_NAME.AND.GTEB) THEN
YRECFM='BEM'
CALL OPEN_AUX_IO_SURF(&
HFILEPGD,HFILEPGDTYPE,'TOWN ')
CALL READ_SURF(&
HFILEPGDTYPE,YRECFM,YBEM,IRESP)
CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
ELSE
YBEM='DEF'
ENDIF
!---------------------------------------------------------------------------------------
!
!* 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='A')
CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
!
!---------------------------------------------------------------------------------------
ELSE
!---------------------------------------------------------------------------------------
!
!* 4. TEB fields are read
! -------------------
!
IF (GTEB) THEN
!
CALL READ_TEB_PATCH(&
HFILEPGD,HFILEPGDTYPE,ITEB_PATCH)
YPATCH=' '
IF (ITEB_PATCH>1) THEN
WRITE(YPATCH,FMT='(A,I1,A)') 'T',MIN(KPATCH,ITEB_PATCH),'_'
END IF
!
CALL OPEN_AUX_IO_SURF(&
HFILEPGD,HFILEPGDTYPE,'TOWN ')
!
!---------------------------------------------------------------------------------------
SELECT CASE(HSURF)
!---------------------------------------------------------------------------------------
!
!* 4.1 Profile of temperatures in roads, roofs or walls
! ------------------------------------------------
!
CASE('T_ROAD','T_ROOF','T_WALLA','T_WALLB','T_FLOOR','T_MASS')
YSURF=HSURF(1:6)
!* reading of number of layers
IF (YSURF=='T_ROAD') YRECFM='ROAD_LAYER'
IF (YSURF=='T_ROOF') YRECFM='ROOF_LAYER'
IF (YSURF=='T_WALL') YRECFM='WALL_LAYER'
IF (YSURF=='T_WALLA') YRECFM='WALL_LAYER'
IF (YSURF=='T_WALLB') YRECFM='WALL_LAYER'
IF (YSURF=='T_FLOO' .OR. YSURF=='T_MASS') THEN
IF (YBEM=='DEF') THEN
YRECFM='ROAD_LAYER'
ELSE
YRECFM='FLOOR_LAYER'
END IF
END IF
CALL READ_SURF(&
HFILEPGDTYPE,YRECFM,ILAYER,IRESP)
CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
!
ALLOCATE(ZD(INI,ILAYER))
IF (YSURF=='T_ROAD') CALL GET_TEB_DEPTHS(&
DTCO, &
HFILEPGD,HFILEPGDTYPE,PD_ROAD=ZD)
IF (YSURF=='T_ROOF') CALL GET_TEB_DEPTHS(&
DTCO, &
HFILEPGD,HFILEPGDTYPE,PD_ROOF=ZD)
IF (YSURF=='T_WALL') CALL GET_TEB_DEPTHS(&
DTCO, &
HFILEPGD,HFILEPGDTYPE,PD_WALL=ZD)
IF (YSURF=='T_WALLA') CALL GET_TEB_DEPTHS(&
DTCO, &
HFILEPGD,HFILEPGDTYPE,PD_WALL=ZD)
IF (YSURF=='T_WALLB') CALL GET_TEB_DEPTHS(&
DTCO, &
HFILEPGD,HFILEPGDTYPE,PD_WALL=ZD)
IF (YSURF=='T_MASS') CALL GET_TEB_DEPTHS(&
DTCO, &
HFILEPGD,HFILEPGDTYPE,PD_FLOOR=ZD)
IF (YSURF=='T_FLOO') CALL GET_TEB_DEPTHS(&
DTCO, &
HFILEPGD,HFILEPGDTYPE,PD_FLOOR=ZD)
!
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)
END IF
IF (YSURF =='T_WALLA' .AND. .NOT. GOLD_NAME) THEN
CALL READ_SURF(&
HFILETYPE,'WALL_OPT',YWALL_OPT,IRESP)
END IF
IF (YSURF =='T_WALLB' .AND. .NOT. GOLD_NAME) THEN
CALL READ_SURF(&
HFILETYPE,'WALL_OPT',YWALL_OPT,IRESP)
END IF
!
!* reading of the profile
ALLOCATE(ZFIELD(INI,ILAYER))
DO JLAYER=1,ILAYER
IF (GOLD_NAME) THEN
WRITE(YRECFM,'(A6,I1.1)') HSURF(1:6),JLAYER
ELSE
WRITE(YRECFM,'(A1,A4,I1.1)') HSURF(1:1),HSURF(3:6),JLAYER
IF (YSURF =='T_WALL' .AND. YWALL_OPT/='UNIF') &
WRITE(YRECFM,'(A1,A5,I1.1)') HSURF(1:1),HSURF(3:7),JLAYER
IF (YSURF =='T_WALLA' .AND. YWALL_OPT/='UNIF') &
WRITE(YRECFM,'(A1,A5,I1.1)') HSURF(1:1),HSURF(3:7),JLAYER
IF (YSURF =='T_WALLB' .AND. YWALL_OPT/='UNIF') &
WRITE(YRECFM,'(A1,A5,I1.1)') HSURF(1:1),HSURF(3:7),JLAYER
IF ((HSURF=='T_FLOOR' .OR. HSURF=='T_MASS') .AND. YBEM=='DEF') THEN
IF (HSURF=='T_FLOOR' .AND. JLAYER>1) THEN
WRITE(YRECFM,'(A5,I1.1)') 'TROAD',JLAYER
ELSE
WRITE(YRECFM,'(A6)') 'TI_BLD'
ENDIF
END IF
END IF
YRECFM=YPATCH//YRECFM
YRECFM=ADJUSTL(YRECFM)
CALL READ_SURF(&
HFILETYPE,YRECFM,ZFIELD(:,JLAYER),IRESP,HDIR='A')
END DO
CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
DO JLAYER=1,SIZE(ZFIELD,2)
WHERE (ZMASK(:)==0.) ZFIELD(:,JLAYER) = XUNDEF
ENDDO
!
!* recovers middle layer depth (from the surface)
ALLOCATE(ZDEPTH (INI,ILAYER))
ALLOCATE(ZDEPTH_TOT(INI))
ZDEPTH (:,1)=ZD(:,1)/2.
ZDEPTH_TOT(:) =ZD(:,1)
DO JLAYER=2,ILAYER
ZDEPTH (:,JLAYER) = ZDEPTH_TOT(:) + ZD(:,JLAYER)/2.
ZDEPTH_TOT(:) = ZDEPTH_TOT(:) + ZD(:,JLAYER)
END DO
!
!* in case of wall or roof, normalizes by total wall or roof thickness
IF (YSURF=='T_ROOF' .OR. YSURF=='T_WALL' .OR. HSURF == 'T_FLOOR' &
&.OR. HSURF == 'T_MASS'.OR. YSURF=='T_WALLA' .OR. HSURF == 'T_WALLB') THEN
DO JLAYER=1,ILAYER
ZDEPTH(:,JLAYER) = ZDEPTH(:,JLAYER) / ZDEPTH_TOT(:)
END DO
END IF
!
!* interpolation on the fine vertical grid
IF (YSURF=='T_ROAD') THEN
ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_ROAD)))
CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_ROAD,PFIELD)
ELSEIF (YSURF=='T_ROOF') THEN
ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_ROOF)))
CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_ROOF,PFIELD)
ELSEIF (YSURF=='T_WALL') THEN
ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_WALL)))
CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_WALL,PFIELD)
ELSEIF (YSURF=='T_WALLA') THEN
ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_WALL)))
CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_WALL,PFIELD)
ELSEIF (YSURF=='T_WALLB') THEN
ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_WALL)))
CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_WALL,PFIELD)
ELSEIF (YSURF=='T_FLOO' .OR. YSURF=='T_MASS') THEN
ALLOCATE(PFIELD(SIZE(ZFIELD,1),SIZE(XGRID_FLOOR)))
CALL INTERP_GRID(ZDEPTH,ZFIELD,XGRID_FLOOR,PFIELD)
END IF
!
!* end
DEALLOCATE(ZD)
DEALLOCATE(ZFIELD)
DEALLOCATE(ZDEPTH)
DEALLOCATE(ZDEPTH_TOT)
!---------------------------------------------------------------------------------------
!
!* 4.2 Internal moisture
! ---------------
!
CASE('QI_BLD ')
ALLOCATE(PFIELD(INI,1))
IF (YBEM=='BEM') THEN
YRECFM='QI_BLD'
YRECFM=YPATCH//YRECFM
YRECFM=ADJUSTL(YRECFM)
CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
CALL OPEN_AUX_IO_SURF(&
HFILE,HFILETYPE,'TOWN ')
CALL READ_SURF(&
HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='A')
CALL CLOSE_AUX_IO_SURF(HFILE,HFILETYPE)
WHERE (ZMASK(:)==0.) PFIELD(:,1) = XUNDEF
ELSE
PFIELD(:,1) = XUNDEF
ENDIF
!
!---------------------------------------------------------------------------------------
!
!* 4.2 Other variables
! ---------------
!
CASE DEFAULT
ALLOCATE(PFIELD(INI,1))
YRECFM=HSURF
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
YRECFM='TI_BLD'
ENDIF
ENDIF
YRECFM=YPATCH//YRECFM
YRECFM=ADJUSTL(YRECFM)
CALL CLOSE_AUX_IO_SURF(HFILEPGD,HFILEPGDTYPE)
CALL OPEN_AUX_IO_SURF(&
HFILE,HFILETYPE,'TOWN ')
CALL READ_SURF(&
HFILETYPE,YRECFM,PFIELD(:,1),IRESP,HDIR='A')
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_ROOF','T_WALL','T_WIN1','T_FLOOR','T_CAN','TI_ROAD','T_WALLA','T_WALLB')
YSURF=HSURF(1:6)
!* reading of the soil surface temperature
CALL OPEN_AUX_IO_SURF(&
HFILEPGD,HFILEPGDTYPE,'NATURE')
CALL READ_SURF(&
HFILEPGDTYPE,'PATCH_NUMBER',IPATCH,IRESP)
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 READ_SURF_FIELD2D(&
HFILETYPE,ZFIELD(:,:),'TG2 ',HDIR='A')
ELSE
CALL READ_SURF_FIELD2D(&
HFILETYPE,ZFIELD(:,:),'TG1 ',HDIR='A')
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_ROOF') ILAYER=SIZE(XGRID_ROOF)
IF (YSURF=='T_WALL') ILAYER=SIZE(XGRID_WALL)
IF (YSURF=='T_WALLA') ILAYER=SIZE(XGRID_WALL)
IF (YSURF=='T_WALLB') ILAYER=SIZE(XGRID_WALL)
IF (YSURF=='T_FLOO') ILAYER=SIZE(XGRID_FLOOR)
IF (YSURF=='T_WIN1' .OR. YSURF=='T_CAN ' .OR. YSURF=='TI_ROA') ILAYER=1
ALLOCATE(PFIELD(INI,ILAYER))
IF (YSURF=='T_FLOO') THEN
!* sets the temperature equal to this deep soil temperature
PFIELD(:,1) = XTI_BLD_DEF
ELSE
PFIELD(:,1) = ZFIELD(:,1)
ENDIF
DO JLAYER=2,ILAYER
PFIELD(:,JLAYER) = ZFIELD(:,1)
END DO
DEALLOCATE(ZFIELD)
CASE('T_MASS','TI_BLD','T_WIN2')
YSURF=HSURF(1:6)
IF (YSURF=='T_MASS') ILAYER = SIZE(XGRID_FLOOR)
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
!* 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