Home

Context Navigation

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

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