Changeset 4245 for palm/trunk/SOURCE/netcdf_data_input_mod.f90
- Timestamp:
- Sep 30, 2019 8:40:37 AM (5 years ago)
- Location:
- palm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE
- Property svn:mergeinfo changed
/palm/branches/resler/SOURCE merged: 4192-4193,4199-4200,4202,4206-4207,4240,4244
- Property svn:mergeinfo changed
-
palm/trunk/SOURCE/netcdf_data_input_mod.f90
r4226 r4245 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Add reading and processing of building_surface_pars 28 ! 29 ! 4226 2019-09-10 17:03:24Z suehring 27 30 ! - Netcdf input routine for dimension length renamed 28 31 ! - Move offline-nesting-specific checks to nesting_offl_mod … … 376 379 REAL(wp) :: fill1 = -9999.9_wp !< fill values for lod = 1 377 380 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_2d !< 2d variable (lod = 1) 381 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: oro_max !< terraing height under particular buildings 378 382 END TYPE build_in 379 383 … … 416 420 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pars_xyz !< respective parameters, level of detail = 2 417 421 END TYPE pars 422 ! 423 !-- Data type for surface parameter lists 424 TYPE pars_surf 425 INTEGER(iwp) :: np !< total number of parameters 426 INTEGER(iwp) :: nsurf !< number of local surfaces 427 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: index_ji !< index for beginning and end of surfaces at (j,i) 428 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: coords !< (k,j,i,norm_z,norm_y,norm_x) 429 !< k,j,i: surface position 430 !< norm_z,norm_y,norm_x: surface normal vector 431 432 LOGICAL :: from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used 433 434 REAL(wp) :: fill = -9999.9_wp !< fill value 435 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pars !< respective parameters per surface 436 END TYPE pars_surf 418 437 ! 419 438 !-- Define type for global file attributes … … 551 570 TYPE(pars) :: water_pars_f !< input variable for water parameters 552 571 572 TYPE(pars_surf) :: building_surface_pars_f !< input variable for building surface parameters 573 553 574 TYPE(chem_emis_att_type) :: chem_emis_att !< Input Information of Chemistry Emission Data from netcdf 554 575 TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:) :: chem_emis !< Input Chemistry Emission Data from netcdf … … 683 704 !-- Public variables 684 705 PUBLIC albedo_pars_f, albedo_type_f, basal_area_density_f, buildings_f, & 685 building_id_f, building_pars_f, building_type_f, & 706 building_id_f, building_pars_f, building_surface_pars_f, & 707 building_type_f, & 686 708 char_fill, & 687 709 char_lod, & … … 1958 1980 building_pars_f%from_file = .FALSE. 1959 1981 ENDIF 1982 ! 1983 !-- Read building surface parameters 1984 IF ( check_existence( var_names, 'building_surface_pars' ) ) THEN 1985 building_surface_pars_f%from_file = .TRUE. 1986 CALL get_attribute( id_surf, char_fill, & 1987 building_surface_pars_f%fill, & 1988 .FALSE., 'building_surface_pars' ) 1989 ! 1990 !-- Read building_surface_pars 1991 CALL get_variable_surf( id_surf, 'building_surface_pars', & 1992 building_surface_pars_f ) 1993 ELSE 1994 building_surface_pars_f%from_file = .FALSE. 1995 ENDIF 1960 1996 1961 1997 ! … … 4302 4338 ENDIF 4303 4339 ! 4304 !-- Extrapolate if actual height is above the highest Inifor level 4340 !-- Extrapolate if actual height is above the highest Inifor level by the last value 4305 4341 ELSE 4306 var(k) = var_file(ku) + ( var_file(ku) - var_file(ku-1) ) / & 4307 ( z_file(ku) - z_file(ku-1) ) * & 4308 ( z_grid(k) - z_file(ku) ) 4309 4342 var(k) = var_file(ku) 4310 4343 ENDIF 4311 4344 … … 4965 4998 #endif 4966 4999 END SUBROUTINE get_variable_pr 5000 5001 5002 !------------------------------------------------------------------------------! 5003 ! Description: 5004 ! ------------ 5005 !> Reads a per-surface pars variable from file. Because all surfaces are stored 5006 !> as flat 1-D array, each PE has to scan the data and find the surface indices 5007 !> belonging to its subdomain. During this scan, it also builds a necessary 5008 !> (j,i) index. 5009 !------------------------------------------------------------------------------! 5010 SUBROUTINE get_variable_surf( id, variable_name, surf ) 5011 #if defined( __netcdf ) 5012 5013 USE pegrid 5014 5015 USE indices, & 5016 ONLY: nxl, nxr, nys, nyn, nzb, nzt 5017 5018 USE control_parameters, & 5019 ONLY: dz, message_string 5020 5021 USE grid_variables, & 5022 ONLY: dx, dy 5023 5024 USE basic_constants_and_equations_mod, & 5025 ONLY: pi 5026 5027 USE arrays_3d, & 5028 ONLY: zu, zw 5029 5030 IMPLICIT NONE 5031 5032 INTEGER, PARAMETER :: nsurf_pars_read = 1024**2 !< read buffer size 5033 5034 CHARACTER(LEN=*) :: variable_name !< variable name 5035 5036 INTEGER(iwp), DIMENSION(6) :: coords !< integer coordinates of surface 5037 INTEGER(iwp) :: i, j 5038 INTEGER(iwp) :: isurf !< netcdf surface index 5039 INTEGER(iwp) :: is !< local surface index 5040 INTEGER(iwp), INTENT(IN) :: id !< file id 5041 INTEGER(iwp), DIMENSION(2) :: id_dim !< dimension ids 5042 INTEGER(iwp) :: id_var !< variable id 5043 INTEGER(iwp) :: id_zs !< zs variable id 5044 INTEGER(iwp) :: id_ys !< ys variable id 5045 INTEGER(iwp) :: id_xs !< xs variable id 5046 INTEGER(iwp) :: id_zenith !< zeith variable id 5047 INTEGER(iwp) :: id_azimuth !< azimuth variable id 5048 INTEGER(iwp) :: is0, isc !< read surface start and count 5049 INTEGER(iwp) :: nsurf !< total number of surfaces in file 5050 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nsurf_ji !< numbers of surfaces by coords 5051 5052 TYPE(pars_surf) :: surf !< parameters variable to be loaded 5053 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pars_read !< read buffer 5054 REAL(wp), DIMENSION(:), ALLOCATABLE :: zs, ys, xs !< read buffer for zs(s), ys, xs 5055 REAL(wp), DIMENSION(:), ALLOCATABLE :: zenith !< read buffer for zenith(s) 5056 REAL(wp), DIMENSION(:), ALLOCATABLE :: azimuth !< read buffer for azimuth(s) 5057 REAL(wp) :: oro_max_l !< maximum terrain height under building 5058 5059 ! 5060 !-- First, inquire variable ID 5061 nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var ) 5062 nc_stat = NF90_INQ_VARID( id, 'zs', id_zs ) 5063 nc_stat = NF90_INQ_VARID( id, 'ys', id_ys ) 5064 nc_stat = NF90_INQ_VARID( id, 'xs', id_xs ) 5065 nc_stat = NF90_INQ_VARID( id, 'zenith', id_zenith ) 5066 nc_stat = NF90_INQ_VARID( id, 'azimuth', id_azimuth ) 5067 ! 5068 !-- Inquire dimension sizes 5069 nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim ) 5070 nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(1), LEN = nsurf ) 5071 nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(2), LEN = surf%np ) 5072 5073 ALLOCATE ( pars_read( nsurf_pars_read, surf%np ), & 5074 zs(nsurf_pars_read), ys(nsurf_pars_read), & 5075 xs(nsurf_pars_read), zenith(nsurf_pars_read), & 5076 azimuth(nsurf_pars_read), & 5077 nsurf_ji(nys:nyn, nxl:nxr) ) 5078 5079 nsurf_ji(:,:) = 0 5080 ! 5081 !-- Scan surface coordinates, count local 5082 is0 = 1 5083 DO 5084 isc = MIN(nsurf_pars_read, nsurf - is0 + 1) 5085 IF ( isc <= 0 ) EXIT 5086 5087 nc_stat = NF90_GET_VAR( id, id_ys, ys, & 5088 start = (/ is0 /), & 5089 count = (/ isc /) ) 5090 nc_stat = NF90_GET_VAR( id, id_xs, xs, & 5091 start = (/ is0 /), & 5092 count = (/ isc /) ) 5093 nc_stat = NF90_GET_VAR( id, id_zenith, zenith, & 5094 start = (/ is0 /), & 5095 count = (/ isc /) ) 5096 nc_stat = NF90_GET_VAR( id, id_azimuth, azimuth, & 5097 start = (/ is0 /), & 5098 count = (/ isc /) ) 5099 CALL handle_error( 'get_variable_surf', 682, 'azimuth' ) 5100 5101 DO isurf = 1, isc 5102 ! 5103 !-- Parse coordinates, detect if belongs to subdomain 5104 coords = transform_coords( xs(isurf), ys(isurf), & 5105 zenith(isurf), azimuth(isurf) ) 5106 IF ( coords(2) < nys .OR. coords(2) > nyn .OR. & 5107 coords(3) < nxl .OR. coords(3) > nxr ) CYCLE 5108 5109 nsurf_ji(coords(2), coords(3)) = nsurf_ji(coords(2), coords(3)) + 1 5110 ENDDO 5111 5112 is0 = is0 + isc 5113 ENDDO 5114 ! 5115 !-- Populate reverse index from surface counts 5116 ALLOCATE ( surf%index_ji( 2, nys:nyn, nxl:nxr ) ) 5117 5118 isurf = 1 5119 DO j = nys, nyn 5120 DO i = nxl, nxr 5121 surf%index_ji(:,j,i) = (/ isurf, isurf + nsurf_ji(j,i) - 1 /) 5122 isurf = isurf + nsurf_ji(j,i) 5123 ENDDO 5124 ENDDO 5125 5126 surf%nsurf = isurf - 1 5127 ALLOCATE( surf%pars( 0:surf%np-1, surf%nsurf ), & 5128 surf%coords( 6, surf%nsurf ) ) 5129 ! 5130 !-- Scan surfaces again, saving pars into allocated structures 5131 nsurf_ji(:,:) = 0 5132 is0 = 1 5133 DO 5134 isc = MIN(nsurf_pars_read, nsurf - is0 + 1) 5135 IF ( isc <= 0 ) EXIT 5136 5137 nc_stat = NF90_GET_VAR( id, id_var, pars_read(1:isc, 1:surf%np), & 5138 start = (/ is0, 1 /), & 5139 count = (/ isc, surf%np /) ) 5140 CALL handle_error( 'get_variable_surf', 683, variable_name ) 5141 5142 nc_stat = NF90_GET_VAR( id, id_zs, zs, & 5143 start = (/ is0 /), & 5144 count = (/ isc /) ) 5145 nc_stat = NF90_GET_VAR( id, id_ys, ys, & 5146 start = (/ is0 /), & 5147 count = (/ isc /) ) 5148 nc_stat = NF90_GET_VAR( id, id_xs, xs, & 5149 start = (/ is0 /), & 5150 count = (/ isc /) ) 5151 nc_stat = NF90_GET_VAR( id, id_zenith, zenith, & 5152 start = (/ is0 /), & 5153 count = (/ isc /) ) 5154 nc_stat = NF90_GET_VAR( id, id_azimuth, azimuth, & 5155 start = (/ is0 /), & 5156 count = (/ isc /) ) 5157 5158 DO isurf = 1, isc 5159 ! 5160 !-- Parse coordinates, detect if belongs to subdomain 5161 coords = transform_coords( xs(isurf), ys(isurf), & 5162 zenith(isurf), azimuth(isurf) ) 5163 IF ( coords(2) < nys .OR. coords(2) > nyn .OR. & 5164 coords(3) < nxl .OR. coords(3) > nxr ) CYCLE 5165 ! 5166 !-- Determine maximum terrain under building (base z-coordinate). Using 5167 !-- normal vector to locate building inner coordinates. 5168 oro_max_l = buildings_f%oro_max(coords(2)-coords(5), coords(3)-coords(6)) 5169 IF ( oro_max_l == buildings_f%fill1 ) THEN 5170 WRITE( message_string, * ) 'Found building surface on ' // & 5171 'non-building coordinates (xs, ys, zenith, azimuth): ', & 5172 xs(isurf), ys(isurf), zenith(isurf), azimuth(isurf) 5173 CALL message( 'get_variable_surf', 'PA0684', 2, 2, 0, 6, 0 ) 5174 ENDIF 5175 ! 5176 !-- Urban layer has no stretching, therefore using dz(1) instead of linear 5177 !-- searching through zu/zw 5178 coords(1) = NINT((zs(isurf) + oro_max_l) / dz(1) + & 5179 0.5_wp + 0.5_wp * coords(4), KIND=iwp) 5180 ! 5181 !-- Save surface entry 5182 is = surf%index_ji(1, coords(2), coords(3)) + nsurf_ji(coords(2), coords(3)) 5183 surf%pars(:,is) = pars_read(isurf,:) 5184 surf%coords(:,is) = coords(:) 5185 5186 nsurf_ji(coords(2), coords(3)) = nsurf_ji(coords(2), coords(3)) + 1 5187 ENDDO 5188 5189 is0 = is0 + isc 5190 ENDDO 5191 5192 DEALLOCATE( pars_read, zs, ys, xs, zenith, azimuth, nsurf_ji ) 5193 5194 CONTAINS 5195 5196 PURE FUNCTION transform_coords( x, y, zenith, azimuth ) 5197 5198 REAL(wp), INTENT(in) :: x, y !< surface centre coordinates in metres from origin 5199 REAL(wp), INTENT(in) :: zenith !< surface normal zenith angle in degrees 5200 REAL(wp), INTENT(in) :: azimuth !< surface normal azimuth angle in degrees 5201 5202 INTEGER(iwp), DIMENSION(6) :: transform_coords !< (k,j,i,norm_z,norm_y,norm_x) 5203 5204 transform_coords(4) = NINT(COS(zenith*pi/180._wp), KIND=iwp) 5205 IF ( transform_coords(4) == 0 ) THEN 5206 transform_coords(5) = NINT(COS(azimuth*pi/180._wp), KIND=iwp) 5207 transform_coords(6) = NINT(SIN(azimuth*pi/180._wp), KIND=iwp) 5208 ELSE 5209 transform_coords(5) = 0._wp 5210 transform_coords(6) = 0._wp 5211 ENDIF 5212 5213 transform_coords(1) = -999._wp ! not calculated here 5214 transform_coords(2) = NINT(y/dy - 0.5_wp + 0.5_wp*transform_coords(5), KIND=iwp) 5215 transform_coords(3) = NINT(x/dx - 0.5_wp + 0.5_wp*transform_coords(6), KIND=iwp) 5216 5217 END FUNCTION transform_coords 5218 5219 #endif 5220 END SUBROUTINE get_variable_surf 4967 5221 4968 5222
Note: See TracChangeset
for help on using the changeset viewer.
