Home

Context Navigation

data_output_3d.f90 @ 1328

Last change on this file since 1328 was 1327, checked in by raasch, 11 years ago

Changed:

-s real64 removed (.mrun.config.hlrnIII)
-r8 removed (.mrun.config.imuk)
deleted: .mrun.config.imuk_ice2_netcdf4 .mrun.config.imuk_hlrn

REAL constants defined as wp-kind in modules

"baroclinicity" renamed "baroclinity", "ocean version" replaced by
"ocean mode"

code parts concerning old output formats "iso2d" and "avs" removed.
netCDF is the only remaining output format.

Errors:

bugfix: duplicate error message 56 removed

Property svn:keywords set to Id

File size: 18.8 KB

Rev	Line
[1]	1	SUBROUTINE data_output_3d( av )
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later 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	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	18	!--------------------------------------------------------------------------------!
	19	!
[254]	20	! Current revisions:
[1106]	21	! ------------------
[1327]	22	! parts concerning avs output removed,
	23	! -netcdf output queries
[1321]	24	!
	25	! Former revisions:
	26	! -----------------
	27	! $Id: data_output_3d.f90 1327 2014-03-21 11:00:16Z maronga $
	28	!
	29	! 1320 2014-03-20 08:40:49Z raasch
[1320]	30	! ONLY-attribute added to USE-statements,
	31	! kind-parameters added to all INTEGER and REAL declaration statements,
	32	! kinds are defined in new module kinds,
	33	! old module precision_kind is removed,
	34	! revision history before 2012 removed,
	35	! comment fields (!:) to be used for variable explanations added to
	36	! all variable declaration statements
[674]	37	!
[1319]	38	! 1318 2014-03-17 13:35:16Z raasch
	39	! barrier argument removed from cpu_log,
	40	! module interfaces removed
	41	!
[1309]	42	! 1308 2014-03-13 14:58:42Z fricke
	43	! Check, if the limit of the time dimension is exceeded for parallel output
	44	! To increase the performance for parallel output, the following is done:
	45	! - Update of time axis is only done by PE0
	46	!
[1245]	47	! 1244 2013-10-31 08:16:56Z raasch
	48	! Bugfix for index bounds in case of 3d-parallel output
	49	!
[1116]	50	! 1115 2013-03-26 18:16:16Z hoffmann
	51	! ql is calculated by calc_liquid_water_content
	52	!
[1107]	53	! 1106 2013-03-04 05:31:38Z raasch
	54	! array_kind renamed precision_kind
	55	!
[1077]	56	! 1076 2012-12-05 08:30:18Z hoffmann
	57	! Bugfix in output of ql
	58	!
[1054]	59	! 1053 2012-11-13 17:11:03Z hoffmann
	60	! +nr, qr, prr, qc and averaged quantities
	61	!
[1037]	62	! 1036 2012-10-22 13:43:42Z raasch
	63	! code put under GPL (PALM 3.9)
	64	!
[1035]	65	! 1031 2012-10-19 14:35:30Z raasch
	66	! netCDF4 without parallel file support implemented
	67	!
[1008]	68	! 1007 2012-09-19 14:30:36Z franke
	69	! Bugfix: missing calculation of ql_vp added
	70	!
[1]	71	! Revision 1.1 1997/09/03 06:29:36 raasch
	72	! Initial revision
	73	!
	74	!
	75	! Description:
	76	! ------------
[1031]	77	! Output of the 3D-arrays in netCDF and/or AVS format.
[1]	78	!------------------------------------------------------------------------------!
	79
[1320]	80	USE arrays_3d, &
	81	ONLY: e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, sa, tend, u, v, &
	82	vpt, w
	83
[1]	84	USE averaging
[1320]	85
	86	USE cloud_parameters, &
	87	ONLY: l_d_cp, prr, pt_d_t
	88
	89	USE control_parameters, &
[1327]	90	ONLY: avs_data_file, cloud_physics, do3d, do3d_avs_n, &
	91	do3d_no, do3d_time_count, io_blocks, io_group, &
	92	message_string, netcdf_data_format, ntdim_3d, &
[1320]	93	nz_do3d, plot_3d_precision, psolver, simulated_time, &
	94	simulated_time_chr, skip_do_avs, time_since_reference_point
	95
	96	USE cpulog, &
	97	ONLY: log_point, cpu_log
	98
	99	USE indices, &
	100	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzt, &
	101	nzb
	102
	103	USE kinds
	104
[1]	105	USE netcdf_control
[1320]	106
	107	USE particle_attributes, &
	108	ONLY: particles, prt_count, prt_start_index
	109
[1]	110	USE pegrid
	111
	112	IMPLICIT NONE
	113
[1320]	114	CHARACTER (LEN=9) :: simulated_time_mod !:
[1]	115
[1320]	116	INTEGER(iwp) :: av !:
	117	INTEGER(iwp) :: i !:
	118	INTEGER(iwp) :: if !:
	119	INTEGER(iwp) :: j !:
	120	INTEGER(iwp) :: k !:
	121	INTEGER(iwp) :: n !:
	122	INTEGER(iwp) :: pos !:
	123	INTEGER(iwp) :: prec !:
	124	INTEGER(iwp) :: psi !:
[1]	125
[1320]	126	LOGICAL :: found !:
	127	LOGICAL :: resorted !:
[1]	128
[1320]	129	REAL(wp) :: mean_r !:
	130	REAL(wp) :: s_r3 !:
	131	REAL(wp) :: s_r4 !:
[1]	132
[1320]	133	REAL(sp), DIMENSION(:,:,:), ALLOCATABLE :: local_pf !:
[1]	134
[1320]	135	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !:
[1]	136
	137	!
	138	!-- Return, if nothing to output
	139	IF ( do3d_no(av) == 0 ) RETURN
[1308]	140	!
	141	!-- For parallel netcdf output the time axis must be limited. Return, if this
	142	!-- limit is exceeded. This could be the case, if the simulated time exceeds
	143	!-- the given end time by the length of the given output interval.
	144	IF ( netcdf_data_format > 4 ) THEN
	145	IF ( do3d_time_count(av) + 1 > ntdim_3d(av) ) THEN
	146	WRITE ( message_string, * ) 'Output of 3d data is not given at t=', &
	147	simulated_time, '&because the maximum ', &
	148	'number of output time levels is ', &
	149	'exceeded.'
	150	CALL message( 'data_output_3d', 'PA0387', 0, 1, 0, 6, 0 )
	151	RETURN
	152	ENDIF
	153	ENDIF
[1]	154
	155	CALL cpu_log (log_point(14),'data_output_3d','start')
	156
	157	!
	158	!-- Open output file.
	159	!-- Also creates coordinate and fld-file for AVS.
[1031]	160	!-- For classic or 64bit netCDF output or output of other (old) data formats,
[493]	161	!-- for a run on more than one PE, each PE opens its own file and
[1]	162	!-- writes the data of its subdomain in binary format (regardless of the format
	163	!-- the user has requested). After the run, these files are combined to one
	164	!-- file by combine_plot_fields in the format requested by the user (netcdf
[493]	165	!-- and/or avs).
[1031]	166	!-- For netCDF4/HDF5 output, data is written in parallel into one file.
[1327]	167	IF ( netcdf_data_format < 5 ) THEN
	168	CALL check_open( 30 )
	169	IF ( myid == 0 ) CALL check_open( 106+av*10 )
[493]	170	ELSE
[1327]	171	CALL check_open( 106+av*10 )
[493]	172	ENDIF
[1]	173
	174	!
	175	!-- Allocate a temporary array with the desired output dimensions.
[667]	176	ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) )
[1]	177
	178	!
[1031]	179	!-- Update the netCDF time axis
[1308]	180	!-- In case of parallel output, this is only done by PE0 to increase the
	181	!-- performance.
[1]	182	#if defined( __netcdf )
[1308]	183	do3d_time_count(av) = do3d_time_count(av) + 1
	184	IF ( myid == 0 ) THEN
[1327]	185	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av), &
	186	(/ time_since_reference_point /), &
	187	start = (/ do3d_time_count(av) /), &
	188	count = (/ 1 /) )
	189	CALL handle_netcdf_error( 'data_output_3d', 376 )
[1]	190	ENDIF
	191	#endif
	192
	193	!
	194	!-- Loop over all variables to be written.
	195	if = 1
	196
	197	DO WHILE ( do3d(av,if)(1:1) /= ' ' )
	198	!
	199	!-- Store the array chosen on the temporary array.
	200	resorted = .FALSE.
	201	SELECT CASE ( TRIM( do3d(av,if) ) )
	202
	203	CASE ( 'e' )
	204	IF ( av == 0 ) THEN
	205	to_be_resorted => e
	206	ELSE
	207	to_be_resorted => e_av
	208	ENDIF
	209
[771]	210	CASE ( 'lpt' )
	211	IF ( av == 0 ) THEN
	212	to_be_resorted => pt
	213	ELSE
	214	to_be_resorted => lpt_av
	215	ENDIF
	216
[1053]	217	CASE ( 'nr' )
	218	IF ( av == 0 ) THEN
	219	to_be_resorted => nr
	220	ELSE
	221	to_be_resorted => nr_av
	222	ENDIF
	223
[1]	224	CASE ( 'p' )
	225	IF ( av == 0 ) THEN
[727]	226	IF ( psolver /= 'sor' ) CALL exchange_horiz( p, nbgp )
[1]	227	to_be_resorted => p
	228	ELSE
[727]	229	IF ( psolver /= 'sor' ) CALL exchange_horiz( p_av, nbgp )
[1]	230	to_be_resorted => p_av
	231	ENDIF
	232
	233	CASE ( 'pc' ) ! particle concentration (requires ghostpoint exchange)
	234	IF ( av == 0 ) THEN
	235	tend = prt_count
[667]	236	CALL exchange_horiz( tend, nbgp )
	237	DO i = nxlg, nxrg
	238	DO j = nysg, nyng
[1]	239	DO k = nzb, nz_do3d
	240	local_pf(i,j,k) = tend(k,j,i)
	241	ENDDO
	242	ENDDO
	243	ENDDO
	244	resorted = .TRUE.
	245	ELSE
[667]	246	CALL exchange_horiz( pc_av, nbgp )
[1]	247	to_be_resorted => pc_av
	248	ENDIF
	249
	250	CASE ( 'pr' ) ! mean particle radius
	251	IF ( av == 0 ) THEN
	252	DO i = nxl, nxr
	253	DO j = nys, nyn
[790]	254	DO k = nzb, nz_do3d
[1]	255	psi = prt_start_index(k,j,i)
	256	s_r3 = 0.0
	257	s_r4 = 0.0
	258	DO n = psi, psi+prt_count(k,j,i)-1
[1320]	259	s_r3 = s_r3 + particles(n)%radius*3 &
[790]	260	particles(n)%weight_factor
[1320]	261	s_r4 = s_r4 + particles(n)%radius*4 &
[790]	262	particles(n)%weight_factor
[1]	263	ENDDO
	264	IF ( s_r3 /= 0.0 ) THEN
	265	mean_r = s_r4 / s_r3
	266	ELSE
	267	mean_r = 0.0
	268	ENDIF
	269	tend(k,j,i) = mean_r
	270	ENDDO
	271	ENDDO
	272	ENDDO
[667]	273	CALL exchange_horiz( tend, nbgp )
	274	DO i = nxlg, nxrg
	275	DO j = nysg, nyng
[790]	276	DO k = nzb, nz_do3d
[1]	277	local_pf(i,j,k) = tend(k,j,i)
	278	ENDDO
	279	ENDDO
	280	ENDDO
	281	resorted = .TRUE.
	282	ELSE
[667]	283	CALL exchange_horiz( pr_av, nbgp )
[1]	284	to_be_resorted => pr_av
	285	ENDIF
	286
[1053]	287	CASE ( 'prr' )
	288	IF ( av == 0 ) THEN
	289	CALL exchange_horiz( prr, nbgp )
	290	DO i = nxlg, nxrg
	291	DO j = nysg, nyng
	292	DO k = nzb, nzt+1
	293	local_pf(i,j,k) = prr(k,j,i)
	294	ENDDO
	295	ENDDO
	296	ENDDO
	297	ELSE
	298	CALL exchange_horiz( prr_av, nbgp )
	299	DO i = nxlg, nxrg
	300	DO j = nysg, nyng
	301	DO k = nzb, nzt+1
	302	local_pf(i,j,k) = prr_av(k,j,i)
	303	ENDDO
	304	ENDDO
	305	ENDDO
	306	ENDIF
	307	resorted = .TRUE.
	308
[1]	309	CASE ( 'pt' )
	310	IF ( av == 0 ) THEN
	311	IF ( .NOT. cloud_physics ) THEN
	312	to_be_resorted => pt
	313	ELSE
[667]	314	DO i = nxlg, nxrg
	315	DO j = nysg, nyng
[1]	316	DO k = nzb, nz_do3d
[1320]	317	local_pf(i,j,k) = pt(k,j,i) + l_d_cp * &
	318	pt_d_t(k) * &
[1]	319	ql(k,j,i)
	320	ENDDO
	321	ENDDO
	322	ENDDO
	323	resorted = .TRUE.
	324	ENDIF
	325	ELSE
	326	to_be_resorted => pt_av
	327	ENDIF
	328
	329	CASE ( 'q' )
	330	IF ( av == 0 ) THEN
	331	to_be_resorted => q
	332	ELSE
	333	to_be_resorted => q_av
	334	ENDIF
[691]	335
[1053]	336	CASE ( 'qc' )
[1]	337	IF ( av == 0 ) THEN
[1115]	338	to_be_resorted => qc
[1]	339	ELSE
[1115]	340	to_be_resorted => qc_av
[1]	341	ENDIF
	342
[1053]	343	CASE ( 'ql' )
	344	IF ( av == 0 ) THEN
[1115]	345	to_be_resorted => ql
[1053]	346	ELSE
[1115]	347	to_be_resorted => ql_av
[1053]	348	ENDIF
	349
[1]	350	CASE ( 'ql_c' )
	351	IF ( av == 0 ) THEN
	352	to_be_resorted => ql_c
	353	ELSE
	354	to_be_resorted => ql_c_av
	355	ENDIF
	356
	357	CASE ( 'ql_v' )
	358	IF ( av == 0 ) THEN
	359	to_be_resorted => ql_v
	360	ELSE
	361	to_be_resorted => ql_v_av
	362	ENDIF
	363
	364	CASE ( 'ql_vp' )
	365	IF ( av == 0 ) THEN
[1007]	366	DO i = nxl, nxr
	367	DO j = nys, nyn
	368	DO k = nzb, nz_do3d
	369	psi = prt_start_index(k,j,i)
	370	DO n = psi, psi+prt_count(k,j,i)-1
[1320]	371	tend(k,j,i) = tend(k,j,i) + &
	372	particles(n)%weight_factor / &
[1007]	373	prt_count(k,j,i)
	374	ENDDO
	375	ENDDO
	376	ENDDO
	377	ENDDO
	378	CALL exchange_horiz( tend, nbgp )
	379	DO i = nxlg, nxrg
	380	DO j = nysg, nyng
	381	DO k = nzb, nz_do3d
	382	local_pf(i,j,k) = tend(k,j,i)
	383	ENDDO
	384	ENDDO
	385	ENDDO
	386	resorted = .TRUE.
[1]	387	ELSE
[1007]	388	CALL exchange_horiz( ql_vp_av, nbgp )
[1]	389	to_be_resorted => ql_vp_av
	390	ENDIF
	391
[1053]	392	CASE ( 'qr' )
	393	IF ( av == 0 ) THEN
	394	to_be_resorted => qr
	395	ELSE
	396	to_be_resorted => qr_av
	397	ENDIF
	398
[1]	399	CASE ( 'qv' )
	400	IF ( av == 0 ) THEN
[667]	401	DO i = nxlg, nxrg
	402	DO j = nysg, nyng
[1]	403	DO k = nzb, nz_do3d
	404	local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
	405	ENDDO
	406	ENDDO
	407	ENDDO
	408	resorted = .TRUE.
	409	ELSE
	410	to_be_resorted => qv_av
	411	ENDIF
	412
[96]	413	CASE ( 'rho' )
	414	IF ( av == 0 ) THEN
	415	to_be_resorted => rho
	416	ELSE
	417	to_be_resorted => rho_av
	418	ENDIF
[691]	419
[1]	420	CASE ( 's' )
	421	IF ( av == 0 ) THEN
	422	to_be_resorted => q
	423	ELSE
[355]	424	to_be_resorted => s_av
[1]	425	ENDIF
[691]	426
[96]	427	CASE ( 'sa' )
	428	IF ( av == 0 ) THEN
	429	to_be_resorted => sa
	430	ELSE
	431	to_be_resorted => sa_av
	432	ENDIF
[691]	433
[1]	434	CASE ( 'u' )
	435	IF ( av == 0 ) THEN
	436	to_be_resorted => u
	437	ELSE
	438	to_be_resorted => u_av
	439	ENDIF
	440
	441	CASE ( 'v' )
	442	IF ( av == 0 ) THEN
	443	to_be_resorted => v
	444	ELSE
	445	to_be_resorted => v_av
	446	ENDIF
	447
	448	CASE ( 'vpt' )
	449	IF ( av == 0 ) THEN
	450	to_be_resorted => vpt
	451	ELSE
	452	to_be_resorted => vpt_av
	453	ENDIF
	454
	455	CASE ( 'w' )
	456	IF ( av == 0 ) THEN
	457	to_be_resorted => w
	458	ELSE
	459	to_be_resorted => w_av
	460	ENDIF
	461
	462	CASE DEFAULT
	463	!
	464	!-- User defined quantity
[1320]	465	CALL user_data_output_3d( av, do3d(av,if), found, local_pf, &
[1]	466	nz_do3d )
	467	resorted = .TRUE.
	468
[254]	469	IF ( .NOT. found ) THEN
[1320]	470	message_string = 'no output available for: ' // &
[274]	471	TRIM( do3d(av,if) )
[254]	472	CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
[1]	473	ENDIF
	474
	475	END SELECT
	476
	477	!
	478	!-- Resort the array to be output, if not done above
	479	IF ( .NOT. resorted ) THEN
[667]	480	DO i = nxlg, nxrg
	481	DO j = nysg, nyng
[1]	482	DO k = nzb, nz_do3d
	483	local_pf(i,j,k) = to_be_resorted(k,j,i)
	484	ENDDO
	485	ENDDO
	486	ENDDO
	487	ENDIF
	488
	489	!
[1327]	490	!-- Output of the 3D-array
	491	#if defined( __parallel )
	492	IF ( netcdf_data_format < 5 ) THEN
[1]	493	!
[1327]	494	!-- Non-parallel netCDF output. Data is output in parallel in
	495	!-- FORTRAN binary format here, and later collected into one file by
	496	!-- combine_plot_fields
	497	IF ( myid == 0 ) THEN
	498	WRITE ( 30 ) time_since_reference_point, &
	499	do3d_time_count(av), av
[1]	500	ENDIF
[1327]	501	DO i = 0, io_blocks-1
	502	IF ( i == io_group ) THEN
	503	WRITE ( 30 ) nxlg, nxrg, nysg, nyng, nzb, nz_do3d
	504	WRITE ( 30 ) local_pf
	505	ENDIF
[1]	506	#if defined( __parallel )
[1327]	507	CALL MPI_BARRIER( comm2d, ierr )
[759]	508	#endif
[1327]	509	ENDDO
[559]	510
[1327]	511	ELSE
[646]	512	#if defined( __netcdf )
[493]	513	!
[1327]	514	!-- Parallel output in netCDF4/HDF5 format.
	515	!-- Do not output redundant ghost point data except for the
	516	!-- boundaries of the total domain.
	517	IF ( nxr == nx .AND. nyn /= ny ) THEN
	518	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
	519	local_pf(nxl:nxr+1,nys:nyn,nzb:nz_do3d), &
	520	start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
	521	count = (/ nxr-nxl+2, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
	522	ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
	523	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
	524	local_pf(nxl:nxr,nys:nyn+1,nzb:nz_do3d), &
	525	start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
	526	count = (/ nxr-nxl+1, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
	527	ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
	528	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
	529	local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d), &
	530	start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
	531	count = (/ nxr-nxl+2, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
	532	ELSE
	533	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
	534	local_pf(nxl:nxr,nys:nyn,nzb:nz_do3d), &
	535	start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
	536	count = (/ nxr-nxl+1, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
	537	ENDIF
	538	CALL handle_netcdf_error( 'data_output_3d', 386 )
[646]	539	#endif
[1327]	540	ENDIF
[1]	541	#else
	542	#if defined( __netcdf )
[1327]	543	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
	544	local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d), &
	545	start = (/ 1, 1, 1, do3d_time_count(av) /), &
	546	count = (/ nx+2, ny+2, nz_do3d-nzb+1, 1 /) )
	547	CALL handle_netcdf_error( 'data_output_3d', 446 )
[1]	548	#endif
	549	#endif
	550
	551	if = if + 1
	552
	553	ENDDO
	554
	555	!
	556	!-- Deallocate temporary array.
	557	DEALLOCATE( local_pf )
	558
	559
[1318]	560	CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
[1]	561
	562	!
	563	!-- Formats.
[1320]	564	3300 FORMAT ('variable ',I4,' file=',A,' filetype=unformatted skip=',I12/ &
[1]	565	'label = ',A,A)
	566
	567	END SUBROUTINE data_output_3d

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

Context Navigation

source: palm/trunk/SOURCE/data_output_3d.f90 @ 1328

Download in other formats: