!> @file calc_obstruction.f90 !------------------------------------------------------------------------------! ! This file is part of the PALM model system. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2021 Leibniz Universitaet Hannover !------------------------------------------------------------------------------! ! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: calc_obstruction.f90 4870 2021-02-08 18:09:41Z suehring $ ! Initial revision ! ! ! Authors: ! -------- ! @author Michael Schrempf (main author), ! Matthias Sühring (first major revision of the code) ! !> @todo - Code formatting !------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculation of sky-obstruction. !------------------------------------------------------------------------------! PROGRAM obstruction #if defined ( __netcdf ) USE netcdf #endif IMPLICIT NONE CHARACTER(LEN=10) :: b(3) CHARACTER(LEN=300) :: input_file CHARACTER(LEN=300) :: path_to_expo_files INTEGER:: spx, spy, counter, cc, sza, rza, azi, str_i, str_j, ii, & vector_length_x, vector_length_y, id, id_var, nc_stat, i, j, k INTEGER(KIND=1) :: bitvar INTEGER, DIMENSION(1:8):: date_time INTEGER, DIMENSION(0:99):: valy, valx INTEGER :: obs_ny = 9 !< number of zenith values INTEGER :: obs_nx = 35 !< number of azimuth values INTEGER :: obs_nz = 44 !< ? INTEGER(KIND=2), DIMENSION(:,:,:), ALLOCATABLE :: obs INTEGER(KIND=2), DIMENSION(:,:,:), ALLOCATABLE :: obsi ! INTEGER(KIND=2), DIMENSION(0:1159,0:1159,0:44):: obs ! INTEGER(KIND=2), DIMENSION(0:959,0:959,0:44):: obsi INTEGER :: nx !< number of grid points in x-direction INTEGER :: ny !< number of grid points in y-direction INTEGER :: nsza REAL :: dx !< grid spacing in x REAL :: dy !< grid spacing in y REAL :: dz !< grid spacing in z REAL :: delta = 10.0, fill !< delta of zenith and azimuth angles (namelist param) REAL:: dtor, pi = 3.141592654, t1 = 0.0, t2 = 0.0 REAL, DIMENSION(0:99):: value_y, value_x, value_z, value_zn, distance_x, distance_y ! REAL, DIMENSION(0:1159,0:1159):: topo ! REAL, DIMENSION(0:959,0:959):: topo_0 REAL, DIMENSION(:,:), ALLOCATABLE :: topo REAL, DIMENSION(:,:), ALLOCATABLE :: topo_0 REAL(kind=8), DIMENSION(0:35,0:9) :: integration REAL, DIMENSION(0:2,0:90) :: irrad REAL, DIMENSION(0:35,0:9) :: tmp REAL, DIMENSION(0:35,0:9,0:90) :: cube_rad REAL, DIMENSION(0:35,0:9,0:2) :: pfs CHARACTER(LEN=15) :: output ='Data_input_all' CALL obstruct_parin nsza = 90 dtor = pi / 180. vector_length_x = INT( 99.0 / dx ) vector_length_y = INT( 99.0 / dy ) ALLOCATE( obs(0:nx-1+(2*vector_length_x),0:ny-1+(2*vector_length_y),0:44) ) ALLOCATE( obsi(0:nx-1,0:ny-1,0:44) ) ALLOCATE( topo(0:nx-1+(2*vector_length_x),0:ny-1+(2*vector_length_y)) ) ALLOCATE( topo_0(0:nx-1,0:ny-1) ) obs(:,:,:) = 9 obsi(:,:,:) = 0 topo(:,:) = 0 call date_and_time(b(1), b(2), b(3), date_time) t1 = date_time(5)+ date_time(6)/60. + (date_time(7)+(date_time(8)/1000.))/3600. distance_x = 0.0 distance_y = 0.0 DO ii = 0, 99 distance_x(ii) = ii * dx distance_y(ii) = ii * dy ENDDO distance_x = distance_x + dx distance_y = distance_y + dy #if defined ( __netcdf ) nc_stat = NF90_OPEN( input_file, NF90_NOWRITE, id ) nc_stat = NF90_INQ_VARID( id, 'buildings_2d', id_var ) nc_stat = NF90_GET_VAR( id, id_var, topo_0, start = (/ 1, 1 /), count = (/ nx, ny /) ) nc_stat = NF90_GET_ATT( id, id_var, '_FillValue', fill ) nc_stat = NF90_CLOSE( id ) #endif ! OPEN(11, FILE='/home/matthias/UNI/MOSAIK/UV/test/Topo_input/topo_Reuter.txt', FORM='FORMATTED') ! DO rza = 0, nx-1 ! READ(11,'(960(F5.1,1X))') topo_0(:,rza) !(nx-1)-rza ! ENDDO ! CLOSE(11) topo(vector_length_x:nx-1+vector_length_x,vector_length_y:ny-1+vector_length_y) = & MERGE( topo_0, 0.0, topo_0 /= fill ) Do spy = vector_length_y, ny-1+vector_length_y DO spx = vector_length_x, nx-1+vector_length_x print*, spy, spx, topo(spx,spy) IF (topo(spx,spy) == 0.0 ) THEN ! print*, "in if" obs(spx,spy,:)=0 cc = 0 bitvar=2 DO sza = 0,90, 10 ! print*, sza DO azi = 0,359, 10 ! print*, azi value_zn = SIN((sza+0.0000001)*dtor) value_x = ANINT((distance_x(:)*SIN(azi*dtor)) *value_zn(:)) value_y = ANINT((distance_y(:)*COS(azi*dtor)) *value_zn(:)) ! print*, value_x ! print*, value_y valx = INT(value_x) valy = INT(value_y) ! print*, "valx", valx !value_z = (distance(:)*COS(sza*dtor)) value_z = COS(sza*dtor) * SQRT( value_x**2 + value_y**2) / value_zn ! print*, sza*dtor, sza !print*, value_z ! stop counter=0 DO str_i = 0, vector_length_x-1 DO str_j = 0, vector_length_y-1 IF ( topo( spx+valx(str_i), spy+valy(str_j) ) >= (value_z(str_i)*dz) ) THEN IF ( topo( spx+valx(str_i), spy+valy(str_j) ) == 0) THEN counter=counter + 0 ELSE counter=counter + 1 ENDIF ENDIF ENDDO ENDDO IF (sza == 90 ) Counter = 1 IF (counter == 0) bitvar = IBSET(bitvar,MOD(cc,8)) IF (counter > 0) bitvar = IBCLR(bitvar,MOD(cc,8)) !IF (counter .GT. 0) bitvar = IBSET(bitvar,MOD(cc,8)) IF ( MOD(cc,8) == 7 ) THEN obs(spx,spy,cc/8) = bitvar bitvar = 2 ENDIF cc = cc+1 ENDDO ENDDO ENDIF ENDDO ENDDO call date_and_time(b(1), b(2), b(3), date_time) t2 = date_time(5)+ date_time(6)/60. + (date_time(7)+(date_time(8)/1000.))/3600. obsi(:,:,:)=obs(vector_length_x:nx-1+vector_length_x,vector_length_y:ny-1+vector_length_y,:) OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Integration_solid_angles.txt', FORM='FORMATTED') DO rza = 0,9 READ(11,'(36(F10.8,2X))') integration(:,rza) ENDDO CLOSE(11) OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Irradiance_ROM.txt', FORM='FORMATTED') Do sza = 0,90 READ(11,'(3(F10.6,2X))') irrad(:,sza) ENDDO CLOSE(11) OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Projection_areas_CL0.txt', FORM='FORMATTED') DO rza = 0,9 READ(11,'(36(F10.6,2X))') tmp(:,rza) ENDDO CLOSE(11) pfs(:,:,0) = tmp(:,:) OPEN(11, FILE=TRIM( path_to_expo_files ) // 'Projection_areas_CL1.txt', FORM='FORMATTED') DO rza = 0,9 READ(11,'(36(F10.6,2X))') tmp(:,rza) ENDDO CLOSE(11) pfs(:,:,1) = tmp(:,:) OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Projection_areas_CL3.txt', FORM='FORMATTED') DO rza = 0,9 READ(11,'(36(F10.6,2X))') tmp(:,rza) ENDDO CLOSE(11) pfs(:,:,2) = tmp(:,:) OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Radiance_ROM.txt', FORM='FORMATTED') Do sza = 0,90 DO rza = 0,9 READ(11,'(36(F10.6,2X))') tmp(:,rza) ENDDO cube_rad(:,:,sza) = tmp(:,:) ENDDO CLOSE(11) !- Finalize output CALL data_output('open', output) CALL data_output('write', output) CALL data_output('close', output) !- End WRITE(*, '(A)') ' --- creation finished! :-) ---' ! open(unit=12, status='replace',file=& ! '/home/matthias/UNI/MOSAIK/UV/test/result_shading/obstruction_Berlin_reuter.bin', FORM='UNFORMATTED') ! write(12) obsi ! close(12) PRINT *,'Time_Start = ',t1 PRINT *,'Time_ENDE = ',t2 PRINT *,'elapsed time used (dateT) = ',(t2-t1)*3600. CONTAINS SUBROUTINE obstruct_parin IMPLICIT NONE INTEGER :: file_id_parin = 90 NAMELIST /uv_parin/ dx, dy, dz, nx, ny, input_file, path_to_expo_files ! !-- Open namelist file. OPEN( file_id_parin, FILE='uv_parin', STATUS='OLD', FORM='FORMATTED') ! !-- Read namelist. READ ( file_id_parin, uv_parin ) ! !-- Close namelist file. CLOSE( file_id_parin ) END SUBROUTINE obstruct_parin SUBROUTINE data_output(action,output) ! Data output to NetCDF file CHARACTER(LEN=*), INTENT(IN) :: action !< flag to steer routine (open/close/write) CHARACTER(LEN=*), INTENT(IN) :: output !< output file name prefix ! CHARACTER(LEN=10) :: output ='Data_input' INTEGER, SAVE :: do_count=0 !< counting output INTEGER, SAVE :: id_dim_time !< ID of dimension time INTEGER, SAVE :: id_dim_x !< ID of dimension x INTEGER, SAVE :: id_dim_xirr !< ID of dimension x INTEGER, SAVE :: id_dim_y !< ID of dimension y INTEGER, SAVE :: id_dim_z !< ID of dimension z INTEGER, SAVE :: id_dim_obsx !< ID of dimension x INTEGER, SAVE :: id_dim_obsy !< ID of dimension y INTEGER, SAVE :: id_dim_obsz !< ID of dimension z INTEGER, SAVE :: id_file !< ID of NetCDF file INTEGER, SAVE :: id_var_int !< ID of geopotential INTEGER, SAVE :: id_var_irr !< ID of geopotential INTEGER, SAVE :: id_var_obs !< ID of geopotential INTEGER, SAVE :: id_var_obsx !< ID of x INTEGER, SAVE :: id_var_obsy !< ID of y INTEGER, SAVE :: id_var_obsz !< ID of y INTEGER, SAVE :: id_var_pfs !< ID of geopotential INTEGER, SAVE :: id_var_rad !< ID of geopotential INTEGER, SAVE :: id_var_x !< ID of x INTEGER, SAVE :: id_var_xirr !< ID of x INTEGER, SAVE :: id_var_y !< ID of y INTEGER, SAVE :: id_var_z !< ID of y INTEGER, SAVE :: nc_stat !< status flag of NetCDF routines REAL, DIMENSION(:), ALLOCATABLE :: netcdf_data_1d !< 1D output array REAL, DIMENSION(:,:), ALLOCATABLE :: netcdf_data_2d !< 2D output array REAL, DIMENSION(:,:,:), ALLOCATABLE :: netcdf_data_3d !< 3D output array SELECT CASE (TRIM(action)) !- Initialize output file CASE ('open') !- Delete any pre-existing output file OPEN(20, FILE=TRIM(output)//'.nc') CLOSE(20, STATUS='DELETE') !- Open file #if defined ( __netcdf ) nc_stat = NF90_CREATE(TRIM(output)//'.nc', NF90_NOCLOBBER, id_file) IF (nc_stat /= NF90_NOERR) PRINT*, '+++ netcdf error' !- Write global attributes nc_stat = NF90_PUT_ATT(id_file, NF90_GLOBAL, & 'Conventions', 'COARDS') nc_stat = NF90_PUT_ATT(id_file, NF90_GLOBAL, & 'title', 'UVEM-model') !- Define time coordinate nc_stat = NF90_DEF_DIM(id_file, 'TIME', NF90_UNLIMITED, & id_dim_time) !- Define spatial coordinates nc_stat = NF90_DEF_DIM(id_file, 'AZI', obs_nx+1, id_dim_x) nc_stat = NF90_DEF_VAR(id_file, 'AZI', NF90_DOUBLE, id_dim_x, id_var_x) nc_stat = NF90_PUT_ATT(id_file, id_var_x, 'units', 'Grad') nc_stat = NF90_DEF_DIM(id_file, 'ZEN', obs_ny+1, id_dim_y) nc_stat = NF90_DEF_VAR(id_file, 'ZEN', NF90_DOUBLE, & id_dim_y, id_var_y) nc_stat = NF90_PUT_ATT(id_file, id_var_y, 'units', 'Grad') nc_stat = NF90_DEF_DIM(id_file, 'sza', nsza+1, id_dim_z) nc_stat = NF90_DEF_VAR(id_file, 'sza', NF90_DOUBLE, & id_dim_z, id_var_z) nc_stat = NF90_PUT_ATT(id_file, id_var_z, 'units', 'Grad') !- Define spatial coordinates nc_stat = NF90_DEF_DIM(id_file, 'COM', 3, id_dim_xirr) nc_stat = NF90_DEF_VAR(id_file, 'COM', NF90_DOUBLE, id_dim_xirr, id_var_xirr) nc_stat = NF90_PUT_ATT(id_file, id_var_xirr, 'units', 'Komponente') !- Define spatial coordinates for obstruction file nc_stat = NF90_DEF_DIM(id_file, 'X', nx-1+1, id_dim_obsx) nc_stat = NF90_DEF_VAR(id_file, 'X', NF90_DOUBLE, id_dim_obsx, id_var_obsx) nc_stat = NF90_PUT_ATT(id_file, id_var_obsx, 'units', 'meters') nc_stat = NF90_DEF_DIM(id_file, 'Y', ny-1+1, id_dim_obsy) nc_stat = NF90_DEF_VAR(id_file, 'Y', NF90_DOUBLE, id_dim_obsy, id_var_obsy) nc_stat = NF90_PUT_ATT(id_file, id_var_obsy, 'units', 'meters') nc_stat = NF90_DEF_DIM(id_file, 'directions', nsza+1, id_dim_obsz) nc_stat = NF90_DEF_VAR(id_file, 'directions', NF90_DOUBLE, id_dim_obsz, id_var_obsz) nc_stat = NF90_PUT_ATT(id_file, id_var_obsz, 'units', 'ibits') !- Define output arrays nc_stat = NF90_DEF_VAR(id_file, 'radiance', NF90_DOUBLE, & (/id_dim_x, id_dim_y, id_dim_z/), id_var_rad) nc_stat = NF90_PUT_ATT(id_file, id_var_rad, & 'long_name', 'vitD weighetd radiance at 300DU') nc_stat = NF90_PUT_ATT(id_file, id_var_rad, & 'short_name', 'vitD radiance') nc_stat = NF90_PUT_ATT(id_file, id_var_rad, 'units', 'W/m2 sr') nc_stat = NF90_DEF_VAR(id_file, 'int_factors', NF90_DOUBLE, & (/id_dim_x, id_dim_y/), id_var_int) nc_stat = NF90_PUT_ATT(id_file, id_var_int, & 'long_name', 'integration factors') nc_stat = NF90_PUT_ATT(id_file, id_var_int, & 'short_name', 'integration factors') ! nc_stat = NF90_PUT_ATT(id_file, id_var_int, 'units', '-') !- Define output arrays nc_stat = NF90_DEF_VAR(id_file, 'projarea', NF90_DOUBLE, & (/id_dim_x, id_dim_y, id_dim_time/), id_var_pfs) nc_stat = NF90_PUT_ATT(id_file, id_var_pfs, & 'long_name', 'projection areas of a human') nc_stat = NF90_PUT_ATT(id_file, id_var_pfs, & 'short_name', 'projection areas') nc_stat = NF90_PUT_ATT(id_file, id_var_pfs, 'units', 'm2') !- Define output arrays nc_stat = NF90_DEF_VAR(id_file, 'irradiance', NF90_DOUBLE, & (/id_dim_xirr, id_dim_z/), id_var_irr) nc_stat = NF90_PUT_ATT(id_file, id_var_irr, & 'long_name', 'global, diffuse and direct solar irradiance') nc_stat = NF90_PUT_ATT(id_file, id_var_irr, & 'short_name', 'irradiance values') nc_stat = NF90_PUT_ATT(id_file, id_var_irr, 'units', 'W/m2') !- Define output arrays for obstruction file nc_stat = NF90_DEF_VAR(id_file, 'obstruction', NF90_DOUBLE, & (/id_dim_obsx, id_dim_obsy, id_dim_obsz/), id_var_obs) nc_stat = NF90_PUT_ATT(id_file, id_var_obs, & 'long_name', 'calculated obstructions in different directions') nc_stat = NF90_PUT_ATT(id_file, id_var_obs, & 'short_name', 'obstructions') nc_stat = NF90_PUT_ATT(id_file, id_var_obs, 'units', '-') !- Leave define mode nc_stat = NF90_ENDDEF(id_file) !- Write axis to file ALLOCATE(netcdf_data_1d(0:nx)) !- x axis DO i = 0, nx netcdf_data_1d(i) = i * delta ENDDO nc_stat = NF90_PUT_VAR(id_file, id_var_x, netcdf_data_1d, & start = (/1/), count = (/nx+1/)) DO i = 0, 2 netcdf_data_1d(i) = i ENDDO nc_stat = NF90_PUT_VAR(id_file, id_var_xirr, netcdf_data_1d, & start = (/1/), count = (/3/)) !- y axis DO j = 0, ny netcdf_data_1d(j) = j * delta ENDDO nc_stat = NF90_PUT_VAR(id_file, id_var_y, netcdf_data_1d, & start = (/1/), count = (/ny+1/)) DEALLOCATE(netcdf_data_1d) !- z axis ALLOCATE(netcdf_data_1d(0:nsza)) DO k = 0, nsza netcdf_data_1d(k) = k ENDDO nc_stat = NF90_PUT_VAR(id_file, id_var_z, netcdf_data_1d, & start = (/1/), count = (/nsza+1/)) DEALLOCATE(netcdf_data_1d) !- define axis for obstruction file !- x axis ALLOCATE(netcdf_data_1d(0:nx-1)) DO i = 0, nx-1 netcdf_data_1d(i) = i * dx ENDDO nc_stat = NF90_PUT_VAR(id_file, id_var_obsx, netcdf_data_1d, & start = (/1/), count = (/nx-1+1/)) !- y axis DO j = 0, ny-1 netcdf_data_1d(j) = j * dy ENDDO nc_stat = NF90_PUT_VAR(id_file, id_var_obsy, netcdf_data_1d, & start = (/1/), count = (/ny-1+1/)) DEALLOCATE(netcdf_data_1d) !- z axis ALLOCATE(netcdf_data_1d(0:obs_nz)) DO k = 0, obs_nz netcdf_data_1d(k) = k ENDDO nc_stat = NF90_PUT_VAR(id_file, id_var_obsz, netcdf_data_1d, & start = (/1/), count = (/obs_nz+1/)) DEALLOCATE(netcdf_data_1d) #endif !- Close NetCDF file CASE ('close') #if defined ( __netcdf ) nc_stat = NF90_CLOSE(id_file) #endif !- Write data arrays to file CASE ('write') WRITE(*, '(A)') ' --- NetCDF Write! :-) ---' ! ALLOCATE(netcdf_data_2d(0:nx,0:ny)) ALLOCATE(netcdf_data_3d(0:nx,0:ny, 0:nsza)) !- Write geopotential DO k = 0, nsza DO i = 0, obs_nx DO j = 0, obs_ny netcdf_data_3d(i,j, k) = cube_rad(i,j, k) !normalweise j,i bei Feld2 ENDDO ENDDO ENDDO #if defined ( __netcdf ) nc_stat = NF90_PUT_VAR(id_file, id_var_rad, netcdf_data_3d, & start = (/1, 1, 1/), & count = (/obs_nx+1, obs_ny+1, nsza+1/)) nc_stat = NF90_PUT_VAR(id_file, id_var_int, integration, & start = (/1, 1, 1/), & count = (/obs_nx+1, obs_ny+1/)) nc_stat = NF90_PUT_VAR(id_file, id_var_pfs, pfs, & start = (/1, 1, 1/), & count = (/obs_nx+1, obs_ny+1, 3/)) nc_stat = NF90_PUT_VAR(id_file, id_var_irr, irrad, & start = (/1, 1, 1/), & count = (/3, nsza+1/)) DEALLOCATE(netcdf_data_3d) ALLOCATE(netcdf_data_3d(0:nx-1,0:ny-1, 0:obs_nz)) !- Write geopotential DO k = 0, obs_nz DO i = 0, nx-1 DO j = 0, ny-1 netcdf_data_3d(i,j, k) = obsi(i,j, k) !normalweise j,i bei Feld2 ENDDO ENDDO ENDDO nc_stat = NF90_PUT_VAR(id_file, id_var_obs, netcdf_data_3d, & start = (/1, 1, 1/), & count = (/nx-1+1, ny-1+1, obs_nz+1/)) DEALLOCATE(netcdf_data_3d) #endif !- Print error message if unknown action selected CASE DEFAULT WRITE(*, '(A)') ' ** data_output: action "'// & TRIM(action)//'"unknown!' END SELECT END SUBROUTINE data_output END PROGRAM obstruction