prep_teb_extern.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 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 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, 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
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=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 :: 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, 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
 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
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.AND.gteb) THEN 
  yrecfm='FRAC_TOWN'
  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))
    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_ROOF','T_WALLA','T_WALLB','T_FLOOR','T_MASS')
      CALL open_aux_io_surf(hfilepgd,hfilepgdtype,'TOWN  ')
      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_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,hdir='-')
      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
        !
        IF (gold_name) THEN
          WRITE(yrecfm,'(A6,I1.1)') hsurf(1:6),jlayer
        ELSE
          !
          IF (ysurf =='T_WALL' .AND. ywall_opt/='UNIF') THEN
            WRITE(yrecfm,'(A1,A5,I1.1)') hsurf(1:1),hsurf(3:7),jlayer
          ELSEIF ((ysurf=='T_FLOO' .OR. ysurf=='T_MASS') .AND. ybem=='DEF') THEN
            IF (ysurf=='T_FLOO' .AND. jlayer>1) THEN 
              WRITE(yrecfm,'(A5,I1.1)') 'TROAD',jlayer
            ELSE
              WRITE(yrecfm,'(A6)') 'TI_BLD'
            ENDIF
          ELSE
            WRITE(yrecfm,'(A1,A4,I1.1)') hsurf(1:1),hsurf(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
      !
      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_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_floor=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
          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_FLOO' .OR. ysurf=='T_MASS') THEN
          ALLOCATE(pfield(SIZE(zfield,1),SIZE(xgrid_floor)))
          CALL interp_grid(zdepth,zfield,xgrid_floor,pfield)
        END IF
        DEALLOCATE(zdepth)        
        !
      ENDIF
      !
      DEALLOCATE(zd)
      DEALLOCATE(zfield)
!---------------------------------------------------------------------------------------
!
!*      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 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
          yrecfm='TI_BLD'
        ENDIF
      ENDIF
      yrecfm=ypatch//yrecfm
      yrecfm=adjustl(yrecfm)
      CALL read_surf(hfiletype,yrecfm,pfield(:,1),iresp,hdir='E')
      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')
      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_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') 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