Home

Context Navigation

source: palm/trunk/SOURCE/data_output_mask.f90 @ 3654

Last change on this file since 3654 was 3632, checked in by kanani, 6 years ago
Fix for masked output, change default chemistry mechanism
Property svn:keywords set to `Id`
File size: 33.5 KB

Line
1	!> @file data_output_mask.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-2018 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: data_output_mask.f90 3632 2018-12-17 14:52:50Z suehring $
27	! Fix output time levels (use time_since_reference_point)
28	!
29	! 3589 2018-11-30 15:09:51Z suehring
30	! Move the control parameter "salsa" from salsa_mod.f90 to control_parameters
31	!
32	! 3582 2018-11-29 19:16:36Z suehring
33	! add variable description
34	!
35	! 3467 2018-10-30 19:05:21Z suehring
36	! Implementation of a new aerosol module salsa.
37	!
38	! 3435 2018-10-26 18:25:44Z gronemeier
39	! Add terrain-following output
40	!
41	! 3421 2018-10-24 18:39:32Z gronemeier
42	! Renamed output variables
43	!
44	! 3419 2018-10-24 17:27:31Z gronemeier
45	! Minor formatting (kanani)
46	! Call for chem_data_output_mask added (basit)
47	!
48	! 3274 2018-09-24 15:42:55Z knoop
49	! Modularization of all bulk cloud physics code components
50	!
51	! 3241 2018-09-12 15:02:00Z raasch
52	! unused variable removed
53	!
54	! 3083 2018-06-19 14:03:12Z gronemeier
55	!
56	!
57	! 3045 2018-05-28 07:55:41Z Giersch
58	! Error messages revised
59	!
60	! 3030 2018-05-23 14:37:00Z raasch
61	! variable if renamed ivar
62	!
63	! 2718 2018-01-02 08:49:38Z maronga
64	! Corrected "Former revisions" section
65	!
66	! 2696 2017-12-14 17:12:51Z kanani
67	! Change in file header (GPL part)
68	!
69	! 2292 2017-06-20 09:51:42Z schwenkel
70	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
71	! includes two more prognostic equations for cloud drop concentration (nc)
72	! and cloud water content (qc).
73	!
74	! 2101 2017-01-05 16:42:31Z suehring
75	!
76	! 2031 2016-10-21 15:11:58Z knoop
77	! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
78	!
79	! 2000 2016-08-20 18:09:15Z knoop
80	! Forced header and separation lines into 80 columns
81	!
82	! 1980 2016-07-29 15:51:57Z suehring
83	! Bugfix, in order to steer user-defined output, setting flag found explicitly
84	! to .F.
85	!
86	! 1976 2016-07-27 13:28:04Z maronga
87	! Output of radiation quantities is now done directly in the respective module
88	!
89	! 1960 2016-07-12 16:34:24Z suehring
90	! Separate humidity and passive scalar
91	!
92	! 2016-03-06 18:36:17Z raasch
93	! name change of netcdf routines and module + related changes,
94	! switch back of netcdf data format moved from time integration routine to here
95	!
96	! 1691 2015-10-26 16:17:44Z maronga
97	! Added output of radiative heating rates for RRTMG
98	!
99	! 1682 2015-10-07 23:56:08Z knoop
100	! Code annotations made doxygen readable
101	!
102	! 1585 2015-04-30 07:05:52Z maronga
103	! Added support for RRTMG
104	!
105	! 1438 2014-07-22 14:14:06Z heinze
106	! +nr, qc, qr
107	!
108	! 1359 2014-04-11 17:15:14Z hoffmann
109	! New particle structure integrated.
110	!
111	! 1353 2014-04-08 15:21:23Z heinze
112	! REAL constants provided with KIND-attribute
113	!
114	! 1327 2014-03-21 11:00:16Z raasch
115	!
116	!
117	! 1320 2014-03-20 08:40:49Z raasch
118	! ONLY-attribute added to USE-statements,
119	! kind-parameters added to all INTEGER and REAL declaration statements,
120	! kinds are defined in new module kinds,
121	! revision history before 2012 removed,
122	! comment fields (!:) to be used for variable explanations added to
123	! all variable declaration statements
124	!
125	! 1318 2014-03-17 13:35:16Z raasch
126	! barrier argument removed from cpu_log,
127	! module interfaces removed
128	!
129	! 1092 2013-02-02 11:24:22Z raasch
130	! unused variables removed
131	!
132	! 1036 2012-10-22 13:43:42Z raasch
133	! code put under GPL (PALM 3.9)
134	!
135	! 1031 2012-10-19 14:35:30Z raasch
136	! netCDF4 without parallel file support implemented
137	!
138	! 1007 2012-09-19 14:30:36Z franke
139	! Bugfix: calculation of pr must depend on the particle weighting factor,
140	! missing calculation of ql_vp added
141	!
142	! 410 2009-12-04 17:05:40Z letzel
143	! Initial version
144	!
145	! Description:
146	! ------------
147	!> Masked data output in netCDF format for current mask (current value of mid).
148	!------------------------------------------------------------------------------!
149	SUBROUTINE data_output_mask( av )
150
151
152
153	#if defined( __netcdf )
154	USE arrays_3d, &
155	ONLY: e, nc, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, &
156	tend, u, v, vpt, w, d_exner
157
158	USE averaging, &
159	ONLY: e_av, lpt_av, nc_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av, &
160	qc_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av, &
161	rho_ocean_av, s_av, sa_av, u_av, v_av, vpt_av, w_av
162
163	USE basic_constants_and_equations_mod, &
164	ONLY: lv_d_cp
165
166	USE chemistry_model_mod, &
167	ONLY: chem_data_output_mask
168
169	USE control_parameters, &
170	ONLY: air_chemistry, domask, domask_no, domask_time_count, mask_i, &
171	mask_j, mask_k, mask_size, mask_size_l, mask_start_l, &
172	mask_surface, &
173	max_masks, message_string, mid, nz_do3d, salsa, &
174	time_since_reference_point
175
176	USE cpulog, &
177	ONLY: cpu_log, log_point
178
179	USE indices, &
180	ONLY: nbgp, nxl, nxr, nyn, nys, nzb, nzt
181
182	USE kinds
183
184	USE bulk_cloud_model_mod, &
185	ONLY: bulk_cloud_model
186
187	USE NETCDF
188
189	USE netcdf_interface, &
190	ONLY: id_set_mask, id_var_domask, id_var_time_mask, nc_stat, &
191	netcdf_data_format, netcdf_handle_error
192
193	USE particle_attributes, &
194	ONLY: grid_particles, number_of_particles, particles, &
195	particle_advection_start, prt_count
196
197	USE pegrid
198
199	USE radiation_model_mod, &
200	ONLY: radiation, radiation_data_output_mask
201
202	USE salsa_mod, &
203	ONLY: salsa_data_output_mask
204
205	USE surface_mod, &
206	ONLY : surf_def_h, surf_lsm_h, surf_usm_h, get_topography_top_index_ji
207
208	IMPLICIT NONE
209
210	CHARACTER(LEN=5) :: grid !< flag to distinquish between staggered grids
211
212	INTEGER(iwp) :: av !< flag for (non-)average output
213	INTEGER(iwp) :: ngp !< number of grid points of an output slice
214	INTEGER(iwp) :: i !< loop index
215	INTEGER(iwp) :: ivar !< variable index
216	INTEGER(iwp) :: j !< loop index
217	INTEGER(iwp) :: k !< loop index
218	INTEGER(iwp) :: kk !< vertical index
219	INTEGER(iwp) :: n !< loop index
220	INTEGER(iwp) :: netcdf_data_format_save !< value of netcdf_data_format
221	INTEGER(iwp) :: sender !< PE id of sending PE
222	INTEGER(iwp) :: topo_top_ind !< k index of highest horizontal surface
223	INTEGER(iwp) :: ind(6) !< index limits (lower/upper bounds) of array 'local_2d'
224
225	LOGICAL :: found !< true if output variable was found
226	LOGICAL :: resorted !< true if variable is resorted
227
228	REAL(wp) :: mean_r !< mean particle radius
229	REAL(wp) :: s_r2 !< sum( particle-radius**2 )
230	REAL(wp) :: s_r3 !< sum( particle-radius**3 )
231
232	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_pf !< output array
233	#if defined( __parallel )
234	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: total_pf !< collected output array
235	#endif
236	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which shall be output
237
238	!
239	!-- Return, if nothing to output
240	IF ( domask_no(mid,av) == 0 ) RETURN
241
242	CALL cpu_log (log_point(49),'data_output_mask','start')
243
244	!
245	!-- Parallel netcdf output is not tested so far for masked data, hence
246	!-- netcdf_data_format is switched back to non-paralell output.
247	netcdf_data_format_save = netcdf_data_format
248	IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
249	IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
250
251	!
252	!-- Open output file.
253	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
254	CALL check_open( 200+mid+av*max_masks )
255	ENDIF
256
257	!
258	!-- Allocate total and local output arrays.
259	#if defined( __parallel )
260	IF ( myid == 0 ) THEN
261	ALLOCATE( total_pf(mask_size(mid,1),mask_size(mid,2),mask_size(mid,3)) )
262	ENDIF
263	#endif
264	ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
265	mask_size_l(mid,3)) )
266
267	!
268	!-- Update the netCDF time axis.
269	domask_time_count(mid,av) = domask_time_count(mid,av) + 1
270	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
271	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_time_mask(mid,av), &
272	(/ time_since_reference_point /), &
273	start = (/ domask_time_count(mid,av) /), &
274	count = (/ 1 /) )
275	CALL netcdf_handle_error( 'data_output_mask', 460 )
276	ENDIF
277
278	!
279	!-- Loop over all variables to be written.
280	ivar = 1
281
282	DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
283	!
284	!-- Reallocate local_pf on PE 0 since its shape changes during MPI exchange
285	IF ( netcdf_data_format < 5 .AND. myid == 0 .AND. ivar > 1 ) THEN
286	DEALLOCATE( local_pf )
287	ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
288	mask_size_l(mid,3)) )
289	ENDIF
290	!
291	!-- Set default grid for terrain-following output
292	grid = 's'
293	!
294	!-- Set flag to steer output of radiation, land-surface, or user-defined
295	!-- quantities
296	found = .FALSE.
297	!
298	!-- Store the variable chosen.
299	resorted = .FALSE.
300	SELECT CASE ( TRIM( domask(mid,av,ivar) ) )
301
302	CASE ( 'e' )
303	IF ( av == 0 ) THEN
304	to_be_resorted => e
305	ELSE
306	to_be_resorted => e_av
307	ENDIF
308
309	CASE ( 'thetal' )
310	IF ( av == 0 ) THEN
311	to_be_resorted => pt
312	ELSE
313	to_be_resorted => lpt_av
314	ENDIF
315
316	CASE ( 'nc' )
317	IF ( av == 0 ) THEN
318	to_be_resorted => nc
319	ELSE
320	to_be_resorted => nc_av
321	ENDIF
322
323	CASE ( 'nr' )
324	IF ( av == 0 ) THEN
325	to_be_resorted => nr
326	ELSE
327	to_be_resorted => nr_av
328	ENDIF
329
330	CASE ( 'p' )
331	IF ( av == 0 ) THEN
332	to_be_resorted => p
333	ELSE
334	to_be_resorted => p_av
335	ENDIF
336
337	CASE ( 'pc' ) ! particle concentration (requires ghostpoint exchange)
338	IF ( av == 0 ) THEN
339	tend = prt_count
340	CALL exchange_horiz( tend, nbgp )
341	IF ( .NOT. mask_surface(mid) ) THEN
342	DO i = 1, mask_size_l(mid,1)
343	DO j = 1, mask_size_l(mid,2)
344	DO k = 1, mask_size_l(mid,3)
345	local_pf(i,j,k) = tend(mask_k(mid,k), &
346	mask_j(mid,j),mask_i(mid,i))
347	ENDDO
348	ENDDO
349	ENDDO
350	ELSE
351	!
352	!-- Terrain-following masked output
353	DO i = 1, mask_size_l(mid,1)
354	DO j = 1, mask_size_l(mid,2)
355	!
356	!-- Get k index of highest horizontal surface
357	topo_top_ind = &
358	get_topography_top_index_ji( mask_j(mid,j), &
359	mask_i(mid,i), &
360	grid )
361	DO k = 1, mask_size_l(mid,3)
362	kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
363	local_pf(i,j,k) = &
364	tend(kk,mask_j(mid,j),mask_i(mid,i))
365	ENDDO
366	ENDDO
367	ENDDO
368	ENDIF
369	resorted = .TRUE.
370	ELSE
371	CALL exchange_horiz( pc_av, nbgp )
372	to_be_resorted => pc_av
373	ENDIF
374
375	CASE ( 'pr' ) ! mean particle radius (effective radius)
376	IF ( av == 0 ) THEN
377	IF ( time_since_reference_point >= particle_advection_start ) THEN
378	DO i = nxl, nxr
379	DO j = nys, nyn
380	DO k = nzb, nz_do3d
381	number_of_particles = prt_count(k,j,i)
382	IF (number_of_particles <= 0) CYCLE
383	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
384	s_r2 = 0.0_wp
385	s_r3 = 0.0_wp
386	DO n = 1, number_of_particles
387	IF ( particles(n)%particle_mask ) THEN
388	s_r2 = s_r2 + grid_particles(k,j,i)%particles(n)%radius*2 &
389	grid_particles(k,j,i)%particles(n)%weight_factor
390	s_r3 = s_r3 + grid_particles(k,j,i)%particles(n)%radius*3 &
391	grid_particles(k,j,i)%particles(n)%weight_factor
392	ENDIF
393	ENDDO
394	IF ( s_r2 > 0.0_wp ) THEN
395	mean_r = s_r3 / s_r2
396	ELSE
397	mean_r = 0.0_wp
398	ENDIF
399	tend(k,j,i) = mean_r
400	ENDDO
401	ENDDO
402	ENDDO
403	CALL exchange_horiz( tend, nbgp )
404	ELSE
405	tend = 0.0_wp
406	ENDIF
407	IF ( .NOT. mask_surface(mid) ) THEN
408	DO i = 1, mask_size_l(mid,1)
409	DO j = 1, mask_size_l(mid,2)
410	DO k = 1, mask_size_l(mid,3)
411	local_pf(i,j,k) = tend(mask_k(mid,k), &
412	mask_j(mid,j),mask_i(mid,i))
413	ENDDO
414	ENDDO
415	ENDDO
416	ELSE
417	!
418	!-- Terrain-following masked output
419	DO i = 1, mask_size_l(mid,1)
420	DO j = 1, mask_size_l(mid,2)
421	!
422	!-- Get k index of highest horizontal surface
423	topo_top_ind = &
424	get_topography_top_index_ji( mask_j(mid,j), &
425	mask_i(mid,i), &
426	grid )
427	DO k = 1, mask_size_l(mid,3)
428	kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
429	local_pf(i,j,k) = &
430	tend(kk,mask_j(mid,j),mask_i(mid,i))
431	ENDDO
432	ENDDO
433	ENDDO
434	ENDIF
435	resorted = .TRUE.
436	ELSE
437	CALL exchange_horiz( pr_av, nbgp )
438	to_be_resorted => pr_av
439	ENDIF
440
441	CASE ( 'theta' )
442	IF ( av == 0 ) THEN
443	IF ( .NOT. bulk_cloud_model ) THEN
444	to_be_resorted => pt
445	ELSE
446	IF ( .NOT. mask_surface(mid) ) THEN
447	DO i = 1, mask_size_l(mid,1)
448	DO j = 1, mask_size_l(mid,2)
449	DO k = 1, mask_size_l(mid,3)
450	local_pf(i,j,k) = &
451	pt(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) &
452	+ lv_d_cp * d_exner(mask_k(mid,k)) * &
453	ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
454	ENDDO
455	ENDDO
456	ENDDO
457	ELSE
458	!
459	!-- Terrain-following masked output
460	DO i = 1, mask_size_l(mid,1)
461	DO j = 1, mask_size_l(mid,2)
462	!
463	!-- Get k index of highest horizontal surface
464	topo_top_ind = &
465	get_topography_top_index_ji( mask_j(mid,j), &
466	mask_i(mid,i), &
467	grid )
468	DO k = 1, mask_size_l(mid,3)
469	kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
470	local_pf(i,j,k) = &
471	pt(kk,mask_j(mid,j),mask_i(mid,i) ) &
472	+ lv_d_cp * d_exner(kk) * &
473	ql(kk,mask_j(mid,j),mask_i(mid,i))
474	ENDDO
475	ENDDO
476	ENDDO
477	ENDIF
478	resorted = .TRUE.
479	ENDIF
480	ELSE
481	to_be_resorted => pt_av
482	ENDIF
483
484	CASE ( 'q' )
485	IF ( av == 0 ) THEN
486	to_be_resorted => q
487	ELSE
488	to_be_resorted => q_av
489	ENDIF
490
491	CASE ( 'qc' )
492	IF ( av == 0 ) THEN
493	to_be_resorted => qc
494	ELSE
495	to_be_resorted => qc_av
496	ENDIF
497
498	CASE ( 'ql' )
499	IF ( av == 0 ) THEN
500	to_be_resorted => ql
501	ELSE
502	to_be_resorted => ql_av
503	ENDIF
504
505	CASE ( 'ql_c' )
506	IF ( av == 0 ) THEN
507	to_be_resorted => ql_c
508	ELSE
509	to_be_resorted => ql_c_av
510	ENDIF
511
512	CASE ( 'ql_v' )
513	IF ( av == 0 ) THEN
514	to_be_resorted => ql_v
515	ELSE
516	to_be_resorted => ql_v_av
517	ENDIF
518
519	CASE ( 'ql_vp' )
520	IF ( av == 0 ) THEN
521	IF ( time_since_reference_point >= particle_advection_start ) THEN
522	DO i = nxl, nxr
523	DO j = nys, nyn
524	DO k = nzb, nz_do3d
525	number_of_particles = prt_count(k,j,i)
526	IF (number_of_particles <= 0) CYCLE
527	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
528	DO n = 1, number_of_particles
529	IF ( particles(n)%particle_mask ) THEN
530	tend(k,j,i) = tend(k,j,i) + &
531	particles(n)%weight_factor / &
532	prt_count(k,j,i)
533	ENDIF
534	ENDDO
535	ENDDO
536	ENDDO
537	ENDDO
538	CALL exchange_horiz( tend, nbgp )
539	ELSE
540	tend = 0.0_wp
541	ENDIF
542	IF ( .NOT. mask_surface(mid) ) THEN
543	DO i = 1, mask_size_l(mid,1)
544	DO j = 1, mask_size_l(mid,2)
545	DO k = 1, mask_size_l(mid,3)
546	local_pf(i,j,k) = tend(mask_k(mid,k), &
547	mask_j(mid,j),mask_i(mid,i))
548	ENDDO
549	ENDDO
550	ENDDO
551	ELSE
552	!
553	!-- Terrain-following masked output
554	DO i = 1, mask_size_l(mid,1)
555	DO j = 1, mask_size_l(mid,2)
556	!
557	!-- Get k index of highest horizontal surface
558	topo_top_ind = &
559	get_topography_top_index_ji( mask_j(mid,j), &
560	mask_i(mid,i), &
561	grid )
562	DO k = 1, mask_size_l(mid,3)
563	kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
564	local_pf(i,j,k) = &
565	tend(kk,mask_j(mid,j),mask_i(mid,i))
566	ENDDO
567	ENDDO
568	ENDDO
569	ENDIF
570	resorted = .TRUE.
571	ELSE
572	CALL exchange_horiz( ql_vp_av, nbgp )
573	to_be_resorted => ql_vp_av
574	ENDIF
575
576	CASE ( 'qv' )
577	IF ( av == 0 ) THEN
578	IF ( .NOT. mask_surface(mid) ) THEN
579	DO i = 1, mask_size_l(mid,1)
580	DO j = 1, mask_size_l(mid,2)
581	DO k = 1, mask_size_l(mid,3)
582	local_pf(i,j,k) = &
583	q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) - &
584	ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
585	ENDDO
586	ENDDO
587	ENDDO
588	ELSE
589	!
590	!-- Terrain-following masked output
591	DO i = 1, mask_size_l(mid,1)
592	DO j = 1, mask_size_l(mid,2)
593	!
594	!-- Get k index of highest horizontal surface
595	topo_top_ind = &
596	get_topography_top_index_ji( mask_j(mid,j), &
597	mask_i(mid,i), &
598	grid )
599	DO k = 1, mask_size_l(mid,3)
600	kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
601	local_pf(i,j,k) = &
602	q(kk,mask_j(mid,j),mask_i(mid,i)) - &
603	ql(kk,mask_j(mid,j),mask_i(mid,i))
604	ENDDO
605	ENDDO
606	ENDDO
607	ENDIF
608	resorted = .TRUE.
609	ELSE
610	to_be_resorted => qv_av
611	ENDIF
612
613	CASE ( 'qr' )
614	IF ( av == 0 ) THEN
615	to_be_resorted => qr
616	ELSE
617	to_be_resorted => qr_av
618	ENDIF
619
620	CASE ( 'rho_sea_water' )
621	IF ( av == 0 ) THEN
622	to_be_resorted => rho_ocean
623	ELSE
624	to_be_resorted => rho_ocean_av
625	ENDIF
626
627	CASE ( 's' )
628	IF ( av == 0 ) THEN
629	to_be_resorted => s
630	ELSE
631	to_be_resorted => s_av
632	ENDIF
633
634	CASE ( 'sa' )
635	IF ( av == 0 ) THEN
636	to_be_resorted => sa
637	ELSE
638	to_be_resorted => sa_av
639	ENDIF
640
641	CASE ( 'u' )
642	IF ( av == 0 ) THEN
643	to_be_resorted => u
644	ELSE
645	to_be_resorted => u_av
646	ENDIF
647
648	CASE ( 'v' )
649	IF ( av == 0 ) THEN
650	to_be_resorted => v
651	ELSE
652	to_be_resorted => v_av
653	ENDIF
654
655	CASE ( 'thetav' )
656	IF ( av == 0 ) THEN
657	to_be_resorted => vpt
658	ELSE
659	to_be_resorted => vpt_av
660	ENDIF
661
662	CASE ( 'w' )
663	grid = 'w'
664	IF ( av == 0 ) THEN
665	to_be_resorted => w
666	ELSE
667	to_be_resorted => w_av
668	ENDIF
669
670	CASE DEFAULT
671
672	!
673	!-- Radiation quantity
674	IF ( radiation ) THEN
675	CALL radiation_data_output_mask(av, domask(mid,av,ivar), found,&
676	local_pf )
677	ENDIF
678
679	IF ( air_chemistry ) THEN
680	CALL chem_data_output_mask(av, domask(mid,av,ivar), found, &
681	local_pf )
682	ENDIF
683	!
684	!-- SALSA quantities
685	IF ( salsa ) THEN
686	CALL salsa_data_output_mask( av, domask(mid,av,ivar), found, &
687	local_pf )
688	ENDIF
689	!
690	!-- User defined quantity
691	IF ( .NOT. found ) THEN
692	CALL user_data_output_mask(av, domask(mid,av,ivar), found, &
693	local_pf )
694	ENDIF
695
696	resorted = .TRUE.
697
698	IF ( .NOT. found ) THEN
699	WRITE ( message_string, * ) 'no masked output available for: ',&
700	TRIM( domask(mid,av,ivar) )
701	CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
702	ENDIF
703
704	END SELECT
705
706	!
707	!-- Resort the array to be output, if not done above
708	IF ( .NOT. resorted ) THEN
709	IF ( .NOT. mask_surface(mid) ) THEN
710	!
711	!-- Default masked output
712	DO i = 1, mask_size_l(mid,1)
713	DO j = 1, mask_size_l(mid,2)
714	DO k = 1, mask_size_l(mid,3)
715	local_pf(i,j,k) = to_be_resorted(mask_k(mid,k), &
716	mask_j(mid,j),mask_i(mid,i))
717	ENDDO
718	ENDDO
719	ENDDO
720
721	ELSE
722	!
723	!-- Terrain-following masked output
724	DO i = 1, mask_size_l(mid,1)
725	DO j = 1, mask_size_l(mid,2)
726	!
727	!-- Get k index of highest horizontal surface
728	topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), &
729	mask_i(mid,i), &
730	grid )
731	!
732	!-- Save output array
733	DO k = 1, mask_size_l(mid,3)
734	local_pf(i,j,k) = to_be_resorted( &
735	MIN( topo_top_ind+mask_k(mid,k), &
736	nzt+1 ), &
737	mask_j(mid,j), &
738	mask_i(mid,i) )
739	ENDDO
740	ENDDO
741	ENDDO
742
743	ENDIF
744	ENDIF
745
746	!
747	!-- I/O block. I/O methods are implemented
748	!-- (1) for parallel execution
749	!-- a. with netCDF 4 parallel I/O-enabled library
750	!-- b. with netCDF 3 library
751	!-- (2) for serial execution.
752	!-- The choice of method depends on the correct setting of preprocessor
753	!-- directives __parallel and __netcdf4_parallel as well as on the parameter
754	!-- netcdf_data_format.
755	#if defined( __parallel )
756	#if defined( __netcdf4_parallel )
757	IF ( netcdf_data_format > 4 ) THEN
758	!
759	!-- (1) a. Parallel I/O using netCDF 4 (not yet tested)
760	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
761	id_var_domask(mid,av,ivar), local_pf, &
762	start = (/ mask_start_l(mid,1), mask_start_l(mid,2), &
763	mask_start_l(mid,3), domask_time_count(mid,av) /), &
764	count = (/ mask_size_l(mid,1), mask_size_l(mid,2), &
765	mask_size_l(mid,3), 1 /) )
766	CALL netcdf_handle_error( 'data_output_mask', 461 )
767	ELSE
768	#endif
769	!
770	!-- (1) b. Conventional I/O only through PE0
771	!-- PE0 receives partial arrays from all processors of the respective mask
772	!-- and outputs them. Here a barrier has to be set, because otherwise
773	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
774	CALL MPI_BARRIER( comm2d, ierr )
775
776	ngp = mask_size_l(mid,1) * mask_size_l(mid,2) * mask_size_l(mid,3)
777	IF ( myid == 0 ) THEN
778	!
779	!-- Local array can be relocated directly.
780	total_pf( &
781	mask_start_l(mid,1):mask_start_l(mid,1)+mask_size_l(mid,1)-1, &
782	mask_start_l(mid,2):mask_start_l(mid,2)+mask_size_l(mid,2)-1, &
783	mask_start_l(mid,3):mask_start_l(mid,3)+mask_size_l(mid,3)-1 ) &
784	= local_pf
785	!
786	!-- Receive data from all other PEs.
787	DO n = 1, numprocs-1
788	!
789	!-- Receive index limits first, then array.
790	!-- Index limits are received in arbitrary order from the PEs.
791	CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0, &
792	comm2d, status, ierr )
793	!
794	!-- Not all PEs have data for the mask
795	IF ( ind(1) /= -9999 ) THEN
796	ngp = ( ind(2)-ind(1)+1 ) * (ind(4)-ind(3)+1 ) * &
797	( ind(6)-ind(5)+1 )
798	sender = status(MPI_SOURCE)
799	DEALLOCATE( local_pf )
800	ALLOCATE(local_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)))
801	CALL MPI_RECV( local_pf(ind(1),ind(3),ind(5)), ngp, &
802	MPI_REAL, sender, 1, comm2d, status, ierr )
803	total_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)) &
804	= local_pf
805	ENDIF
806	ENDDO
807
808	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
809	id_var_domask(mid,av,ivar), total_pf, &
810	start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
811	count = (/ mask_size(mid,1), mask_size(mid,2), &
812	mask_size(mid,3), 1 /) )
813	CALL netcdf_handle_error( 'data_output_mask', 462 )
814
815	ELSE
816	!
817	!-- If at least part of the mask resides on the PE, send the index
818	!-- limits for the target array, otherwise send -9999 to PE0.
819	IF ( mask_size_l(mid,1) > 0 .AND. mask_size_l(mid,2) > 0 .AND. &
820	mask_size_l(mid,3) > 0 ) &
821	THEN
822	ind(1) = mask_start_l(mid,1)
823	ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
824	ind(3) = mask_start_l(mid,2)
825	ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
826	ind(5) = mask_start_l(mid,3)
827	ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
828	ELSE
829	ind(1) = -9999; ind(2) = -9999
830	ind(3) = -9999; ind(4) = -9999
831	ind(5) = -9999; ind(6) = -9999
832	ENDIF
833	CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
834	!
835	!-- If applicable, send data to PE0.
836	IF ( ind(1) /= -9999 ) THEN
837	CALL MPI_SEND( local_pf(1,1,1), ngp, MPI_REAL, 0, 1, comm2d, &
838	ierr )
839	ENDIF
840	ENDIF
841	!
842	!-- A barrier has to be set, because otherwise some PEs may proceed too
843	!-- fast so that PE0 may receive wrong data on tag 0.
844	CALL MPI_BARRIER( comm2d, ierr )
845	#if defined( __netcdf4_parallel )
846	ENDIF
847	#endif
848	#else
849	!
850	!-- (2) For serial execution of PALM, the single processor (PE0) holds all
851	!-- data and writes them directly to file.
852	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
853	id_var_domask(mid,av,ivar), local_pf, &
854	start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
855	count = (/ mask_size_l(mid,1), mask_size_l(mid,2),&
856	mask_size_l(mid,3), 1 /) )
857	CALL netcdf_handle_error( 'data_output_mask', 463 )
858	#endif
859
860	ivar = ivar + 1
861
862	ENDDO
863
864	!
865	!-- Deallocate temporary arrays.
866	DEALLOCATE( local_pf )
867	#if defined( __parallel )
868	IF ( myid == 0 ) THEN
869	DEALLOCATE( total_pf )
870	ENDIF
871	#endif
872
873	!
874	!-- Switch back to original format given by user (see beginning of this routine)
875	netcdf_data_format = netcdf_data_format_save
876
877	CALL cpu_log( log_point(49), 'data_output_mask', 'stop' )
878	#endif
879
880
881	END SUBROUTINE data_output_mask

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |