Home

Context Navigation

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

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