Home

Context Navigation

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

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