Newer
Older
!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
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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
!
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 ')
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),'_'
!
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)
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
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)
!
!---------------------------------------------------------------------------------------
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')
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)
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
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
!-------------------------------------------------------------------------------------
!
!* 6. End of IO
! ---------
!
IF (LHOOK) CALL DR_HOOK('PREP_TEB_EXTERN',1,ZHOOK_HANDLE)
!
!---------------------------------------------------------------------------------------
!
END SUBROUTINE PREP_TEB_EXTERN