Home

Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

data_output_mask.f90

Last change on this file was 4896, checked in by raasch, 3 years ago
small re-formatting to follow the coding standard, typo in file appendix removed, more meaningful variable names assigned, redundant code removed
Property svn:keywords set to `Id`
File size: 33.5 KB

Rev	Line
[4246]	1	!> @file data_output_mask.f90
[4559]	2	!--------------------------------------------------------------------------------------------------!
[4246]	3	! This file is part of the PALM model system.
	4	!
[4559]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[4246]	8	!
[4559]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[4246]	12	!
[4559]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
	15
[4246]	16	!
[4828]	17	! Copyright 1997-2021 Leibniz Universitaet Hannover
[4559]	18	!--------------------------------------------------------------------------------------------------!
[4246]	19	!
	20	! Current revisions:
	21	! -----------------
[4895]	22	!
	23	!
[4246]	24	! Former revisions:
	25	! -----------------
	26	! $Id: data_output_mask.f90 4896 2021-03-03 16:10:18Z banzhafs $
[4896]	27	! small re-formatting to follow the coding standard
	28	!
	29	! 4895 2021-03-03 15:39:08Z suehring
[4895]	30	! Remove offset in terrain-following masked output
	31	!
	32	! 4828 2021-01-05 11:21:41Z Giersch
[4559]	33	! file re-formatted to follow the PALM coding standard
	34	!
	35	! 4457 2020-03-11 14:20:43Z raasch
[4457]	36	! use statement for exchange horiz added
[4559]	37	!
[4457]	38	! 4444 2020-03-05 15:59:50Z raasch
[4444]	39	! bugfix: cpp-directives for serial mode added
	40	!
	41	! 4377 2020-01-15 11:10:51Z gronemeier
[4559]	42	! bugfix: set fill value for output according to wall_flags_total_0 for non-terrain following output
[4377]	43	!
	44	! 4360 2020-01-07 11:25:50Z suehring
[4559]	45	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
	46	! information used in wall_flags_static_0
[4377]	47	!
[4346]	48	! 4331 2019-12-10 18:25:02Z suehring
[4331]	49	! Formatting adjustment
[4377]	50	!
[4331]	51	! 4329 2019-12-10 15:46:36Z motisi
[4329]	52	! Renamed wall_flags_0 to wall_flags_static_0
[4377]	53	!
[4329]	54	! 4246 2019-09-30 09:27:52Z pavelkrc
[4182]	55	! Corrected "Former revisions" section
[4377]	56	!
[4182]	57	! 4168 2019-08-16 13:50:17Z suehring
[4246]	58	! Remove variable grid
[4377]	59	!
[4246]	60	! 4167 2019-08-16 11:01:48Z suehring
[4559]	61	! Changed behaviour of masked output over surface to follow terrain and ignore buildings
	62	! (J.Resler, T.Gronemeier)
[4377]	63	!
[4246]	64	! 4069 2019-07-01 14:05:51Z Giersch
[4559]	65	! Masked output running index mid has been introduced as a local variable to avoid runtime error
	66	! (Loop variable has been modified) in time_integration
[4377]	67	!
[4246]	68	! 4039 2019-06-18 10:32:41Z suehring
	69	! Modularize diagnostic output
[4377]	70	!
[4246]	71	! 3994 2019-05-22 18:08:09Z suehring
	72	! output of turbulence intensity added
[4377]	73	!
[4246]	74	! 3665 2019-01-10 08:28:24Z raasch
	75	! unused variables removed
[4377]	76	!
[4246]	77	! 3655 2019-01-07 16:51:22Z knoop
	78	! Fix output time levels (use time_since_reference_point)
	79	!
	80	! 410 2009-12-04 17:05:40Z letzel
	81	! Initial version
	82	!
	83	! Description:
	84	! ------------
	85	!> Masked data output in netCDF format for current mask (current value of mid).
[4559]	86	!--------------------------------------------------------------------------------------------------!
[4246]	87	SUBROUTINE data_output_mask( av, mid )
	88
	89
[4377]	90
[4246]	91	#if defined( __netcdf )
[4559]	92	USE arrays_3d, &
	93	ONLY: d_exner, e, nc, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, tend, u, v, &
	94	vpt, w
[4246]	95
[4559]	96	USE averaging, &
	97	ONLY: e_av, lpt_av, nc_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av, qc_av, ql_av, ql_c_av, &
	98	ql_v_av, ql_vp_av, qv_av, qr_av, rho_ocean_av, s_av, sa_av, u_av, v_av, vpt_av, w_av
[4246]	99
[4559]	100	USE basic_constants_and_equations_mod, &
[4246]	101	ONLY: lv_d_cp
	102
[4559]	103	USE chemistry_model_mod, &
[4246]	104	ONLY: chem_data_output_mask
	105
[4559]	106	USE control_parameters, &
	107	ONLY: air_chemistry, domask, domask_no, domask_time_count, mask_i, mask_j, mask_k, &
	108	mask_size_l, mask_surface, max_masks, message_string, nz_do3d, salsa, &
[4246]	109	time_since_reference_point
	110
[4444]	111	#if defined( __parallel )
[4559]	112	USE control_parameters, &
[4444]	113	ONLY: mask_size, mask_start_l
	114	#endif
	115
[4559]	116	USE cpulog, &
[4457]	117	ONLY: cpu_log, log_point
	118
[4559]	119	USE diagnostic_output_quantities_mod, &
[4246]	120	ONLY: doq_output_mask
[4377]	121
[4559]	122	USE exchange_horiz_mod, &
[4457]	123	ONLY: exchange_horiz
[4246]	124
[4559]	125	USE indices, &
[4346]	126	ONLY: nbgp, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
[4246]	127
	128	USE kinds
	129
[4559]	130	USE bulk_cloud_model_mod, &
[4246]	131	ONLY: bulk_cloud_model
[4377]	132
[4246]	133	USE NETCDF
[4377]	134
[4559]	135	USE netcdf_interface, &
	136	ONLY: fill_value, id_set_mask, id_var_domask, id_var_time_mask, nc_stat, &
	137	netcdf_data_format, netcdf_handle_error
[4377]	138
[4559]	139	USE particle_attributes, &
	140	ONLY: grid_particles, number_of_particles, particles, particle_advection_start, prt_count
[4377]	141
[4246]	142	USE pegrid
	143
[4559]	144	USE radiation_model_mod, &
[4246]	145	ONLY: radiation, radiation_data_output_mask
	146
[4559]	147	USE salsa_mod, &
[4377]	148	ONLY: salsa_data_output_mask
[4246]	149
	150
	151	IMPLICIT NONE
	152
	153	INTEGER(iwp) :: av !< flag for (non-)average output
[4377]	154	INTEGER(iwp) :: flag_nr !< number of masking flag
[4246]	155	INTEGER(iwp) :: i !< loop index
[4559]	156	INTEGER(iwp) :: im !< loop index for masked variables
[4246]	157	INTEGER(iwp) :: ivar !< variable index
	158	INTEGER(iwp) :: j !< loop index
[4559]	159	INTEGER(iwp) :: jm !< loop index for masked variables
[4246]	160	INTEGER(iwp) :: k !< loop index
	161	INTEGER(iwp) :: kk !< vertical index
[4895]	162	INTEGER(iwp) :: ktt !< k index of lowest non-terrain surface
[4246]	163	INTEGER(iwp) :: mid !< masked output running index
	164	INTEGER(iwp) :: n !< loop index
	165	INTEGER(iwp) :: netcdf_data_format_save !< value of netcdf_data_format
[4444]	166	#if defined( __parallel )
[4559]	167	INTEGER(iwp) :: ind(6) !< index limits (lower/upper bounds) of array 'local_2d'
[4444]	168	INTEGER(iwp) :: ngp !< number of grid points of an output slice
[4246]	169	INTEGER(iwp) :: sender !< PE id of sending PE
[4444]	170	#endif
[4246]	171
	172	LOGICAL :: found !< true if output variable was found
	173	LOGICAL :: resorted !< true if variable is resorted
	174
	175	REAL(wp) :: mean_r !< mean particle radius
	176	REAL(wp) :: s_r2 !< sum( particle-radius**2 )
	177	REAL(wp) :: s_r3 !< sum( particle-radius**3 )
[4377]	178
[4246]	179	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_pf !< output array
	180	#if defined( __parallel )
	181	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: total_pf !< collected output array
	182	#endif
	183	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which shall be output
	184
[4559]	185
[4246]	186	!
	187	!-- Return, if nothing to output
	188	IF ( domask_no(mid,av) == 0 ) RETURN
	189
	190	CALL cpu_log (log_point(49),'data_output_mask','start')
	191
	192	!
[4559]	193	!-- Parallel netcdf output is not tested so far for masked data, hence netcdf_data_format is
	194	!-- switched back to non-paralell output.
[4246]	195	netcdf_data_format_save = netcdf_data_format
	196	IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
	197	IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
	198
	199	!
	200	!-- Open output file.
	201	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
	202	CALL check_open( 200+mid+av*max_masks )
[4377]	203	ENDIF
[4246]	204
	205	!
	206	!-- Allocate total and local output arrays.
	207	#if defined( __parallel )
	208	IF ( myid == 0 ) THEN
[4559]	209	ALLOCATE( total_pf( mask_size(mid,1),mask_size(mid,2),mask_size(mid,3) ) )
[4246]	210	ENDIF
	211	#endif
[4559]	212	ALLOCATE( local_pf( mask_size_l(mid,1),mask_size_l(mid,2), mask_size_l(mid,3) ) )
[4246]	213
	214	!
	215	!-- Update the netCDF time axis.
	216	domask_time_count(mid,av) = domask_time_count(mid,av) + 1
	217	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
[4559]	218	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_time_mask(mid,av), &
	219	(/ time_since_reference_point /), &
	220	start = (/ domask_time_count(mid,av) /), &
[4246]	221	count = (/ 1 /) )
	222	CALL netcdf_handle_error( 'data_output_mask', 460 )
	223	ENDIF
	224
	225	!
	226	!-- Loop over all variables to be written.
	227	ivar = 1
	228
	229	DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
	230	!
	231	!-- Reallocate local_pf on PE 0 since its shape changes during MPI exchange
[4559]	232	IF ( netcdf_data_format < 5 .AND. myid == 0 .AND. ivar > 1 ) THEN
[4246]	233	DEALLOCATE( local_pf )
[4559]	234	ALLOCATE( local_pf( mask_size_l(mid,1),mask_size_l(mid,2), mask_size_l(mid,3) ) )
[4246]	235	ENDIF
	236	!
[4377]	237	!-- Set masking flag for topography for not resorted arrays
	238	flag_nr = 0
	239	!
[4246]	240	!-- Store the variable chosen.
	241	resorted = .FALSE.
	242	SELECT CASE ( TRIM( domask(mid,av,ivar) ) )
	243
	244	CASE ( 'e' )
	245	IF ( av == 0 ) THEN
	246	to_be_resorted => e
	247	ELSE
	248	to_be_resorted => e_av
	249	ENDIF
	250
	251	CASE ( 'thetal' )
	252	IF ( av == 0 ) THEN
	253	to_be_resorted => pt
	254	ELSE
	255	to_be_resorted => lpt_av
	256	ENDIF
	257
	258	CASE ( 'nc' )
	259	IF ( av == 0 ) THEN
	260	to_be_resorted => nc
	261	ELSE
	262	to_be_resorted => nc_av
	263	ENDIF
	264
	265	CASE ( 'nr' )
	266	IF ( av == 0 ) THEN
	267	to_be_resorted => nr
	268	ELSE
	269	to_be_resorted => nr_av
	270	ENDIF
	271
	272	CASE ( 'p' )
	273	IF ( av == 0 ) THEN
	274	to_be_resorted => p
	275	ELSE
	276	to_be_resorted => p_av
	277	ENDIF
	278
	279	CASE ( 'pc' ) ! particle concentration (requires ghostpoint exchange)
	280	IF ( av == 0 ) THEN
	281	tend = prt_count
	282	CALL exchange_horiz( tend, nbgp )
	283	IF ( .NOT. mask_surface(mid) ) THEN
	284	DO i = 1, mask_size_l(mid,1)
	285	DO j = 1, mask_size_l(mid,2)
	286	DO k = 1, mask_size_l(mid,3)
[4559]	287	local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )
[4246]	288	ENDDO
	289	ENDDO
	290	ENDDO
	291	ELSE
	292	!
	293	!-- Terrain-following masked output
	294	DO i = 1, mask_size_l(mid,1)
	295	DO j = 1, mask_size_l(mid,2)
[4895]	296	!
	297	!-- Get k index of the lowest non-terrain grid point
[4246]	298	im = mask_i(mid,i)
	299	jm = mask_j(mid,j)
[4559]	300	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ), &
[4246]	301	DIM = 1 ) - 1
	302	DO k = 1, mask_size_l(mid,3)
[4895]	303	kk = MIN( ktt + mask_k(mid,k) - 1, nzt+1 )
	304	!-- Set value if not in building, else set fill value
[4346]	305	IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) ) THEN
[4246]	306	local_pf(i,j,k) = fill_value
	307	ELSE
	308	local_pf(i,j,k) = tend(kk,jm,im)
	309	ENDIF
	310	ENDDO
	311	ENDDO
	312	ENDDO
	313	ENDIF
	314	resorted = .TRUE.
	315	ELSE
	316	CALL exchange_horiz( pc_av, nbgp )
	317	to_be_resorted => pc_av
	318	ENDIF
	319
	320	CASE ( 'pr' ) ! mean particle radius (effective radius)
	321	IF ( av == 0 ) THEN
	322	IF ( time_since_reference_point >= particle_advection_start ) THEN
	323	DO i = nxl, nxr
	324	DO j = nys, nyn
	325	DO k = nzb, nz_do3d
	326	number_of_particles = prt_count(k,j,i)
[4559]	327	IF ( number_of_particles <= 0 ) CYCLE
[4246]	328	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	329	s_r2 = 0.0_wp
	330	s_r3 = 0.0_wp
	331	DO n = 1, number_of_particles
	332	IF ( particles(n)%particle_mask ) THEN
[4559]	333	s_r2 = s_r2 + grid_particles(k,j,i)%particles(n)%radius*2 &
[4246]	334	grid_particles(k,j,i)%particles(n)%weight_factor
[4559]	335	s_r3 = s_r3 + grid_particles(k,j,i)%particles(n)%radius*3 &
[4246]	336	grid_particles(k,j,i)%particles(n)%weight_factor
	337	ENDIF
	338	ENDDO
	339	IF ( s_r2 > 0.0_wp ) THEN
	340	mean_r = s_r3 / s_r2
	341	ELSE
	342	mean_r = 0.0_wp
	343	ENDIF
	344	tend(k,j,i) = mean_r
	345	ENDDO
	346	ENDDO
	347	ENDDO
	348	CALL exchange_horiz( tend, nbgp )
	349	ELSE
	350	tend = 0.0_wp
	351	ENDIF
	352	IF ( .NOT. mask_surface(mid) ) THEN
	353	DO i = 1, mask_size_l(mid,1)
	354	DO j = 1, mask_size_l(mid,2)
	355	DO k = 1, mask_size_l(mid,3)
[4559]	356	local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )
[4246]	357	ENDDO
	358	ENDDO
	359	ENDDO
	360	ELSE
	361	!
	362	!-- Terrain-following masked output
	363	DO i = 1, mask_size_l(mid,1)
	364	DO j = 1, mask_size_l(mid,2)
[4895]	365	!
	366	!-- Get k index of the lowest non-terrain grid point
[4246]	367	im = mask_i(mid,i)
	368	jm = mask_j(mid,j)
[4559]	369	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
	370	DIM = 1 ) - 1
[4246]	371	DO k = 1, mask_size_l(mid,3)
[4895]	372	kk = MIN( ktt + mask_k(mid,k) - 1, nzt+1 )
	373	!-- Set value if not in building, else set fill value
[4346]	374	IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) ) THEN
[4246]	375	local_pf(i,j,k) = fill_value
	376	ELSE
	377	local_pf(i,j,k) = tend(kk,jm,im)
	378	ENDIF
	379	ENDDO
	380	ENDDO
	381	ENDDO
	382	ENDIF
	383	resorted = .TRUE.
	384	ELSE
	385	CALL exchange_horiz( pr_av, nbgp )
	386	to_be_resorted => pr_av
	387	ENDIF
	388
	389	CASE ( 'theta' )
	390	IF ( av == 0 ) THEN
	391	IF ( .NOT. bulk_cloud_model ) THEN
	392	to_be_resorted => pt
	393	ELSE
	394	IF ( .NOT. mask_surface(mid) ) THEN
	395	DO i = 1, mask_size_l(mid,1)
	396	DO j = 1, mask_size_l(mid,2)
	397	DO k = 1, mask_size_l(mid,3)
[4559]	398	local_pf(i,j,k) = pt( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) ) &
	399	+ lv_d_cp * d_exner( mask_k(mid,k) ) * &
	400	ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i) )
[4246]	401	ENDDO
	402	ENDDO
	403	ENDDO
	404	ELSE
	405	!
	406	!-- Terrain-following masked output
	407	DO i = 1, mask_size_l(mid,1)
	408	DO j = 1, mask_size_l(mid,2)
[4895]	409	!
	410	!-- Get k index of the lowest non-terrain grid point
[4246]	411	im = mask_i(mid,i)
	412	jm = mask_j(mid,j)
[4559]	413	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
	414	DIM = 1 ) - 1
[4246]	415	DO k = 1, mask_size_l(mid,3)
[4895]	416	kk = MIN( ktt + mask_k(mid,k) - 1, nzt+1 )
	417	!
	418	!-- Set value if not in building, else set fill value
[4346]	419	IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) ) THEN
[4246]	420	local_pf(i,j,k) = fill_value
	421	ELSE
[4559]	422	local_pf(i,j,k) = pt(kk,jm,im) + &
	423	lv_d_cp * d_exner(kk) * ql(kk,jm,im)
[4246]	424	ENDIF
	425	ENDDO
	426	ENDDO
	427	ENDDO
[4377]	428	ENDIF
[4246]	429	resorted = .TRUE.
	430	ENDIF
	431	ELSE
	432	to_be_resorted => pt_av
	433	ENDIF
	434
	435	CASE ( 'q' )
	436	IF ( av == 0 ) THEN
	437	to_be_resorted => q
	438	ELSE
	439	to_be_resorted => q_av
	440	ENDIF
	441
	442	CASE ( 'qc' )
	443	IF ( av == 0 ) THEN
	444	to_be_resorted => qc
	445	ELSE
	446	to_be_resorted => qc_av
	447	ENDIF
	448
	449	CASE ( 'ql' )
	450	IF ( av == 0 ) THEN
	451	to_be_resorted => ql
	452	ELSE
	453	to_be_resorted => ql_av
	454	ENDIF
	455
	456	CASE ( 'ql_c' )
	457	IF ( av == 0 ) THEN
	458	to_be_resorted => ql_c
	459	ELSE
	460	to_be_resorted => ql_c_av
	461	ENDIF
	462
	463	CASE ( 'ql_v' )
	464	IF ( av == 0 ) THEN
	465	to_be_resorted => ql_v
	466	ELSE
	467	to_be_resorted => ql_v_av
	468	ENDIF
	469
	470	CASE ( 'ql_vp' )
	471	IF ( av == 0 ) THEN
	472	IF ( time_since_reference_point >= particle_advection_start ) THEN
	473	DO i = nxl, nxr
	474	DO j = nys, nyn
	475	DO k = nzb, nz_do3d
	476	number_of_particles = prt_count(k,j,i)
[4559]	477	IF ( number_of_particles <= 0 ) CYCLE
[4246]	478	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	479	DO n = 1, number_of_particles
	480	IF ( particles(n)%particle_mask ) THEN
	481	tend(k,j,i) = tend(k,j,i) + &
[4559]	482	particles(n)%weight_factor / prt_count(k,j,i)
[4246]	483	ENDIF
	484	ENDDO
	485	ENDDO
	486	ENDDO
	487	ENDDO
	488	CALL exchange_horiz( tend, nbgp )
	489	ELSE
	490	tend = 0.0_wp
	491	ENDIF
	492	IF ( .NOT. mask_surface(mid) ) THEN
	493	DO i = 1, mask_size_l(mid,1)
	494	DO j = 1, mask_size_l(mid,2)
	495	DO k = 1, mask_size_l(mid,3)
[4559]	496	local_pf(i,j,k) = tend( mask_k(mid,k), mask_j(mid,j), mask_i(mid,i) )
[4246]	497	ENDDO
	498	ENDDO
	499	ENDDO
	500	ELSE
	501	!
	502	!-- Terrain-following masked output
	503	DO i = 1, mask_size_l(mid,1)
	504	DO j = 1, mask_size_l(mid,2)
[4895]	505	!
	506	!-- Get k index of the lowest non-terrain grid point
[4246]	507	im = mask_i(mid,i)
	508	jm = mask_j(mid,j)
[4559]	509	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
	510	DIM = 1 ) - 1
[4246]	511	DO k = 1, mask_size_l(mid,3)
[4895]	512	kk = MIN( ktt + mask_k(mid,k) - 1, nzt+1 )
	513	!
	514	!-- Set value if not in building, else set fill value
[4346]	515	IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) ) THEN
[4246]	516	local_pf(i,j,k) = fill_value
	517	ELSE
	518	local_pf(i,j,k) = tend(kk,jm,im)
	519	ENDIF
	520	ENDDO
	521	ENDDO
	522	ENDDO
	523	ENDIF
	524	resorted = .TRUE.
	525	ELSE
	526	CALL exchange_horiz( ql_vp_av, nbgp )
	527	to_be_resorted => ql_vp_av
	528	ENDIF
	529
	530	CASE ( 'qv' )
	531	IF ( av == 0 ) THEN
	532	IF ( .NOT. mask_surface(mid) ) THEN
	533	DO i = 1, mask_size_l(mid,1)
	534	DO j = 1, mask_size_l(mid,2)
	535	DO k = 1, mask_size_l(mid,3)
[4896]	536	local_pf(i,j,k) = q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) - &
	537	ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
[4246]	538	ENDDO
	539	ENDDO
	540	ENDDO
	541	ELSE
	542	!
	543	!-- Terrain-following masked output
	544	DO i = 1, mask_size_l(mid,1)
	545	DO j = 1, mask_size_l(mid,2)
[4895]	546	!
	547	!-- Get k index of the lowest non-terrain grid point
[4246]	548	im = mask_i(mid,i)
	549	jm = mask_j(mid,j)
[4559]	550	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
	551	DIM = 1 ) - 1
[4246]	552	DO k = 1, mask_size_l(mid,3)
[4895]	553	kk = MIN( ktt + mask_k(mid,k) - 1, nzt+1 )
	554	!-- Set value if not in building, else set fill value
[4346]	555	IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) ) THEN
[4246]	556	local_pf(i,j,k) = fill_value
	557	ELSE
	558	local_pf(i,j,k) = q(kk,jm,im) - ql(kk,jm,im)
	559	ENDIF
	560	ENDDO
	561	ENDDO
	562	ENDDO
	563	ENDIF
	564	resorted = .TRUE.
	565	ELSE
	566	to_be_resorted => qv_av
	567	ENDIF
	568
	569	CASE ( 'qr' )
	570	IF ( av == 0 ) THEN
	571	to_be_resorted => qr
	572	ELSE
	573	to_be_resorted => qr_av
	574	ENDIF
	575
	576	CASE ( 'rho_sea_water' )
	577	IF ( av == 0 ) THEN
	578	to_be_resorted => rho_ocean
	579	ELSE
	580	to_be_resorted => rho_ocean_av
	581	ENDIF
	582
	583	CASE ( 's' )
	584	IF ( av == 0 ) THEN
	585	to_be_resorted => s
	586	ELSE
	587	to_be_resorted => s_av
	588	ENDIF
	589
	590	CASE ( 'sa' )
	591	IF ( av == 0 ) THEN
	592	to_be_resorted => sa
	593	ELSE
	594	to_be_resorted => sa_av
	595	ENDIF
	596
	597	CASE ( 'u' )
[4377]	598	flag_nr = 1
[4246]	599	IF ( av == 0 ) THEN
	600	to_be_resorted => u
	601	ELSE
	602	to_be_resorted => u_av
	603	ENDIF
	604
	605	CASE ( 'v' )
[4377]	606	flag_nr = 2
[4246]	607	IF ( av == 0 ) THEN
	608	to_be_resorted => v
	609	ELSE
	610	to_be_resorted => v_av
	611	ENDIF
	612
	613	CASE ( 'thetav' )
	614	IF ( av == 0 ) THEN
	615	to_be_resorted => vpt
	616	ELSE
	617	to_be_resorted => vpt_av
	618	ENDIF
	619
	620	CASE ( 'w' )
[4377]	621	flag_nr = 3
[4246]	622	IF ( av == 0 ) THEN
	623	to_be_resorted => w
	624	ELSE
	625	to_be_resorted => w_av
	626	ENDIF
	627
	628	CASE DEFAULT
	629	!
[4559]	630	!-- Set flag to steer output of radiation, land-surface, or user-defined quantities
[4246]	631	found = .FALSE.
	632	!
	633	!-- Radiation quantity
	634	IF ( .NOT. found .AND. radiation ) THEN
[4559]	635	CALL radiation_data_output_mask( av, domask(mid,av,ivar), found, local_pf, mid )
[4246]	636	ENDIF
	637
	638	IF ( .NOT. found .AND. air_chemistry ) THEN
[4559]	639	CALL chem_data_output_mask( av, domask(mid,av,ivar), found, local_pf, mid )
[4246]	640	ENDIF
	641	!
	642	!-- Check for diagnostic quantities
	643	IF ( .NOT. found ) THEN
[4559]	644	CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf, mid )
[4246]	645	ENDIF
	646	!
	647	!-- SALSA quantities
[4559]	648	IF ( .NOT. found .AND. salsa ) THEN
	649	CALL salsa_data_output_mask( av, domask(mid,av,ivar), found, local_pf, mid )
[4377]	650	ENDIF
[4246]	651	!
	652	!-- User defined quantity
	653	IF ( .NOT. found ) THEN
[4559]	654	CALL user_data_output_mask( av, domask(mid,av,ivar), found, local_pf, mid )
[4246]	655	ENDIF
	656
	657	resorted = .TRUE.
	658
	659	IF ( .NOT. found ) THEN
[4559]	660	WRITE ( message_string, * ) 'no masked output available for: ', &
[4246]	661	TRIM( domask(mid,av,ivar) )
	662	CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
	663	ENDIF
	664
	665	END SELECT
	666
	667	!
	668	!-- Resort the array to be output, if not done above
	669	IF ( .NOT. resorted ) THEN
	670	IF ( .NOT. mask_surface(mid) ) THEN
	671	!
	672	!-- Default masked output
	673	DO i = 1, mask_size_l(mid,1)
	674	DO j = 1, mask_size_l(mid,2)
	675	DO k = 1, mask_size_l(mid,3)
[4559]	676	local_pf(i,j,k) = MERGE( to_be_resorted( mask_k(mid,k), &
	677	mask_j(mid,j), &
	678	mask_i(mid,i)), &
	679	REAL( fill_value, KIND = wp ), &
	680	BTEST( wall_flags_total_0( mask_k(mid,k), &
	681	mask_j(mid,j), &
	682	mask_i(mid,i) ), &
	683	flag_nr ) &
	684	)
[4246]	685	ENDDO
	686	ENDDO
	687	ENDDO
	688
	689	ELSE
	690	!
	691	!-- Terrain-following masked output
	692	DO i = 1, mask_size_l(mid,1)
	693	DO j = 1, mask_size_l(mid,2)
[4895]	694	!
	695	!-- Get k index of the lowest non-terrain grid point
[4246]	696	im = mask_i(mid,i)
	697	jm = mask_j(mid,j)
[4895]	698	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ), &
[4559]	699	DIM = 1 ) - 1
[4246]	700	DO k = 1, mask_size_l(mid,3)
[4895]	701	kk = MIN( ktt + mask_k(mid,k) - 1, nzt+1 )
	702	!
	703	!-- Set value if not in building, else set fill value
[4346]	704	IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) ) THEN
[4246]	705	local_pf(i,j,k) = fill_value
	706	ELSE
	707	local_pf(i,j,k) = to_be_resorted(kk,jm,im)
	708	ENDIF
	709	ENDDO
	710	ENDDO
	711	ENDDO
	712
	713	ENDIF
	714	ENDIF
	715
	716	!
	717	!-- I/O block. I/O methods are implemented
	718	!-- (1) for parallel execution
	719	!-- a. with netCDF 4 parallel I/O-enabled library
	720	!-- b. with netCDF 3 library
	721	!-- (2) for serial execution.
	722	!-- The choice of method depends on the correct setting of preprocessor
	723	!-- directives __parallel and __netcdf4_parallel as well as on the parameter
	724	!-- netcdf_data_format.
	725	#if defined( __parallel )
	726	#if defined( __netcdf4_parallel )
	727	IF ( netcdf_data_format > 4 ) THEN
	728	!
	729	!-- (1) a. Parallel I/O using netCDF 4 (not yet tested)
[4559]	730	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
	731	id_var_domask(mid,av,ivar), local_pf, &
	732	start = (/ mask_start_l(mid,1), mask_start_l(mid,2), &
	733	mask_start_l(mid,3), domask_time_count(mid,av) /), &
	734	count = (/ mask_size_l(mid,1), mask_size_l(mid,2), &
	735	mask_size_l(mid,3), 1 /) )
[4246]	736	CALL netcdf_handle_error( 'data_output_mask', 461 )
	737	ELSE
	738	#endif
	739	!
	740	!-- (1) b. Conventional I/O only through PE0
[4559]	741	!-- PE0 receives partial arrays from all processors of the respective mask and outputs them.
	742	!-- Here a barrier has to be set, because otherwise "-MPI- FATAL: Remote protocol queue full"
	743	!-- may occur.
[4246]	744	CALL MPI_BARRIER( comm2d, ierr )
	745
	746	ngp = mask_size_l(mid,1) * mask_size_l(mid,2) * mask_size_l(mid,3)
	747	IF ( myid == 0 ) THEN
	748	!
	749	!-- Local array can be relocated directly.
[4559]	750	total_pf( mask_start_l(mid,1):mask_start_l(mid,1)+mask_size_l(mid,1)-1, &
	751	mask_start_l(mid,2):mask_start_l(mid,2)+mask_size_l(mid,2)-1, &
	752	mask_start_l(mid,3):mask_start_l(mid,3)+mask_size_l(mid,3)-1 ) &
	753	= local_pf
[4246]	754	!
	755	!-- Receive data from all other PEs.
	756	DO n = 1, numprocs-1
	757	!
	758	!-- Receive index limits first, then array.
	759	!-- Index limits are received in arbitrary order from the PEs.
[4559]	760	CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0, comm2d, status, ierr )
[4246]	761	!
	762	!-- Not all PEs have data for the mask
	763	IF ( ind(1) /= -9999 ) THEN
[4559]	764	ngp = ( ind(2)-ind(1)+1 ) * (ind(4)-ind(3)+1 ) * ( ind(6)-ind(5)+1 )
[4246]	765	sender = status(MPI_SOURCE)
	766	DEALLOCATE( local_pf )
[4559]	767	ALLOCATE( local_pf( ind(1):ind(2),ind(3):ind(4),ind(5):ind(6) ) )
	768	CALL MPI_RECV( local_pf(ind(1),ind(3),ind(5)), ngp, MPI_REAL, sender, 1, comm2d,&
	769	status, ierr )
	770	total_pf( ind(1):ind(2),ind(3):ind(4),ind(5):ind(6) ) = local_pf
[4246]	771	ENDIF
	772	ENDDO
	773
[4559]	774	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
	775	id_var_domask(mid,av,ivar), total_pf, &
	776	start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
	777	count = (/ mask_size(mid,1), mask_size(mid,2), &
	778	mask_size(mid,3), 1 /) )
[4246]	779	CALL netcdf_handle_error( 'data_output_mask', 462 )
	780
	781	ELSE
	782	!
[4559]	783	!-- If at least part of the mask resides on the PE, send the index limits for the target
	784	!-- array, otherwise send -9999 to PE0.
	785	IF ( mask_size_l(mid,1) > 0 .AND. mask_size_l(mid,2) > 0 .AND. &
	786	mask_size_l(mid,3) > 0 ) THEN
[4246]	787	ind(1) = mask_start_l(mid,1)
	788	ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
	789	ind(3) = mask_start_l(mid,2)
	790	ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
	791	ind(5) = mask_start_l(mid,3)
	792	ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
	793	ELSE
	794	ind(1) = -9999; ind(2) = -9999
	795	ind(3) = -9999; ind(4) = -9999
	796	ind(5) = -9999; ind(6) = -9999
	797	ENDIF
	798	CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
	799	!
	800	!-- If applicable, send data to PE0.
	801	IF ( ind(1) /= -9999 ) THEN
[4559]	802	CALL MPI_SEND( local_pf(1,1,1), ngp, MPI_REAL, 0, 1, comm2d, ierr )
[4246]	803	ENDIF
	804	ENDIF
	805	!
[4559]	806	!-- A barrier has to be set, because otherwise some PEs may proceed too fast so that PE0 may
	807	!-- receive wrong data on tag 0.
[4246]	808	CALL MPI_BARRIER( comm2d, ierr )
	809	#if defined( __netcdf4_parallel )
	810	ENDIF
	811	#endif
	812	#else
	813	!
[4559]	814	!-- (2) For serial execution of PALM, the single processor (PE0) holds all data and writes them
	815	!-- directly to file.
	816	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
	817	id_var_domask(mid,av,ivar), local_pf, &
	818	start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
	819	count = (/ mask_size_l(mid,1), mask_size_l(mid,2), &
	820	mask_size_l(mid,3), 1 /) )
[4246]	821	CALL netcdf_handle_error( 'data_output_mask', 463 )
	822	#endif
	823
	824	ivar = ivar + 1
	825
	826	ENDDO
	827
	828	!
	829	!-- Deallocate temporary arrays.
	830	DEALLOCATE( local_pf )
	831	#if defined( __parallel )
	832	IF ( myid == 0 ) THEN
	833	DEALLOCATE( total_pf )
	834	ENDIF
	835	#endif
	836
	837	!
	838	!-- Switch back to original format given by user (see beginning of this routine)
	839	netcdf_data_format = netcdf_data_format_save
	840
	841	CALL cpu_log( log_point(49), 'data_output_mask', 'stop' )
	842	#endif
	843
	844
[4895]	845	END SUBROUTINE data_output_mask

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

Context Navigation

source: palm/trunk/SOURCE/data_output_mask.f90

Download in other formats: