! radsurf_properties.f90 - Derived type for surface properties ! (C) Copyright 2017- ECMWF. ! This software is licensed under the terms of the Apache Licence Version 2.0 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. ! In applying this licence, ECMWF does not waive the privileges and immunities ! granted to it by virtue of its status as an intergovernmental organisation ! nor does it submit to any jurisdiction. ! Author: Robin Hogan ! Email: r.j.hogan@ecmwf.int ! module radsurf_properties use parkind1, only : jpim, jprb implicit none public ! Number of tile types integer(kind=jpim), parameter :: NTileTypes = 3 ! Codes for the different type of tile enum, bind(c) enumerator :: ITileFlat = 1, & & ITileVegetation, & & ITileUrban3D end enum character(len=*), parameter :: TileRepresentationName(NTileTypes) & & = (/ 'Flat ', & & 'HomogeneousVegetation', & & 'Urban3D ' /) ! Number of surface facets and regions integer(kind=jpim), parameter :: NTileFacets (NTileTypes) = (/ 1, 1, 3 /) integer(kind=jpim), parameter :: NTileRegions(NTileTypes) = (/ 0, 1, 1 /) !--------------------------------------------------------------------- ! Derived type storing a physical description of the properties of ! the surface tiles needed to compute the boundary conditions for ! the radiation scheme. type surface_type ! Skin temperature (ncol,nfacet) real(kind=jprb), allocatable :: skin_temperature(:,:) ! Air temperature in canopy (ncol,ntile) real(kind=jprb), allocatable :: canopy_temperature(:,:) ! Shortwave albedo: if sw_albedo_direct is not allocated then ! sw_albedo will be used for both direct and diffuse solar ! radiation; otherwise, sw_albedo will be used for diffuse ! radiation and sw_albedo_direct for direct ! radiation. Dimensioning is (ncol,nalbedobands,nfacet). real(kind=jprb), allocatable, dimension(:,:,:) :: & & sw_albedo, sw_albedo_direct ! Longwave emissivity: dimensioning is (ncol,nemissbands,nfacet) real(kind=jprb), allocatable, dimension(:,:,:) :: & & lw_emissivity ! Representation codes (ITileFlat etc) for each tile: dimensioning ! is (ntile) integer(kind=jpim), allocatable, dimension(:) :: & & i_representation ! Indices to the facets: dimensioning is (ncol,ntile) integer(kind=jpim), allocatable, dimension(:,:) :: & & i_ground_facet, & & i_roof_facet, & & i_wall_facet ! Indices to the regions: dimensioning is (ncol,ntile) integer(kind=jpim), allocatable, dimension(:,:) :: & & i_region_1 !, i_region_2 ! Physical properties dimensioned (ncol,ntile) real(kind=jprb), allocatable, dimension(:,:) :: & & tile_fraction, & ! Fraction of column occupied by each tile & canopy_depth ! Depth of canopy (m) ! Building properties in urban tile dimensioned (ncol,ntile) real(kind=jprb), allocatable, dimension(:,:) :: & & building_fraction, & ! Fraction of canopy containing buildings & building_normalized_perimeter ! Perimeter length of buildings per unit area (m) ! Vegetation physical properties dimensioned (ncol,ntile) real(kind=jprb), allocatable, dimension(:,:) :: & & vegetation_optical_depth, & ! Same as LAI? & vegetation_fractional_std ! Fractional standard deviation ! Vegetation optical properties real(kind=jprb), allocatable, dimension(:,:,:) :: & & vegetation_sw_albedo, & ! (ncol,nalbedobands,ntile) & vegetation_lw_emissivity! (ncol,nemissbands,ntile) ! When 3D effects in vegetation canopies are represented, these will ! be available ! & vegetation_normalized_perimeter ! & vegetation_fraction, & ! Number of tiles, columns, facets, regions integer(kind=jpim) :: ntile, ncol, nfacet, nregion ! Number of albedo and emissivity bands integer(kind=jpim) :: nalbedobands, nemissbands ! Is this a simple surface containing one "flat" tile? logical :: is_simple = .true. contains procedure :: allocate => allocate_surface procedure :: deallocate => deallocate_surface procedure :: read => read_from_netcdf procedure :: set_facet_indices end type surface_type contains !--------------------------------------------------------------------- subroutine allocate_surface(this, ncol, nalbedobands, nemissbands, & & i_representation, & & use_sw_albedo_direct) use yomhook, only : lhook, dr_hook class(surface_type), intent(inout) :: this integer(kind=jpim), intent(in) :: ncol, nalbedobands, nemissbands integer(kind=jpim), intent(in), optional :: i_representation(:) logical(kind=jpim), intent(in), optional :: use_sw_albedo_direct ! Number of facets and regions integer(kind=jpim) :: nfacet, nregion, ntile real(jprb) :: hook_handle if (lhook) call dr_hook('radiation_surface:allocate',0,hook_handle) call this%deallocate() this%ncol = ncol if (present(i_representation)) then nfacet = sum(NTileFacets (i_representation)) nregion = sum(NTileRegions(i_representation)) ntile = size(i_representation) this%is_simple = .false. else nfacet = 1 nregion = 0 ntile = 1 this%is_simple = .true. end if this%nfacet = nfacet this%nregion= nregion this%ntile = ntile this%nalbedobands = nalbedobands this%nemissbands = nemissbands allocate(this%skin_temperature(ncol,nfacet)) allocate(this%sw_albedo(ncol, nalbedobands, nfacet)) if (present(use_sw_albedo_direct)) then if (use_sw_albedo_direct) then allocate(this%sw_albedo_direct(ncol, nalbedobands, nfacet)) end if end if allocate(this%lw_emissivity(ncol, nemissbands, nfacet)) if (this%ntile > 1) then allocate(this%tile_fraction(ncol,ntile)) end if if (nregion > 0) then allocate(this%canopy_depth(ncol,ntile)) allocate(this%canopy_temperature(ncol,ntile)) this%canopy_depth = 0.0_jprb this%canopy_temperature = -1.0_jprb end if if (.not. this%is_simple) then allocate(this%i_representation(ntile)) this%i_representation = i_representation if (any(i_representation == ITileUrban3D)) then allocate(this%building_fraction(ncol,ntile)) allocate(this%building_normalized_perimeter(ncol,ntile)) ! Initialize to default values this%building_fraction = 0.0_jprb this%building_normalized_perimeter = 0.0_jprb end if if (any(i_representation == ITileVegetation)) then allocate(this%vegetation_optical_depth(ncol,ntile)) allocate(this%vegetation_fractional_std(ncol,ntile)) allocate(this%vegetation_sw_albedo(ncol,nalbedobands,ntile)) allocate(this%vegetation_lw_emissivity(ncol,nemissbands,ntile)) ! Initialize to default values this%vegetation_optical_depth = 0.0_jprb this%vegetation_fractional_std = 0.0_jprb this%vegetation_sw_albedo = 0.0_jprb this%vegetation_lw_emissivity = 1.0_jprb end if call this%set_facet_indices end if if (lhook) call dr_hook('radiation_surface:allocate',1,hook_handle) end subroutine allocate_surface !--------------------------------------------------------------------- ! Set the indices to facets subroutine set_facet_indices(this) use radiation_io, only : nulerr, radiation_abort class(surface_type), intent(inout) :: this ! Indices to tiles and facets integer(kind=jpim) :: ifacet, iregion, jtile if (.not. this%is_simple) then allocate(this%i_ground_facet(this%ncol,this%ntile)) this%i_ground_facet = 0 if (any(this%i_representation == ITileUrban3D)) then allocate(this%i_roof_facet(this%ncol,this%ntile)) allocate(this%i_wall_facet(this%ncol,this%ntile)) ! Initialize to default values this%i_roof_facet = 0 this%i_wall_facet = 0 end if if (this%nregion > 0) then allocate(this%i_region_1(this%ncol,this%ntile)) this%i_region_1 = 0 end if ! Set the indices to the various facets ifacet = 1 iregion = 1 DO jtile = 1,this%ntile this%i_ground_facet(:,jtile) = ifacet ifacet = ifacet + 1 if (this%i_representation(jtile) == ITileVegetation) then this%i_region_1(:,jtile) = iregion iregion = iregion+1 else if (this%i_representation(jtile) == ITileUrban3D) then this%i_roof_facet(:,jtile) = ifacet this%i_wall_facet(:,jtile) = ifacet+1 ifacet = ifacet+2 this%i_region_1(:,jtile) = iregion iregion = iregion+1 else if (this%i_representation(jtile) /= ITileFlat) then write(nulerr,'(a,i0,a)') '*** Error: tile representation ', & & this%i_representation(jtile), ' not understood' call radiation_abort() end if end do end if end subroutine set_facet_indices !--------------------------------------------------------------------- subroutine deallocate_surface(this) use yomhook, only : lhook, dr_hook class(surface_type), intent(inout) :: this real(jprb) :: hook_handle if (lhook) call dr_hook('radiation_surface:deallocate',0,hook_handle) if (allocated(this%skin_temperature)) then deallocate(this%skin_temperature) end if if (allocated(this%sw_albedo)) then deallocate(this%sw_albedo) end if if (allocated(this%sw_albedo_direct)) then deallocate(this%sw_albedo_direct) end if if (allocated(this%lw_emissivity)) then deallocate(this%lw_emissivity) end if if (allocated(this%i_representation)) then deallocate(this%i_representation) end if if (allocated(this%i_ground_facet)) then deallocate(this%i_ground_facet) end if if (allocated(this%i_roof_facet)) then deallocate(this%i_roof_facet) end if if (allocated(this%i_wall_facet)) then deallocate(this%i_wall_facet) end if if (allocated(this%tile_fraction)) then deallocate(this%tile_fraction) end if if (allocated(this%canopy_depth)) then deallocate(this%canopy_depth) end if if (allocated(this%canopy_temperature)) then deallocate(this%canopy_temperature) end if if (allocated(this%building_fraction)) then deallocate(this%building_fraction) end if if (allocated(this%building_normalized_perimeter)) then deallocate(this%building_normalized_perimeter) end if if (allocated(this%vegetation_optical_depth)) then deallocate(this%vegetation_optical_depth) end if if (allocated(this%vegetation_fractional_std)) then deallocate(this%vegetation_fractional_std) end if if (allocated(this%vegetation_sw_albedo)) then deallocate(this%vegetation_sw_albedo) end if if (allocated(this%vegetation_lw_emissivity)) then deallocate(this%vegetation_lw_emissivity) end if this%ntile = 0 this%ncol = 0 this%nfacet = 0 if (lhook) call dr_hook('radiation_surface:deallocate',1,hook_handle) end subroutine deallocate_surface !--------------------------------------------------------------------- ! Print a description of the surface tile types subroutine print_surface_representation(i_representation) use radiation_io, only : nulout integer(kind=jpim), dimension(:), allocatable, intent(in) :: i_representation integer :: ntile, jtile, ifacet, iregion write(nulout,'(a)') 'Surface tile representation:' if (.not. allocated(i_representation)) then write(nulout,'(a)') ' Simple (one flat tile)' else ntile = size(i_representation,1) ifacet = 1 iregion = 1; DO jtile = 1,ntile write(nulout,'(a,i0,a,a)') ' Tile ', jtile, ': ', trim(TileRepresentationName(i_representation(jtile))) write(nulout,'(a,i0,a)') ' Facet ', ifacet, ': ground' ifacet = ifacet+1 select case (i_representation(jtile)) case (ITileVegetation) write(nulout,'(a,i0,a)') ' Region ', iregion, ': vegetation canopy' iregion = iregion+1 case (ITileUrban3D) write(nulout,'(a,i0,a)') ' Facet ', ifacet, ': roof' write(nulout,'(a,i0,a)') ' Facet ', ifacet+1, ': wall' ifacet = ifacet+1 write(nulout,'(a,i0,a)') ' Region ', iregion, ': street canyon' iregion = iregion+1 end select end do end if end subroutine print_surface_representation !--------------------------------------------------------------------- subroutine read_from_netcdf(this, file) use parkind1, only : jprb, jpim use easy_netcdf, only : netcdf_file implicit none type(netcdf_file), intent(in) :: file class(surface_type), intent(inout) :: this real(kind=jprb), allocatable :: data_1d(:) integer(kind=jpim), parameter :: ipermute(3) = [2,3,1] call this%deallocate this%is_simple = .false. call file%get('skin_temperature', this%skin_temperature, do_transp=.true.) call file%get('canopy_temperature', this%canopy_temperature, do_transp=.true.) call file%get('sw_albedo', this%sw_albedo, ipermute=ipermute) if (file%exists('sw_albedo_direct')) then call file%get('sw_albedo', this%sw_albedo_direct, ipermute=ipermute) end if call file%get('lw_emissivity', this%lw_emissivity, ipermute=ipermute) call file%get('tile_representation', data_1d) this%ntile = size(data_1d) allocate(this%i_representation(this%ntile)) this%i_representation = int(data_1d) call file%get('tile_fraction', this%tile_fraction, do_transp=.true.) call file%get('canopy_depth', this%canopy_depth, do_transp=.true.) call file%get('building_fraction', this%building_fraction, do_transp=.true.) if (file%exists('building_normalized_perimeter')) then call file%get('building_normalized_perimeter', & & this%building_normalized_perimeter, do_transp=.true.) else ! Convert building scale to normalized perimeter call file%get('building_scale', this%building_normalized_perimeter, do_transp=.true.) this%building_normalized_perimeter = 4.0_jprb * this%building_fraction & & * (1.0_jprb-this%building_fraction) / max(1.0e-8_jprb,this%building_normalized_perimeter) end if call file%get('vegetation_optical_depth', this%vegetation_optical_depth, do_transp=.true.) call file%get('vegetation_sw_albedo', this%vegetation_sw_albedo, ipermute=ipermute) call file%get('vegetation_lw_emissivity', this%vegetation_lw_emissivity, ipermute=ipermute) this%nregion = sum(NTileRegions(this%i_representation)) this%ncol = size(this%skin_temperature,1) this%nfacet = size(this%skin_temperature,2) this%nalbedobands = size(this%sw_albedo,2) this%nemissbands = size(this%lw_emissivity,2) call this%set_facet_indices end subroutine read_from_netcdf end module radsurf_properties