Home

Context Navigation

calc_obstruction.f90

Last change on this file was 4870, checked in by suehring, 4 years ago
preprocessor to calculate UV-obstruction added
Property svn:keywords set to `Id`
File size: 21.1 KB

Rev	Line
[4870]	1	!> @file calc_obstruction.f90
	2	!------------------------------------------------------------------------------!
	3	! This file is part of the PALM model system.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 1997-2021 Leibniz Universitaet Hannover
	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! -----------------
	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: calc_obstruction.f90 4870 2021-02-08 18:09:41Z banzhafs $
	27	! Initial revision
	28	!
	29	!
	30	! Authors:
	31	! --------
	32	! @author Michael Schrempf (main author),
	33	! Matthias SÃŒhring (first major revision of the code)
	34	!
	35	!> @todo - Code formatting
	36	!------------------------------------------------------------------------------!
	37	! Description:
	38	! ------------
	39	!> Calculation of sky-obstruction.
	40	!------------------------------------------------------------------------------!
	41
	42	PROGRAM obstruction
	43	#if defined ( __netcdf )
	44	USE netcdf
	45	#endif
	46
	47	IMPLICIT NONE
	48
	49	CHARACTER(LEN=10) :: b(3)
	50	CHARACTER(LEN=300) :: input_file
	51	CHARACTER(LEN=300) :: path_to_expo_files
	52
	53	INTEGER:: spx, spy, counter, cc, sza, rza, azi, str_i, str_j, ii, &
	54	vector_length_x, vector_length_y, id, id_var, nc_stat, i, j, k
	55	INTEGER(KIND=1) :: bitvar
	56	INTEGER, DIMENSION(1:8):: date_time
	57	INTEGER, DIMENSION(0:99):: valy, valx
	58
	59	INTEGER :: obs_ny = 9 !< number of zenith values
	60	INTEGER :: obs_nx = 35 !< number of azimuth values
	61	INTEGER :: obs_nz = 44 !< ?
	62
	63	INTEGER(KIND=2), DIMENSION(:,:,:), ALLOCATABLE :: obs
	64	INTEGER(KIND=2), DIMENSION(:,:,:), ALLOCATABLE :: obsi
	65
	66	! INTEGER(KIND=2), DIMENSION(0:1159,0:1159,0:44):: obs
	67	! INTEGER(KIND=2), DIMENSION(0:959,0:959,0:44):: obsi
	68
	69	INTEGER :: nx !< number of grid points in x-direction
	70	INTEGER :: ny !< number of grid points in y-direction
	71	INTEGER :: nsza
	72	REAL :: dx !< grid spacing in x
	73	REAL :: dy !< grid spacing in y
	74	REAL :: dz !< grid spacing in z
	75
	76	REAL :: delta = 10.0, fill !< delta of zenith and azimuth angles (namelist param)
	77
	78	REAL:: dtor, pi = 3.141592654, t1 = 0.0, t2 = 0.0
	79	REAL, DIMENSION(0:99):: value_y, value_x, value_z, value_zn, distance_x, distance_y
	80	! REAL, DIMENSION(0:1159,0:1159):: topo
	81	! REAL, DIMENSION(0:959,0:959):: topo_0
	82	REAL, DIMENSION(:,:), ALLOCATABLE :: topo
	83	REAL, DIMENSION(:,:), ALLOCATABLE :: topo_0
	84
	85	REAL(kind=8), DIMENSION(0:35,0:9) :: integration
	86	REAL, DIMENSION(0:2,0:90) :: irrad
	87	REAL, DIMENSION(0:35,0:9) :: tmp
	88	REAL, DIMENSION(0:35,0:9,0:90) :: cube_rad
	89	REAL, DIMENSION(0:35,0:9,0:2) :: pfs
	90
	91	CHARACTER(LEN=15) :: output ='Data_input_all'
	92
	93	CALL obstruct_parin
	94
	95	nsza = 90
	96
	97	dtor = pi / 180.
	98	vector_length_x = INT( 99.0 / dx )
	99	vector_length_y = INT( 99.0 / dy )
	100
	101
	102	ALLOCATE( obs(0:nx-1+(2vector_length_x),0:ny-1+(2vector_length_y),0:44) )
	103	ALLOCATE( obsi(0:nx-1,0:ny-1,0:44) )
	104	ALLOCATE( topo(0:nx-1+(2vector_length_x),0:ny-1+(2vector_length_y)) )
	105	ALLOCATE( topo_0(0:nx-1,0:ny-1) )
	106
	107	obs(:,:,:) = 9
	108	obsi(:,:,:) = 0
	109	topo(:,:) = 0
	110
	111	call date_and_time(b(1), b(2), b(3), date_time)
	112	t1 = date_time(5)+ date_time(6)/60. + (date_time(7)+(date_time(8)/1000.))/3600.
	113
	114	distance_x = 0.0
	115	distance_y = 0.0
	116
	117	DO ii = 0, 99
	118	distance_x(ii) = ii * dx
	119	distance_y(ii) = ii * dy
	120	ENDDO
	121	distance_x = distance_x + dx
	122	distance_y = distance_y + dy
	123
	124	#if defined ( __netcdf )
	125	nc_stat = NF90_OPEN( input_file, NF90_NOWRITE, id )
	126	nc_stat = NF90_INQ_VARID( id, 'buildings_2d', id_var )
	127	nc_stat = NF90_GET_VAR( id, id_var, topo_0, start = (/ 1, 1 /), count = (/ nx, ny /) )
	128	nc_stat = NF90_GET_ATT( id, id_var, '_FillValue', fill )
	129	nc_stat = NF90_CLOSE( id )
	130	#endif
	131
	132	! OPEN(11, FILE='/home/matthias/UNI/MOSAIK/UV/test/Topo_input/topo_Reuter.txt', FORM='FORMATTED')
	133	! DO rza = 0, nx-1
	134	! READ(11,'(960(F5.1,1X))') topo_0(:,rza) !(nx-1)-rza
	135	! ENDDO
	136	! CLOSE(11)
	137
	138	topo(vector_length_x:nx-1+vector_length_x,vector_length_y:ny-1+vector_length_y) = &
	139	MERGE( topo_0, 0.0, topo_0 /= fill )
	140
	141
	142	Do spy = vector_length_y, ny-1+vector_length_y
	143	DO spx = vector_length_x, nx-1+vector_length_x
	144	print*, spy, spx, topo(spx,spy)
	145	IF (topo(spx,spy) == 0.0 ) THEN
	146	! print*, "in if"
	147	obs(spx,spy,:)=0
	148	cc = 0
	149	bitvar=2
	150	DO sza = 0,90, 10
	151	! print*, sza
	152	DO azi = 0,359, 10
	153	! print*, azi
	154	value_zn = SIN((sza+0.0000001)*dtor)
	155
	156	value_x = ANINT((distance_x(:)SIN(azidtor)) *value_zn(:))
	157	value_y = ANINT((distance_y(:)COS(azidtor)) *value_zn(:))
	158	! print*, value_x
	159	! print*, value_y
	160	valx = INT(value_x)
	161	valy = INT(value_y)
	162	! print*, "valx", valx
	163	!value_z = (distance(:)COS(szadtor))
	164	value_z = COS(szadtor) SQRT( value_x2 + value_y2) / value_zn
	165	! print, szadtor, sza
	166	!print*, value_z
	167	! stop
	168	counter=0
	169	DO str_i = 0, vector_length_x-1
	170	DO str_j = 0, vector_length_y-1
	171	IF ( topo( spx+valx(str_i), spy+valy(str_j) ) >= (value_z(str_i)*dz) ) THEN
	172	IF ( topo( spx+valx(str_i), spy+valy(str_j) ) == 0) THEN
	173	counter=counter + 0
	174	ELSE
	175	counter=counter + 1
	176	ENDIF
	177	ENDIF
	178	ENDDO
	179	ENDDO
	180	IF (sza == 90 ) Counter = 1
	181	IF (counter == 0) bitvar = IBSET(bitvar,MOD(cc,8))
	182	IF (counter > 0) bitvar = IBCLR(bitvar,MOD(cc,8))
	183	!IF (counter .GT. 0) bitvar = IBSET(bitvar,MOD(cc,8))
	184
	185
	186	IF ( MOD(cc,8) == 7 ) THEN
	187	obs(spx,spy,cc/8) = bitvar
	188	bitvar = 2
	189	ENDIF
	190	cc = cc+1
	191	ENDDO
	192	ENDDO
	193	ENDIF
	194	ENDDO
	195	ENDDO
	196
	197	call date_and_time(b(1), b(2), b(3), date_time)
	198	t2 = date_time(5)+ date_time(6)/60. + (date_time(7)+(date_time(8)/1000.))/3600.
	199
	200	obsi(:,:,:)=obs(vector_length_x:nx-1+vector_length_x,vector_length_y:ny-1+vector_length_y,:)
	201
	202
	203
	204	OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Integration_solid_angles.txt', FORM='FORMATTED')
	205	DO rza = 0,9
	206	READ(11,'(36(F10.8,2X))') integration(:,rza)
	207	ENDDO
	208	CLOSE(11)
	209
	210	OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Irradiance_ROM.txt', FORM='FORMATTED')
	211	Do sza = 0,90
	212	READ(11,'(3(F10.6,2X))') irrad(:,sza)
	213	ENDDO
	214	CLOSE(11)
	215
	216
	217
	218	OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Projection_areas_CL0.txt', FORM='FORMATTED')
	219	DO rza = 0,9
	220	READ(11,'(36(F10.6,2X))') tmp(:,rza)
	221	ENDDO
	222	CLOSE(11)
	223	pfs(:,:,0) = tmp(:,:)
	224
	225	OPEN(11, FILE=TRIM( path_to_expo_files ) // 'Projection_areas_CL1.txt', FORM='FORMATTED')
	226	DO rza = 0,9
	227	READ(11,'(36(F10.6,2X))') tmp(:,rza)
	228	ENDDO
	229	CLOSE(11)
	230	pfs(:,:,1) = tmp(:,:)
	231
	232	OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Projection_areas_CL3.txt', FORM='FORMATTED')
	233	DO rza = 0,9
	234	READ(11,'(36(F10.6,2X))') tmp(:,rza)
	235	ENDDO
	236	CLOSE(11)
	237	pfs(:,:,2) = tmp(:,:)
	238
	239
	240	OPEN(11, FILE= TRIM( path_to_expo_files ) // 'Radiance_ROM.txt', FORM='FORMATTED')
	241	Do sza = 0,90
	242	DO rza = 0,9
	243	READ(11,'(36(F10.6,2X))') tmp(:,rza)
	244	ENDDO
	245	cube_rad(:,:,sza) = tmp(:,:)
	246	ENDDO
	247	CLOSE(11)
	248
	249	!- Finalize output
	250	CALL data_output('open', output)
	251	CALL data_output('write', output)
	252	CALL data_output('close', output)
	253
	254
	255	!- End
	256	WRITE(*, '(A)') ' --- creation finished! :-) ---'
	257
	258
	259	! open(unit=12, status='replace',file=&
	260	! '/home/matthias/UNI/MOSAIK/UV/test/result_shading/obstruction_Berlin_reuter.bin', FORM='UNFORMATTED')
	261	! write(12) obsi
	262	! close(12)
	263
	264	PRINT *,'Time_Start = ',t1
	265	PRINT *,'Time_ENDE = ',t2
	266	PRINT ,'elapsed time used (dateT) = ',(t2-t1)3600.
	267
	268
	269	CONTAINS
	270
	271	SUBROUTINE obstruct_parin
	272
	273	IMPLICIT NONE
	274
	275	INTEGER :: file_id_parin = 90
	276
	277	NAMELIST /uv_parin/ dx, dy, dz, nx, ny, input_file, path_to_expo_files
	278	!
	279	!-- Open namelist file.
	280	OPEN( file_id_parin, FILE='uv_parin', STATUS='OLD', FORM='FORMATTED')
	281	!
	282	!-- Read namelist.
	283	READ ( file_id_parin, uv_parin )
	284	!
	285	!-- Close namelist file.
	286	CLOSE( file_id_parin )
	287
	288	END SUBROUTINE obstruct_parin
	289
	290
	291	SUBROUTINE data_output(action,output)
	292	! Data output to NetCDF file
	293
	294	CHARACTER(LEN=*), INTENT(IN) :: action !< flag to steer routine (open/close/write)
	295	CHARACTER(LEN=*), INTENT(IN) :: output !< output file name prefix
	296	! CHARACTER(LEN=10) :: output ='Data_input'
	297
	298	INTEGER, SAVE :: do_count=0 !< counting output
	299	INTEGER, SAVE :: id_dim_time !< ID of dimension time
	300	INTEGER, SAVE :: id_dim_x !< ID of dimension x
	301	INTEGER, SAVE :: id_dim_xirr !< ID of dimension x
	302	INTEGER, SAVE :: id_dim_y !< ID of dimension y
	303	INTEGER, SAVE :: id_dim_z !< ID of dimension z
	304	INTEGER, SAVE :: id_dim_obsx !< ID of dimension x
	305	INTEGER, SAVE :: id_dim_obsy !< ID of dimension y
	306	INTEGER, SAVE :: id_dim_obsz !< ID of dimension z
	307	INTEGER, SAVE :: id_file !< ID of NetCDF file
	308	INTEGER, SAVE :: id_var_int !< ID of geopotential
	309	INTEGER, SAVE :: id_var_irr !< ID of geopotential
	310	INTEGER, SAVE :: id_var_obs !< ID of geopotential
	311	INTEGER, SAVE :: id_var_obsx !< ID of x
	312	INTEGER, SAVE :: id_var_obsy !< ID of y
	313	INTEGER, SAVE :: id_var_obsz !< ID of y
	314
	315	INTEGER, SAVE :: id_var_pfs !< ID of geopotential
	316	INTEGER, SAVE :: id_var_rad !< ID of geopotential
	317	INTEGER, SAVE :: id_var_x !< ID of x
	318	INTEGER, SAVE :: id_var_xirr !< ID of x
	319	INTEGER, SAVE :: id_var_y !< ID of y
	320	INTEGER, SAVE :: id_var_z !< ID of y
	321	INTEGER, SAVE :: nc_stat !< status flag of NetCDF routines
	322
	323	REAL, DIMENSION(:), ALLOCATABLE :: netcdf_data_1d !< 1D output array
	324
	325	REAL, DIMENSION(:,:), ALLOCATABLE :: netcdf_data_2d !< 2D output array
	326
	327	REAL, DIMENSION(:,:,:), ALLOCATABLE :: netcdf_data_3d !< 3D output array
	328
	329
	330	SELECT CASE (TRIM(action))
	331
	332	!- Initialize output file
	333	CASE ('open')
	334
	335	!- Delete any pre-existing output file
	336	OPEN(20, FILE=TRIM(output)//'.nc')
	337	CLOSE(20, STATUS='DELETE')
	338
	339	!- Open file
	340	#if defined ( __netcdf )
	341	nc_stat = NF90_CREATE(TRIM(output)//'.nc', NF90_NOCLOBBER, id_file)
	342	IF (nc_stat /= NF90_NOERR) PRINT*, '+++ netcdf error'
	343
	344	!- Write global attributes
	345	nc_stat = NF90_PUT_ATT(id_file, NF90_GLOBAL, &
	346	'Conventions', 'COARDS')
	347	nc_stat = NF90_PUT_ATT(id_file, NF90_GLOBAL, &
	348	'title', 'UVEM-model')
	349
	350	!- Define time coordinate
	351	nc_stat = NF90_DEF_DIM(id_file, 'TIME', NF90_UNLIMITED, &
	352	id_dim_time)
	353
	354
	355	!- Define spatial coordinates
	356	nc_stat = NF90_DEF_DIM(id_file, 'AZI', obs_nx+1, id_dim_x)
	357	nc_stat = NF90_DEF_VAR(id_file, 'AZI', NF90_DOUBLE, id_dim_x, id_var_x)
	358	nc_stat = NF90_PUT_ATT(id_file, id_var_x, 'units', 'Grad')
	359
	360	nc_stat = NF90_DEF_DIM(id_file, 'ZEN', obs_ny+1, id_dim_y)
	361	nc_stat = NF90_DEF_VAR(id_file, 'ZEN', NF90_DOUBLE, &
	362	id_dim_y, id_var_y)
	363	nc_stat = NF90_PUT_ATT(id_file, id_var_y, 'units', 'Grad')
	364
	365	nc_stat = NF90_DEF_DIM(id_file, 'sza', nsza+1, id_dim_z)
	366	nc_stat = NF90_DEF_VAR(id_file, 'sza', NF90_DOUBLE, &
	367	id_dim_z, id_var_z)
	368	nc_stat = NF90_PUT_ATT(id_file, id_var_z, 'units', 'Grad')
	369
	370	!- Define spatial coordinates
	371	nc_stat = NF90_DEF_DIM(id_file, 'COM', 3, id_dim_xirr)
	372	nc_stat = NF90_DEF_VAR(id_file, 'COM', NF90_DOUBLE, id_dim_xirr, id_var_xirr)
	373	nc_stat = NF90_PUT_ATT(id_file, id_var_xirr, 'units', 'Komponente')
	374
	375
	376	!- Define spatial coordinates for obstruction file
	377	nc_stat = NF90_DEF_DIM(id_file, 'X', nx-1+1, id_dim_obsx)
	378	nc_stat = NF90_DEF_VAR(id_file, 'X', NF90_DOUBLE, id_dim_obsx, id_var_obsx)
	379	nc_stat = NF90_PUT_ATT(id_file, id_var_obsx, 'units', 'meters')
	380
	381	nc_stat = NF90_DEF_DIM(id_file, 'Y', ny-1+1, id_dim_obsy)
	382	nc_stat = NF90_DEF_VAR(id_file, 'Y', NF90_DOUBLE, id_dim_obsy, id_var_obsy)
	383	nc_stat = NF90_PUT_ATT(id_file, id_var_obsy, 'units', 'meters')
	384
	385	nc_stat = NF90_DEF_DIM(id_file, 'directions', nsza+1, id_dim_obsz)
	386	nc_stat = NF90_DEF_VAR(id_file, 'directions', NF90_DOUBLE, id_dim_obsz, id_var_obsz)
	387	nc_stat = NF90_PUT_ATT(id_file, id_var_obsz, 'units', 'ibits')
	388
	389
	390
	391	!- Define output arrays
	392	nc_stat = NF90_DEF_VAR(id_file, 'radiance', NF90_DOUBLE, &
	393	(/id_dim_x, id_dim_y, id_dim_z/), id_var_rad)
	394	nc_stat = NF90_PUT_ATT(id_file, id_var_rad, &
	395	'long_name', 'vitD weighetd radiance at 300DU')
	396	nc_stat = NF90_PUT_ATT(id_file, id_var_rad, &
	397	'short_name', 'vitD radiance')
	398	nc_stat = NF90_PUT_ATT(id_file, id_var_rad, 'units', 'W/m2 sr')
	399
	400
	401
	402	nc_stat = NF90_DEF_VAR(id_file, 'int_factors', NF90_DOUBLE, &
	403	(/id_dim_x, id_dim_y/), id_var_int)
	404	nc_stat = NF90_PUT_ATT(id_file, id_var_int, &
	405	'long_name', 'integration factors')
	406	nc_stat = NF90_PUT_ATT(id_file, id_var_int, &
	407	'short_name', 'integration factors')
	408	! nc_stat = NF90_PUT_ATT(id_file, id_var_int, 'units', '-')
	409
	410
	411	!- Define output arrays
	412	nc_stat = NF90_DEF_VAR(id_file, 'projarea', NF90_DOUBLE, &
	413	(/id_dim_x, id_dim_y, id_dim_time/), id_var_pfs)
	414	nc_stat = NF90_PUT_ATT(id_file, id_var_pfs, &
	415	'long_name', 'projection areas of a human')
	416	nc_stat = NF90_PUT_ATT(id_file, id_var_pfs, &
	417	'short_name', 'projection areas')
	418	nc_stat = NF90_PUT_ATT(id_file, id_var_pfs, 'units', 'm2')
	419
	420
	421	!- Define output arrays
	422	nc_stat = NF90_DEF_VAR(id_file, 'irradiance', NF90_DOUBLE, &
	423	(/id_dim_xirr, id_dim_z/), id_var_irr)
	424	nc_stat = NF90_PUT_ATT(id_file, id_var_irr, &
	425	'long_name', 'global, diffuse and direct solar irradiance')
	426	nc_stat = NF90_PUT_ATT(id_file, id_var_irr, &
	427	'short_name', 'irradiance values')
	428	nc_stat = NF90_PUT_ATT(id_file, id_var_irr, 'units', 'W/m2')
	429
	430
	431	!- Define output arrays for obstruction file
	432	nc_stat = NF90_DEF_VAR(id_file, 'obstruction', NF90_DOUBLE, &
	433	(/id_dim_obsx, id_dim_obsy, id_dim_obsz/), id_var_obs)
	434	nc_stat = NF90_PUT_ATT(id_file, id_var_obs, &
	435	'long_name', 'calculated obstructions in different directions')
	436	nc_stat = NF90_PUT_ATT(id_file, id_var_obs, &
	437	'short_name', 'obstructions')
	438	nc_stat = NF90_PUT_ATT(id_file, id_var_obs, 'units', '-')
	439
	440
	441
	442
	443	!- Leave define mode
	444	nc_stat = NF90_ENDDEF(id_file)
	445
	446	!- Write axis to file
	447	ALLOCATE(netcdf_data_1d(0:nx))
	448
	449	!- x axis
	450	DO i = 0, nx
	451	netcdf_data_1d(i) = i * delta
	452	ENDDO
	453	nc_stat = NF90_PUT_VAR(id_file, id_var_x, netcdf_data_1d, &
	454	start = (/1/), count = (/nx+1/))
	455	DO i = 0, 2
	456	netcdf_data_1d(i) = i
	457	ENDDO
	458	nc_stat = NF90_PUT_VAR(id_file, id_var_xirr, netcdf_data_1d, &
	459	start = (/1/), count = (/3/))
	460
	461	!- y axis
	462	DO j = 0, ny
	463	netcdf_data_1d(j) = j * delta
	464	ENDDO
	465	nc_stat = NF90_PUT_VAR(id_file, id_var_y, netcdf_data_1d, &
	466	start = (/1/), count = (/ny+1/))
	467	DEALLOCATE(netcdf_data_1d)
	468
	469	!- z axis
	470	ALLOCATE(netcdf_data_1d(0:nsza))
	471	DO k = 0, nsza
	472	netcdf_data_1d(k) = k
	473	ENDDO
	474	nc_stat = NF90_PUT_VAR(id_file, id_var_z, netcdf_data_1d, &
	475	start = (/1/), count = (/nsza+1/))
	476	DEALLOCATE(netcdf_data_1d)
	477
	478
	479
	480	!- define axis for obstruction file
	481	!- x axis
	482	ALLOCATE(netcdf_data_1d(0:nx-1))
	483	DO i = 0, nx-1
	484	netcdf_data_1d(i) = i * dx
	485	ENDDO
	486	nc_stat = NF90_PUT_VAR(id_file, id_var_obsx, netcdf_data_1d, &
	487	start = (/1/), count = (/nx-1+1/))
	488	!- y axis
	489	DO j = 0, ny-1
	490	netcdf_data_1d(j) = j * dy
	491	ENDDO
	492	nc_stat = NF90_PUT_VAR(id_file, id_var_obsy, netcdf_data_1d, &
	493	start = (/1/), count = (/ny-1+1/))
	494	DEALLOCATE(netcdf_data_1d)
	495
	496	!- z axis
	497	ALLOCATE(netcdf_data_1d(0:obs_nz))
	498	DO k = 0, obs_nz
	499	netcdf_data_1d(k) = k
	500	ENDDO
	501	nc_stat = NF90_PUT_VAR(id_file, id_var_obsz, netcdf_data_1d, &
	502	start = (/1/), count = (/obs_nz+1/))
	503	DEALLOCATE(netcdf_data_1d)
	504	#endif
	505
	506
	507	!- Close NetCDF file
	508	CASE ('close')
	509	#if defined ( __netcdf )
	510	nc_stat = NF90_CLOSE(id_file)
	511	#endif
	512
	513	!- Write data arrays to file
	514	CASE ('write')
	515	WRITE(*, '(A)') ' --- NetCDF Write! :-) ---'
	516
	517	! ALLOCATE(netcdf_data_2d(0:nx,0:ny))
	518	ALLOCATE(netcdf_data_3d(0:nx,0:ny, 0:nsza))
	519
	520	!- Write geopotential
	521	DO k = 0, nsza
	522	DO i = 0, obs_nx
	523	DO j = 0, obs_ny
	524	netcdf_data_3d(i,j, k) = cube_rad(i,j, k) !normalweise j,i bei Feld2
	525	ENDDO
	526	ENDDO
	527	ENDDO
	528	#if defined ( __netcdf )
	529	nc_stat = NF90_PUT_VAR(id_file, id_var_rad, netcdf_data_3d, &
	530	start = (/1, 1, 1/), &
	531	count = (/obs_nx+1, obs_ny+1, nsza+1/))
	532
	533	nc_stat = NF90_PUT_VAR(id_file, id_var_int, integration, &
	534	start = (/1, 1, 1/), &
	535	count = (/obs_nx+1, obs_ny+1/))
	536
	537	nc_stat = NF90_PUT_VAR(id_file, id_var_pfs, pfs, &
	538	start = (/1, 1, 1/), &
	539	count = (/obs_nx+1, obs_ny+1, 3/))
	540
	541	nc_stat = NF90_PUT_VAR(id_file, id_var_irr, irrad, &
	542	start = (/1, 1, 1/), &
	543	count = (/3, nsza+1/))
	544
	545	DEALLOCATE(netcdf_data_3d)
	546
	547
	548	ALLOCATE(netcdf_data_3d(0:nx-1,0:ny-1, 0:obs_nz))
	549
	550	!- Write geopotential
	551	DO k = 0, obs_nz
	552	DO i = 0, nx-1
	553	DO j = 0, ny-1
	554	netcdf_data_3d(i,j, k) = obsi(i,j, k) !normalweise j,i bei Feld2
	555	ENDDO
	556	ENDDO
	557	ENDDO
	558	nc_stat = NF90_PUT_VAR(id_file, id_var_obs, netcdf_data_3d, &
	559	start = (/1, 1, 1/), &
	560	count = (/nx-1+1, ny-1+1, obs_nz+1/))
	561
	562	DEALLOCATE(netcdf_data_3d)
	563	#endif
	564
	565
	566	!- Print error message if unknown action selected
	567	CASE DEFAULT
	568
	569	WRITE(, '(A)') ' * data_output: action "'// &
	570	TRIM(action)//'"unknown!'
	571
	572	END SELECT
	573
	574	END SUBROUTINE data_output
	575
	576	END PROGRAM obstruction

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/UTIL/uv_obstruction_preprocessor/calc_obstruction.f90

Download in other formats: