Home

Context Navigation

source: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90 @ 4159

Last change on this file since 4159 was 4157, checked in by suehring, 5 years ago
chem_emissions: Replace global arrays also in mode_emis branch; diagnostic output: restructure initialization in order to work also when data output during spin-up is enabled; radiation: give informative message on raytracing distance only by core zero not by all cores
Property svn:keywords set to `Id`
File size: 38.1 KB

Line
1	!> @file diagnostic_output_quantities_mod.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-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diagnostic_output_quantities_mod.f90 4157 2019-08-14 09:19:12Z suehring $
27	! Initialization restructured, in order to work also when data output during
28	! spin-up is enabled.
29	!
30	! 4132 2019-08-02 12:34:17Z suehring
31	! Bugfix in masked data output
32	!
33	! 4069 2019-07-01 14:05:51Z Giersch
34	! Masked output running index mid has been introduced as a local variable to
35	! avoid runtime error (Loop variable has been modified) in time_integration
36	!
37	! 4039 2019-06-18 10:32:41Z suehring
38	! - Add output of uu, vv, ww to enable variance calculation according temporal
39	! EC method
40	! - Allocate arrays only when they are required
41	! - Formatting adjustment
42	! - Rename subroutines
43	! - Further modularization
44	!
45	! 3998 2019-05-23 13:38:11Z suehring
46	! Bugfix in gathering all output strings
47	!
48	! 3995 2019-05-22 18:59:54Z suehring
49	! Avoid compiler warnings about unused variable and fix string operation which
50	! is not allowed with PGI compiler
51	!
52	! 3994 2019-05-22 18:08:09Z suehring
53	! Initial revision
54	!
55	!
56	! @author Farah Kanani-Suehring
57	!
58	! Description:
59	! ------------
60	!> ...
61	!------------------------------------------------------------------------------!
62	MODULE diagnostic_output_quantities_mod
63
64	USE arrays_3d, &
65	ONLY: ddzu, u, v, w
66	!
67	! USE averaging
68	!
69	! USE basic_constants_and_equations_mod, &
70	! ONLY: c_p, lv_d_cp, l_v
71	!
72	! USE bulk_cloud_model_mod, &
73	! ONLY: bulk_cloud_model
74	!
75	USE control_parameters, &
76	ONLY: current_timestep_number, varnamelength
77	!
78	! USE cpulog, &
79	! ONLY: cpu_log, log_point
80
81	USE grid_variables, &
82	ONLY: ddx, ddy
83	!
84	USE indices, &
85	ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
86	!
87	USE kinds
88	!
89	! USE land_surface_model_mod, &
90	! ONLY: zs
91	!
92	! USE module_interface, &
93	! ONLY: module_interface_data_output_2d
94	!
95	! #if defined( __netcdf )
96	! USE NETCDF
97	! #endif
98	!
99	! USE netcdf_interface, &
100	! ONLY: fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d, &
101	! id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat, &
102	! netcdf_data_format, netcdf_handle_error
103	!
104	! USE particle_attributes, &
105	! ONLY: grid_particles, number_of_particles, particle_advection_start, &
106	! particles, prt_count
107	!
108	! USE pegrid
109	!
110	! USE surface_mod, &
111	! ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h, &
112	! surf_lsm_h, surf_usm_h
113	!
114	! USE turbulence_closure_mod, &
115	! ONLY: tcm_data_output_2d
116
117
118	IMPLICIT NONE
119
120	CHARACTER(LEN=varnamelength), DIMENSION(500) :: do_all = ' '
121
122	INTEGER(iwp) :: timestep_number_at_prev_calc = 0 !< ...at previous diagnostic output calculation
123
124	LOGICAL :: initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
125	LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is p
126
127	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity
128	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity
129
130	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu !< uu
131	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu
132
133	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv !< vv
134	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv
135
136	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww !< ww
137	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww_av !< mean of ww
138
139
140	SAVE
141
142	PRIVATE
143
144	!
145	!-- Public variables
146	PUBLIC do_all, &
147	initialized_diagnostic_output_quantities, &
148	prepared_diagnostic_output_quantities, &
149	timestep_number_at_prev_calc, &
150	ti_av, &
151	uu_av, &
152	vv_av, &
153	ww_av
154	!
155	!-- Public routines
156	PUBLIC doq_3d_data_averaging, &
157	doq_calculate, &
158	doq_check_data_output, &
159	doq_define_netcdf_grid, &
160	doq_output_2d, &
161	doq_output_3d, &
162	doq_output_mask, &
163	doq_init, &
164	doq_wrd_local
165	! doq_rrd_local, &
166
167
168	INTERFACE doq_3d_data_averaging
169	MODULE PROCEDURE doq_3d_data_averaging
170	END INTERFACE doq_3d_data_averaging
171
172	INTERFACE doq_calculate
173	MODULE PROCEDURE doq_calculate
174	END INTERFACE doq_calculate
175
176	INTERFACE doq_check_data_output
177	MODULE PROCEDURE doq_check_data_output
178	END INTERFACE doq_check_data_output
179
180	INTERFACE doq_define_netcdf_grid
181	MODULE PROCEDURE doq_define_netcdf_grid
182	END INTERFACE doq_define_netcdf_grid
183
184	INTERFACE doq_output_2d
185	MODULE PROCEDURE doq_output_2d
186	END INTERFACE doq_output_2d
187
188	INTERFACE doq_output_3d
189	MODULE PROCEDURE doq_output_3d
190	END INTERFACE doq_output_3d
191
192	INTERFACE doq_output_mask
193	MODULE PROCEDURE doq_output_mask
194	END INTERFACE doq_output_mask
195
196	INTERFACE doq_init
197	MODULE PROCEDURE doq_init
198	END INTERFACE doq_init
199
200	INTERFACE doq_prepare
201	MODULE PROCEDURE doq_prepare
202	END INTERFACE doq_prepare
203
204	! INTERFACE doq_rrd_local
205	! MODULE PROCEDURE doq_rrd_local
206	! END INTERFACE doq_rrd_local
207
208	INTERFACE doq_wrd_local
209	MODULE PROCEDURE doq_wrd_local
210	END INTERFACE doq_wrd_local
211
212
213	CONTAINS
214
215	!------------------------------------------------------------------------------!
216	! Description:
217	! ------------
218	!> Sum up and time-average diagnostic output quantities as well as allocate
219	!> the array necessary for storing the average.
220	!------------------------------------------------------------------------------!
221	SUBROUTINE doq_3d_data_averaging( mode, variable )
222
223	USE control_parameters, &
224	ONLY: average_count_3d
225
226	CHARACTER (LEN=*) :: mode !<
227	CHARACTER (LEN=*) :: variable !<
228
229	INTEGER(iwp) :: i !<
230	INTEGER(iwp) :: j !<
231	INTEGER(iwp) :: k !<
232
233	IF ( mode == 'allocate' ) THEN
234
235	SELECT CASE ( TRIM( variable ) )
236
237	CASE ( 'ti' )
238	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
239	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
240	ENDIF
241	ti_av = 0.0_wp
242
243	CASE ( 'uu' )
244	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
245	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
246	ENDIF
247	uu_av = 0.0_wp
248
249	CASE ( 'vv' )
250	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
251	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
252	ENDIF
253	vv_av = 0.0_wp
254
255	CASE ( 'ww' )
256	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
257	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
258	ENDIF
259	ww_av = 0.0_wp
260
261	CASE DEFAULT
262	CONTINUE
263
264	END SELECT
265
266	ELSEIF ( mode == 'sum' ) THEN
267
268	SELECT CASE ( TRIM( variable ) )
269
270	CASE ( 'ti' )
271	IF ( ALLOCATED( ti_av ) ) THEN
272	DO i = nxl, nxr
273	DO j = nys, nyn
274	DO k = nzb, nzt+1
275	ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
276	ENDDO
277	ENDDO
278	ENDDO
279	ENDIF
280
281	CASE ( 'uu' )
282	IF ( ALLOCATED( uu_av ) ) THEN
283	DO i = nxl, nxr
284	DO j = nys, nyn
285	DO k = nzb, nzt+1
286	uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
287	ENDDO
288	ENDDO
289	ENDDO
290	ENDIF
291
292	CASE ( 'vv' )
293	IF ( ALLOCATED( vv_av ) ) THEN
294	DO i = nxl, nxr
295	DO j = nys, nyn
296	DO k = nzb, nzt+1
297	vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
298	ENDDO
299	ENDDO
300	ENDDO
301	ENDIF
302
303	CASE ( 'ww' )
304	IF ( ALLOCATED( ww_av ) ) THEN
305	DO i = nxl, nxr
306	DO j = nys, nyn
307	DO k = nzb, nzt+1
308	ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
309	ENDDO
310	ENDDO
311	ENDDO
312	ENDIF
313
314	CASE DEFAULT
315	CONTINUE
316
317	END SELECT
318
319	ELSEIF ( mode == 'average' ) THEN
320
321	SELECT CASE ( TRIM( variable ) )
322
323	CASE ( 'ti' )
324	IF ( ALLOCATED( ti_av ) ) THEN
325	DO i = nxl, nxr
326	DO j = nys, nyn
327	DO k = nzb, nzt+1
328	ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
329	ENDDO
330	ENDDO
331	ENDDO
332	ENDIF
333
334	CASE ( 'uu' )
335	IF ( ALLOCATED( uu_av ) ) THEN
336	DO i = nxl, nxr
337	DO j = nys, nyn
338	DO k = nzb, nzt+1
339	uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
340	ENDDO
341	ENDDO
342	ENDDO
343	ENDIF
344
345	CASE ( 'vv' )
346	IF ( ALLOCATED( vv_av ) ) THEN
347	DO i = nxl, nxr
348	DO j = nys, nyn
349	DO k = nzb, nzt+1
350	vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
351	ENDDO
352	ENDDO
353	ENDDO
354	ENDIF
355
356	CASE ( 'ww' )
357	IF ( ALLOCATED( ww_av ) ) THEN
358	DO i = nxl, nxr
359	DO j = nys, nyn
360	DO k = nzb, nzt+1
361	ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
362	ENDDO
363	ENDDO
364	ENDDO
365	ENDIF
366
367	END SELECT
368
369	ENDIF
370
371
372	END SUBROUTINE doq_3d_data_averaging
373
374	!------------------------------------------------------------------------------!
375	! Description:
376	! ------------
377	!> Check data output for diagnostic output
378	!------------------------------------------------------------------------------!
379	SUBROUTINE doq_check_data_output( var, unit )
380
381	IMPLICIT NONE
382
383	CHARACTER (LEN=*) :: unit !<
384	CHARACTER (LEN=*) :: var !<
385
386	SELECT CASE ( TRIM( var ) )
387
388	CASE ( 'ti' )
389	unit = '1/s'
390
391	CASE ( 'uu' )
392	unit = 'm2/s2'
393
394	CASE ( 'vv' )
395	unit = 'm2/s2'
396
397	CASE ( 'ww' )
398	unit = 'm2/s2'
399
400	CASE DEFAULT
401	unit = 'illegal'
402
403	END SELECT
404
405
406	END SUBROUTINE doq_check_data_output
407
408	!------------------------------------------------------------------------------!
409	!
410	! Description:
411	! ------------
412	!> Subroutine defining appropriate grid for netcdf variables.
413	!------------------------------------------------------------------------------!
414	SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
415
416	IMPLICIT NONE
417
418	CHARACTER (LEN=*), INTENT(IN) :: variable !<
419	LOGICAL, INTENT(OUT) :: found !<
420	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
421	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
422	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
423
424	found = .TRUE.
425
426	SELECT CASE ( TRIM( variable ) )
427	!
428	!-- s grid
429	CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz' )
430
431	grid_x = 'x'
432	grid_y = 'y'
433	grid_z = 'zu'
434	!
435	!-- u grid
436	CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
437
438	grid_x = 'xu'
439	grid_y = 'y'
440	grid_z = 'zu'
441	!
442	!-- v grid
443	CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz' )
444
445	grid_x = 'x'
446	grid_y = 'yv'
447	grid_z = 'zu'
448	!
449	!-- w grid
450	CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz' )
451
452	grid_x = 'x'
453	grid_y = 'y'
454	grid_z = 'zw'
455
456	CASE DEFAULT
457	found = .FALSE.
458	grid_x = 'none'
459	grid_y = 'none'
460	grid_z = 'none'
461
462	END SELECT
463
464
465	END SUBROUTINE doq_define_netcdf_grid
466
467	!------------------------------------------------------------------------------!
468	!
469	! Description:
470	! ------------
471	!> Subroutine defining 2D output variables
472	!------------------------------------------------------------------------------!
473	SUBROUTINE doq_output_2d( av, variable, found, grid, &
474	mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
475
476
477	IMPLICIT NONE
478
479	CHARACTER (LEN=*) :: grid !<
480	CHARACTER (LEN=*) :: mode !<
481	CHARACTER (LEN=*) :: variable !<
482
483	INTEGER(iwp) :: av !< value indicating averaged or non-averaged output
484	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
485	INTEGER(iwp) :: i !< grid index x-direction
486	INTEGER(iwp) :: j !< grid index y-direction
487	INTEGER(iwp) :: k !< grid index z-direction
488	INTEGER(iwp) :: nzb_do !<
489	INTEGER(iwp) :: nzt_do !<
490
491	LOGICAL :: found !< true if variable is in list
492	LOGICAL :: resorted !< true if array is resorted
493	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
494
495	REAL(wp) :: fill_value !< value for the _FillValue attribute
496
497	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
498	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
499
500	flag_nr = 0
501	found = .TRUE.
502	resorted = .FALSE.
503	two_d = .FALSE.
504
505	SELECT CASE ( TRIM( variable ) )
506
507	CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
508	IF ( av == 0 ) THEN
509	to_be_resorted => ti
510	ELSE
511	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
512	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
513	ti_av = REAL( fill_value, KIND = wp )
514	ENDIF
515	to_be_resorted => ti_av
516	ENDIF
517	flag_nr = 0
518
519	IF ( mode == 'xy' ) grid = 'zu'
520
521	CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
522	IF ( av == 0 ) THEN
523	to_be_resorted => uu
524	ELSE
525	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
526	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
527	uu_av = REAL( fill_value, KIND = wp )
528	ENDIF
529	to_be_resorted => uu_av
530	ENDIF
531	flag_nr = 1
532
533	IF ( mode == 'xy' ) grid = 'zu'
534
535	CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
536	IF ( av == 0 ) THEN
537	to_be_resorted => vv
538	ELSE
539	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
540	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
541	vv_av = REAL( fill_value, KIND = wp )
542	ENDIF
543	to_be_resorted => vv_av
544	ENDIF
545	flag_nr = 2
546
547	IF ( mode == 'xy' ) grid = 'zu'
548
549	CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
550	IF ( av == 0 ) THEN
551	to_be_resorted => ww
552	ELSE
553	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
554	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
555	ww_av = REAL( fill_value, KIND = wp )
556	ENDIF
557	to_be_resorted => ww_av
558	ENDIF
559	flag_nr = 3
560
561	IF ( mode == 'xy' ) grid = 'zw'
562
563	CASE DEFAULT
564	found = .FALSE.
565	grid = 'none'
566
567	END SELECT
568
569	IF ( found .AND. .NOT. resorted ) THEN
570	DO i = nxl, nxr
571	DO j = nys, nyn
572	DO k = nzb_do, nzt_do
573	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
574	REAL( fill_value, KIND = wp ), &
575	BTEST( wall_flags_0(k,j,i), flag_nr ) )
576	ENDDO
577	ENDDO
578	ENDDO
579	ENDIF
580
581	END SUBROUTINE doq_output_2d
582
583
584	!------------------------------------------------------------------------------!
585	!
586	! Description:
587	! ------------
588	!> Subroutine defining 3D output variables
589	!------------------------------------------------------------------------------!
590	SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, &
591	nzt_do )
592
593	IMPLICIT NONE
594
595	CHARACTER (LEN=*) :: variable !<
596
597	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
598	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
599	INTEGER(iwp) :: i !< index variable along x-direction
600	INTEGER(iwp) :: j !< index variable along y-direction
601	INTEGER(iwp) :: k !< index variable along z-direction
602	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
603	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
604
605	LOGICAL :: found !< true if variable is in list
606	LOGICAL :: resorted !< true if array is resorted
607
608	REAL(wp) :: fill_value !< value for the _FillValue attribute
609
610	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
611	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
612
613	flag_nr = 0
614	found = .TRUE.
615	resorted = .FALSE.
616
617	SELECT CASE ( TRIM( variable ) )
618
619	CASE ( 'ti' )
620	IF ( av == 0 ) THEN
621	to_be_resorted => ti
622	ELSE
623	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
624	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
625	ti_av = REAL( fill_value, KIND = wp )
626	ENDIF
627	to_be_resorted => ti_av
628	ENDIF
629	flag_nr = 0
630
631	CASE ( 'uu' )
632	IF ( av == 0 ) THEN
633	to_be_resorted => uu
634	ELSE
635	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
636	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
637	uu_av = REAL( fill_value, KIND = wp )
638	ENDIF
639	to_be_resorted => uu_av
640	ENDIF
641	flag_nr = 1
642
643	CASE ( 'vv' )
644	IF ( av == 0 ) THEN
645	to_be_resorted => vv
646	ELSE
647	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
648	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
649	vv_av = REAL( fill_value, KIND = wp )
650	ENDIF
651	to_be_resorted => vv_av
652	ENDIF
653	flag_nr = 2
654
655	CASE ( 'ww' )
656	IF ( av == 0 ) THEN
657	to_be_resorted => ww
658	ELSE
659	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
660	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
661	ww_av = REAL( fill_value, KIND = wp )
662	ENDIF
663	to_be_resorted => ww_av
664	ENDIF
665	flag_nr = 3
666
667	CASE DEFAULT
668	found = .FALSE.
669
670	END SELECT
671
672	IF ( found .AND. .NOT. resorted ) THEN
673	DO i = nxl, nxr
674	DO j = nys, nyn
675	DO k = nzb_do, nzt_do
676	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
677	REAL( fill_value, KIND = wp ), &
678	BTEST( wall_flags_0(k,j,i), flag_nr ) )
679	ENDDO
680	ENDDO
681	ENDDO
682	ENDIF
683
684	END SUBROUTINE doq_output_3d
685
686	! Description:
687	! ------------
688	!> Resorts the user-defined output quantity with indices (k,j,i) to a
689	!> temporary array with indices (i,j,k) for masked data output.
690	!------------------------------------------------------------------------------!
691	SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
692
693	USE control_parameters
694
695	USE indices
696
697	USE surface_mod, &
698	ONLY: get_topography_top_index_ji
699
700	IMPLICIT NONE
701
702	CHARACTER (LEN=*) :: variable !<
703	CHARACTER (LEN=5) :: grid !< flag to distinquish between staggered grids
704
705	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
706	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
707	INTEGER(iwp) :: i !< index variable along x-direction
708	INTEGER(iwp) :: j !< index variable along y-direction
709	INTEGER(iwp) :: k !< index variable along z-direction
710	INTEGER(iwp) :: mid !< masked output running index
711	INTEGER(iwp) :: topo_top_ind !< k index of highest horizontal surface
712
713	LOGICAL :: found !< true if variable is in list
714	LOGICAL :: resorted !< true if array is resorted
715
716	REAL(wp), &
717	DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
718	local_pf !<
719	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
720
721	flag_nr = 0
722	found = .TRUE.
723	resorted = .FALSE.
724	grid = 's'
725
726	SELECT CASE ( TRIM( variable ) )
727
728	CASE ( 'ti' )
729	IF ( av == 0 ) THEN
730	to_be_resorted => ti
731	ELSE
732	to_be_resorted => ti_av
733	ENDIF
734	grid = 's'
735	flag_nr = 0
736
737	CASE ( 'uu' )
738	IF ( av == 0 ) THEN
739	to_be_resorted => uu
740	ELSE
741	to_be_resorted => uu_av
742	ENDIF
743	grid = 'u'
744	flag_nr = 1
745
746	CASE ( 'vv' )
747	IF ( av == 0 ) THEN
748	to_be_resorted => vv
749	ELSE
750	to_be_resorted => vv_av
751	ENDIF
752	grid = 'v'
753	flag_nr = 2
754
755	CASE ( 'ww' )
756	IF ( av == 0 ) THEN
757	to_be_resorted => ww
758	ELSE
759	to_be_resorted => ww_av
760	ENDIF
761	grid = 'w'
762	flag_nr = 3
763
764
765	CASE DEFAULT
766	found = .FALSE.
767
768	END SELECT
769
770	IF ( found .AND. .NOT. resorted ) THEN
771	IF ( .NOT. mask_surface(mid) ) THEN
772	!
773	!-- Default masked output
774	DO i = 1, mask_size_l(mid,1)
775	DO j = 1, mask_size_l(mid,2)
776	DO k = 1, mask_size_l(mid,3)
777	local_pf(i,j,k) = to_be_resorted(mask_k(mid,k), &
778	mask_j(mid,j), &
779	mask_i(mid,i))
780	ENDDO
781	ENDDO
782	ENDDO
783
784	ELSE
785	!
786	!-- Terrain-following masked output
787	DO i = 1, mask_size_l(mid,1)
788	DO j = 1, mask_size_l(mid,2)
789	!
790	!-- Get k index of highest horizontal surface
791	topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), &
792	mask_i(mid,i), &
793	grid )
794	!
795	!-- Save output array
796	DO k = 1, mask_size_l(mid,3)
797	local_pf(i,j,k) = to_be_resorted( &
798	MIN( topo_top_ind+mask_k(mid,k), &
799	nzt+1 ), &
800	mask_j(mid,j), &
801	mask_i(mid,i) )
802	ENDDO
803	ENDDO
804	ENDDO
805
806	ENDIF
807	ENDIF
808
809	END SUBROUTINE doq_output_mask
810
811	!------------------------------------------------------------------------------!
812	! Description:
813	! ------------
814	!> Allocate required arrays
815	!------------------------------------------------------------------------------!
816	SUBROUTINE doq_init
817
818	IMPLICIT NONE
819
820	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
821
822	!
823	!-- Next line is to avoid compiler warnings about unused variables
824	IF ( timestep_number_at_prev_calc == 0 ) CONTINUE
825	!
826	!-- Preparatory steps and initialization of output arrays
827	IF ( .NOT. prepared_diagnostic_output_quantities ) CALL doq_prepare
828
829	initialized_diagnostic_output_quantities = .FALSE.
830
831	ivar = 1
832
833	DO WHILE ( ivar <= SIZE( do_all ) )
834
835	SELECT CASE ( TRIM( do_all(ivar) ) )
836	!
837	!-- Allocate array for 'turbulence intensity'
838	CASE ( 'ti' )
839	IF ( .NOT. ALLOCATED( ti ) ) THEN
840	ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
841	ti = 0.0_wp
842	ENDIF
843	!
844	!-- Allocate array for uu
845	CASE ( 'uu' )
846	IF ( .NOT. ALLOCATED( uu ) ) THEN
847	ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) )
848	uu = 0.0_wp
849	ENDIF
850
851	!
852	!-- Allocate array for vv
853	CASE ( 'vv' )
854	IF ( .NOT. ALLOCATED( vv ) ) THEN
855	ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) )
856	vv = 0.0_wp
857	ENDIF
858
859	!
860	!-- Allocate array for ww
861	CASE ( 'ww' )
862	IF ( .NOT. ALLOCATED( ww ) ) THEN
863	ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
864	ww = 0.0_wp
865	ENDIF
866
867	END SELECT
868
869	ivar = ivar + 1
870	ENDDO
871
872	initialized_diagnostic_output_quantities = .TRUE.
873
874	END SUBROUTINE doq_init
875
876
877	!--------------------------------------------------------------------------------------------------!
878	! Description:
879	! ------------
880	!> Calculate diagnostic quantities
881	!--------------------------------------------------------------------------------------------------!
882	SUBROUTINE doq_calculate
883
884	IMPLICIT NONE
885
886	INTEGER(iwp) :: i, j, k !< grid loop index in x-, y-, z-direction
887	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
888
889
890	! CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
891
892	!
893	!-- Save timestep number to check in time_integration if doq_calculate
894	!-- has been called already, since the CALL occurs at two locations, but the calculations need to be
895	!-- done only once per timestep.
896	timestep_number_at_prev_calc = current_timestep_number
897
898	ivar = 1
899
900	DO WHILE ( ivar <= SIZE( do_all ) )
901
902	SELECT CASE ( TRIM( do_all(ivar) ) )
903	!
904	!-- Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
905	CASE ( 'ti' )
906	DO i = nxl, nxr
907	DO j = nys, nyn
908	DO k = nzb+1, nzt
909	ti(k,j,i) = 0.25_wp * SQRT( &
910	( ( w(k,j+1,i) + w(k-1,j+1,i) &
911	- w(k,j-1,i) - w(k-1,j-1,i) ) * ddy &
912	- ( v(k+1,j,i) + v(k+1,j+1,i) &
913	- v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2 &
914	+ ( ( u(k+1,j,i) + u(k+1,j,i+1) &
915	- u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k) &
916	- ( w(k,j,i+1) + w(k-1,j,i+1) &
917	- w(k,j,i-1) - w(k-1,j,i-1) ) * ddx )**2 &
918	+ ( ( v(k,j,i+1) + v(k,j+1,i+1) &
919	- v(k,j,i-1) - v(k,j+1,i-1) ) * ddx &
920	- ( u(k,j+1,i) + u(k,j+1,i+1) &
921	- u(k,j-1,i) - u(k,j-1,i+1) ) * ddy )**2 ) &
922	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) )
923	ENDDO
924	ENDDO
925	ENDDO
926	!
927	!-- uu
928	CASE ( 'uu' )
929	DO i = nxl, nxr
930	DO j = nys, nyn
931	DO k = nzb+1, nzt
932	uu(k,j,i) = u(k,j,i) * u(k,j,i) &
933	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1) )
934	ENDDO
935	ENDDO
936	ENDDO
937	!
938	!-- vv
939	CASE ( 'vv' )
940	DO i = nxl, nxr
941	DO j = nys, nyn
942	DO k = nzb+1, nzt
943	vv(k,j,i) = v(k,j,i) * v(k,j,i) &
944	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2) )
945	ENDDO
946	ENDDO
947	ENDDO
948	!
949	!-- ww
950	CASE ( 'ww' )
951	DO i = nxl, nxr
952	DO j = nys, nyn
953	DO k = nzb+1, nzt-1
954	ww(k,j,i) = w(k,j,i) * w(k,j,i) &
955	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3) )
956	ENDDO
957	ENDDO
958	ENDDO
959
960	END SELECT
961
962	ivar = ivar + 1
963	ENDDO
964
965	! CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
966
967	END SUBROUTINE doq_calculate
968
969
970	!------------------------------------------------------------------------------!
971	! Description:
972	! ------------
973	!> ...
974	!------------------------------------------------------------------------------!
975	SUBROUTINE doq_prepare
976
977
978	USE control_parameters, &
979	ONLY: do2d, do3d, domask, masks
980
981	IMPLICIT NONE
982
983	CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) :: do2d_var = ' ' !<
984	!< label array for 2d output quantities
985
986	INTEGER(iwp) :: av !< index defining type of output, av=0 instantaneous, av=1 averaged
987	INTEGER(iwp) :: ivar !< loop index
988	INTEGER(iwp) :: ivar_all !< loop index
989	INTEGER(iwp) :: l !< index for cutting string
990	INTEGER(iwp) :: mid !< masked output running index
991
992	prepared_diagnostic_output_quantities = .FALSE.
993
994	ivar = 1
995	ivar_all = 1
996
997	DO av = 0, 1
998	!
999	!-- Remove _xy, _xz, or _yz from string
1000	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1001	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1002	!
1003	!-- Gather 2d output quantity names.
1004	!-- Check for double occurrence of output quantity, e.g. by _xy,
1005	!-- _yz, _xz.
1006	DO WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
1007	IF ( .NOT. ANY( do_all == do2d_var(av,ivar) ) ) THEN
1008	do_all(ivar_all) = do2d_var(av,ivar)
1009	ENDIF
1010	ivar = ivar + 1
1011	ivar_all = ivar_all + 1
1012	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1013	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1014	ENDDO
1015
1016	ivar = 1
1017	!
1018	!-- Gather 3d output quantity names
1019	DO WHILE ( do3d(av,ivar)(1:1) /= ' ' )
1020	do_all(ivar_all) = do3d(av,ivar)
1021	ivar = ivar + 1
1022	ivar_all = ivar_all + 1
1023	ENDDO
1024
1025	ivar = 1
1026	!
1027	!-- Gather masked output quantity names. Also check for double output
1028	!-- e.g. by different masks.
1029	DO mid = 1, masks
1030	DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
1031	IF ( .NOT. ANY( do_all == domask(mid,av,ivar) ) ) THEN
1032	do_all(ivar_all) = domask(mid,av,ivar)
1033	ENDIF
1034
1035	ivar = ivar + 1
1036	ivar_all = ivar_all + 1
1037	ENDDO
1038	ivar = 1
1039	ENDDO
1040
1041	ENDDO
1042
1043	prepared_diagnostic_output_quantities = .TRUE.
1044
1045	END SUBROUTINE doq_prepare
1046
1047	!------------------------------------------------------------------------------!
1048	! Description:
1049	! ------------
1050	!> Subroutine reads local (subdomain) restart data
1051	!> Note: With the current structure reading of non-standard array is not
1052	!> possible
1053	!------------------------------------------------------------------------------!
1054	! SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
1055	! nxr_on_file, nynf, nync, nyn_on_file, nysf, &
1056	! nysc, nys_on_file, tmp_3d_non_standard, found )
1057	!
1058	!
1059	! USE control_parameters
1060	!
1061	! USE indices
1062	!
1063	! USE kinds
1064	!
1065	! USE pegrid
1066	!
1067	!
1068	! IMPLICIT NONE
1069	!
1070	! INTEGER(iwp) :: k !<
1071	! INTEGER(iwp) :: nxlc !<
1072	! INTEGER(iwp) :: nxlf !<
1073	! INTEGER(iwp) :: nxl_on_file !<
1074	! INTEGER(iwp) :: nxrc !<
1075	! INTEGER(iwp) :: nxrf !<
1076	! INTEGER(iwp) :: nxr_on_file !<
1077	! INTEGER(iwp) :: nync !<
1078	! INTEGER(iwp) :: nynf !<
1079	! INTEGER(iwp) :: nyn_on_file !<
1080	! INTEGER(iwp) :: nysc !<
1081	! INTEGER(iwp) :: nysf !<
1082	! INTEGER(iwp) :: nys_on_file !<
1083	!
1084	! LOGICAL, INTENT(OUT) :: found
1085	!
1086	! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions
1087	! !
1088	! !-- If temporary non-standard array for reading is already allocated,
1089	! !-- deallocate it.
1090	! IF ( ALLOCATED( tmp_3d_non_standard ) ) DEALLOCATE( tmp_3d_non_standard )
1091	!
1092	! found = .TRUE.
1093	!
1094	! SELECT CASE ( restart_string(1:length) )
1095	!
1096	! CASE ( 'ti_av' )
1097	! IF ( .NOT. ALLOCATED( ti_av ) ) THEN
1098	! ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1099	! ENDIF
1100	! IF ( k == 1 ) THEN
1101	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
1102	! nxl_on_file:nxr_on_file) )
1103	! READ ( 13 ) tmp_3d_non_standard
1104	! ENDIF
1105	! ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1106	!
1107	! CASE ( 'uu_av' )
1108	! IF ( .NOT. ALLOCATED( uu_av ) ) THEN
1109	! ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1110	! ENDIF
1111	! IF ( k == 1 ) THEN
1112	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
1113	! nxl_on_file:nxr_on_file) )
1114	! READ ( 13 ) tmp_3d_non_standard
1115	! ENDIF
1116	! uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1117	!
1118	! CASE ( 'vv_av' )
1119	! IF ( .NOT. ALLOCATED( vv_av ) ) THEN
1120	! ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1121	! ENDIF
1122	! IF ( k == 1 ) THEN
1123	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
1124	! nxl_on_file:nxr_on_file) )
1125	! READ ( 13 ) tmp_3d_non_standard
1126	! ENDIF
1127	! vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1128	!
1129	! CASE ( 'ww_av' )
1130	! IF ( .NOT. ALLOCATED( ww_av ) ) THEN
1131	! ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1132	! ENDIF
1133	! IF ( k == 1 ) THEN
1134	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
1135	! nxl_on_file:nxr_on_file) )
1136	! READ ( 13 ) tmp_3d_non_standard
1137	! ENDIF
1138	! ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1139	!
1140	!
1141	! CASE DEFAULT
1142	!
1143	! found = .FALSE.
1144	!
1145	! END SELECT
1146	!
1147	! END SUBROUTINE doq_rrd_local
1148
1149	!------------------------------------------------------------------------------!
1150	! Description:
1151	! ------------
1152	!> Subroutine writes local (subdomain) restart data
1153	!------------------------------------------------------------------------------!
1154	SUBROUTINE doq_wrd_local
1155
1156
1157	IMPLICIT NONE
1158
1159	IF ( ALLOCATED( ti_av ) ) THEN
1160	CALL wrd_write_string( 'ti_av' )
1161	WRITE ( 14 ) ti_av
1162	ENDIF
1163
1164	IF ( ALLOCATED( uu_av ) ) THEN
1165	CALL wrd_write_string( 'uu_av' )
1166	WRITE ( 14 ) uu_av
1167	ENDIF
1168
1169	IF ( ALLOCATED( vv_av ) ) THEN
1170	CALL wrd_write_string( 'vv_av' )
1171	WRITE ( 14 ) vv_av
1172	ENDIF
1173
1174	IF ( ALLOCATED( ww_av ) ) THEN
1175	CALL wrd_write_string( 'ww_av' )
1176	WRITE ( 14 ) ww_av
1177	ENDIF
1178
1179
1180	END SUBROUTINE doq_wrd_local
1181
1182
1183
1184	END MODULE diagnostic_output_quantities_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |