Home

Context Navigation

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

Last change on this file since 4343 was 4331, checked in by suehring, 5 years ago
New diagnostic output for 10-m wind speed; Diagnostic output of 2-m potential temperature moved to diagnostic output
Property svn:keywords set to `Id`
File size: 50.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 4331 2019-12-10 18:25:02Z oliver.maas $
27	! - Modularize 2-m potential temperature output
28	! - New output for 10-m wind speed
29	!
30	! 4329 2019-12-10 15:46:36Z motisi
31	! Renamed wall_flags_0 to wall_flags_static_0
32	!
33	! 4182 2019-08-22 15:20:23Z scharf
34	! Corrected "Former revisions" section
35	!
36	! 4167 2019-08-16 11:01:48Z suehring
37	! Changed behaviour of masked output over surface to follow terrain and ignore
38	! buildings (J.Resler, T.Gronemeier)
39	!
40	! 4157 2019-08-14 09:19:12Z suehring
41	! Initialization restructured, in order to work also when data output during
42	! spin-up is enabled.
43	!
44	! 4132 2019-08-02 12:34:17Z suehring
45	! Bugfix in masked data output
46	!
47	! 4069 2019-07-01 14:05:51Z Giersch
48	! Masked output running index mid has been introduced as a local variable to
49	! avoid runtime error (Loop variable has been modified) in time_integration
50	!
51	! 4039 2019-06-18 10:32:41Z suehring
52	! - Add output of uu, vv, ww to enable variance calculation according temporal
53	! EC method
54	! - Allocate arrays only when they are required
55	! - Formatting adjustment
56	! - Rename subroutines
57	! - Further modularization
58	!
59	! 3998 2019-05-23 13:38:11Z suehring
60	! Bugfix in gathering all output strings
61	!
62	! 3995 2019-05-22 18:59:54Z suehring
63	! Avoid compiler warnings about unused variable and fix string operation which
64	! is not allowed with PGI compiler
65	!
66	! 3994 2019-05-22 18:08:09Z suehring
67	! Initial revision
68	!
69	! Authors:
70	! --------
71	! @author Farah Kanani-Suehring
72	!
73	!
74	! Description:
75	! ------------
76	!> ...
77	!------------------------------------------------------------------------------!
78	MODULE diagnostic_output_quantities_mod
79
80	USE arrays_3d, &
81	ONLY: ddzu, &
82	pt, &
83	u, &
84	v, &
85	w, &
86	zu, &
87	zw
88
89	USE basic_constants_and_equations_mod, &
90	ONLY: kappa
91
92	USE control_parameters, &
93	ONLY: current_timestep_number, &
94	data_output, &
95	message_string, &
96	varnamelength
97	!
98	! USE cpulog, &
99	! ONLY: cpu_log, log_point
100
101	USE grid_variables, &
102	ONLY: ddx, ddy
103
104	USE indices, &
105	ONLY: nbgp, &
106	nxl, &
107	nxlg, &
108	nxr, &
109	nxrg, &
110	nyn, &
111	nyng, &
112	nys, &
113	nysg, &
114	nzb, &
115	nzt, &
116	wall_flags_static_0
117
118	USE kinds
119
120	USE surface_mod, &
121	ONLY: surf_def_h, &
122	surf_lsm_h, &
123	surf_type, &
124	surf_usm_h
125
126
127	IMPLICIT NONE
128
129	CHARACTER(LEN=varnamelength), DIMENSION(500) :: do_all = ' '
130
131	INTEGER(iwp) :: timestep_number_at_prev_calc = 0 !< ...at previous diagnostic output calculation
132
133	LOGICAL :: initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
134	LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is p
135
136	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m !< 2-m air potential temperature
137	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m_av !< averaged 2-m air potential temperature
138	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m !< horizontal wind speed at 10m
139	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m_av !< averaged horizontal wind speed at 10m
140
141	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity
142	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity
143	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu !< uu
144	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu
145	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv !< vv
146	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv
147	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww !< ww
148	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww_av !< mean of ww
149
150
151	SAVE
152
153	PRIVATE
154
155	!
156	!-- Public variables
157	PUBLIC do_all, &
158	initialized_diagnostic_output_quantities, &
159	prepared_diagnostic_output_quantities, &
160	timestep_number_at_prev_calc, &
161	pt_2m_av, &
162	ti_av, &
163	uu_av, &
164	uv_10m_av, &
165	vv_av, &
166	ww_av
167	!
168	!-- Public routines
169	PUBLIC doq_3d_data_averaging, &
170	doq_calculate, &
171	doq_check_data_output, &
172	doq_define_netcdf_grid, &
173	doq_output_2d, &
174	doq_output_3d, &
175	doq_output_mask, &
176	doq_init, &
177	doq_wrd_local
178	! doq_rrd_local, &
179
180
181	INTERFACE doq_3d_data_averaging
182	MODULE PROCEDURE doq_3d_data_averaging
183	END INTERFACE doq_3d_data_averaging
184
185	INTERFACE doq_calculate
186	MODULE PROCEDURE doq_calculate
187	END INTERFACE doq_calculate
188
189	INTERFACE doq_check_data_output
190	MODULE PROCEDURE doq_check_data_output
191	END INTERFACE doq_check_data_output
192
193	INTERFACE doq_define_netcdf_grid
194	MODULE PROCEDURE doq_define_netcdf_grid
195	END INTERFACE doq_define_netcdf_grid
196
197	INTERFACE doq_output_2d
198	MODULE PROCEDURE doq_output_2d
199	END INTERFACE doq_output_2d
200
201	INTERFACE doq_output_3d
202	MODULE PROCEDURE doq_output_3d
203	END INTERFACE doq_output_3d
204
205	INTERFACE doq_output_mask
206	MODULE PROCEDURE doq_output_mask
207	END INTERFACE doq_output_mask
208
209	INTERFACE doq_init
210	MODULE PROCEDURE doq_init
211	END INTERFACE doq_init
212
213	INTERFACE doq_prepare
214	MODULE PROCEDURE doq_prepare
215	END INTERFACE doq_prepare
216
217	! INTERFACE doq_rrd_local
218	! MODULE PROCEDURE doq_rrd_local
219	! END INTERFACE doq_rrd_local
220
221	INTERFACE doq_wrd_local
222	MODULE PROCEDURE doq_wrd_local
223	END INTERFACE doq_wrd_local
224
225
226	CONTAINS
227
228	!------------------------------------------------------------------------------!
229	! Description:
230	! ------------
231	!> Sum up and time-average diagnostic output quantities as well as allocate
232	!> the array necessary for storing the average.
233	!------------------------------------------------------------------------------!
234	SUBROUTINE doq_3d_data_averaging( mode, variable )
235
236	USE control_parameters, &
237	ONLY: average_count_3d
238
239	CHARACTER (LEN=*) :: mode !<
240	CHARACTER (LEN=*) :: variable !<
241
242	INTEGER(iwp) :: i !<
243	INTEGER(iwp) :: j !<
244	INTEGER(iwp) :: k !<
245
246	IF ( mode == 'allocate' ) THEN
247
248	SELECT CASE ( TRIM( variable ) )
249
250	CASE ( 'ti' )
251	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
252	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
253	ENDIF
254	ti_av = 0.0_wp
255
256	CASE ( 'uu' )
257	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
258	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
259	ENDIF
260	uu_av = 0.0_wp
261
262	CASE ( 'vv' )
263	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
264	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
265	ENDIF
266	vv_av = 0.0_wp
267
268	CASE ( 'ww' )
269	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
270	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
271	ENDIF
272	ww_av = 0.0_wp
273
274	CASE ( 'theta_2m*' )
275	IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
276	ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
277	ENDIF
278	pt_2m_av = 0.0_wp
279
280	CASE ( 'wspeed_10m*' )
281	IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
282	ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
283	ENDIF
284	uv_10m_av = 0.0_wp
285
286	CASE DEFAULT
287	CONTINUE
288
289	END SELECT
290
291	ELSEIF ( mode == 'sum' ) THEN
292
293	SELECT CASE ( TRIM( variable ) )
294
295	CASE ( 'ti' )
296	IF ( ALLOCATED( ti_av ) ) THEN
297	DO i = nxl, nxr
298	DO j = nys, nyn
299	DO k = nzb, nzt+1
300	ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
301	ENDDO
302	ENDDO
303	ENDDO
304	ENDIF
305
306	CASE ( 'uu' )
307	IF ( ALLOCATED( uu_av ) ) THEN
308	DO i = nxl, nxr
309	DO j = nys, nyn
310	DO k = nzb, nzt+1
311	uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
312	ENDDO
313	ENDDO
314	ENDDO
315	ENDIF
316
317	CASE ( 'vv' )
318	IF ( ALLOCATED( vv_av ) ) THEN
319	DO i = nxl, nxr
320	DO j = nys, nyn
321	DO k = nzb, nzt+1
322	vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
323	ENDDO
324	ENDDO
325	ENDDO
326	ENDIF
327
328	CASE ( 'ww' )
329	IF ( ALLOCATED( ww_av ) ) THEN
330	DO i = nxl, nxr
331	DO j = nys, nyn
332	DO k = nzb, nzt+1
333	ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
334	ENDDO
335	ENDDO
336	ENDDO
337	ENDIF
338
339	CASE ( 'theta_2m*' )
340	IF ( ALLOCATED( pt_2m_av ) ) THEN
341	DO i = nxl, nxr
342	DO j = nys, nyn
343	pt_2m_av(j,i) = pt_2m_av(j,i) + pt_2m(j,i)
344	ENDDO
345	ENDDO
346	ENDIF
347
348	CASE ( 'wspeed_10m*' )
349	IF ( ALLOCATED( uv_10m_av ) ) THEN
350	DO i = nxl, nxr
351	DO j = nys, nyn
352	uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i)
353	ENDDO
354	ENDDO
355	ENDIF
356
357	CASE DEFAULT
358	CONTINUE
359
360	END SELECT
361
362	ELSEIF ( mode == 'average' ) THEN
363
364	SELECT CASE ( TRIM( variable ) )
365
366	CASE ( 'ti' )
367	IF ( ALLOCATED( ti_av ) ) THEN
368	DO i = nxl, nxr
369	DO j = nys, nyn
370	DO k = nzb, nzt+1
371	ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
372	ENDDO
373	ENDDO
374	ENDDO
375	ENDIF
376
377	CASE ( 'uu' )
378	IF ( ALLOCATED( uu_av ) ) THEN
379	DO i = nxl, nxr
380	DO j = nys, nyn
381	DO k = nzb, nzt+1
382	uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
383	ENDDO
384	ENDDO
385	ENDDO
386	ENDIF
387
388	CASE ( 'vv' )
389	IF ( ALLOCATED( vv_av ) ) THEN
390	DO i = nxl, nxr
391	DO j = nys, nyn
392	DO k = nzb, nzt+1
393	vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
394	ENDDO
395	ENDDO
396	ENDDO
397	ENDIF
398
399	CASE ( 'ww' )
400	IF ( ALLOCATED( ww_av ) ) THEN
401	DO i = nxl, nxr
402	DO j = nys, nyn
403	DO k = nzb, nzt+1
404	ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
405	ENDDO
406	ENDDO
407	ENDDO
408	ENDIF
409
410	CASE ( 'theta_2m*' )
411	IF ( ALLOCATED( pt_2m_av ) ) THEN
412	DO i = nxlg, nxrg
413	DO j = nysg, nyng
414	pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
415	ENDDO
416	ENDDO
417	CALL exchange_horiz_2d( pt_2m_av, nbgp )
418	ENDIF
419
420	CASE ( 'wspeed_10m*' )
421	IF ( ALLOCATED( uv_10m_av ) ) THEN
422	DO i = nxlg, nxrg
423	DO j = nysg, nyng
424	uv_10m_av(j,i) = uv_10m_av(j,i) / REAL( average_count_3d, KIND=wp )
425	ENDDO
426	ENDDO
427	CALL exchange_horiz_2d( uv_10m_av, nbgp )
428	ENDIF
429
430	END SELECT
431
432	ENDIF
433
434
435	END SUBROUTINE doq_3d_data_averaging
436
437	!------------------------------------------------------------------------------!
438	! Description:
439	! ------------
440	!> Check data output for diagnostic output
441	!------------------------------------------------------------------------------!
442	SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
443
444	IMPLICIT NONE
445
446	CHARACTER (LEN=*) :: unit !<
447	CHARACTER (LEN=*) :: var !<
448
449	INTEGER(iwp), OPTIONAL, INTENT(IN) :: i !< Current element of data_output
450	INTEGER(iwp), OPTIONAL, INTENT(IN) :: ilen !< Length of current entry in data_output
451	INTEGER(iwp), OPTIONAL, INTENT(IN) :: k !< Output is xy mode? 0 = no, 1 = yes
452
453	SELECT CASE ( TRIM( var ) )
454
455	CASE ( 'ti' )
456	unit = '1/s'
457
458	CASE ( 'uu' )
459	unit = 'm2/s2'
460
461	CASE ( 'vv' )
462	unit = 'm2/s2'
463
464	CASE ( 'ww' )
465	unit = 'm2/s2'
466	!
467	!-- Treat horizotal cross-section output quanatities
468	CASE ( 'theta_2m', 'wspeed_10m' )
469	!
470	!-- Check if output quantity is _xy only.
471	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
472	message_string = 'illegal value for data_output: "' // &
473	TRIM( var ) // '" & only 2d-horizontal ' // &
474	'cross sections are allowed for this value'
475	CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
476	ENDIF
477
478	IF ( TRIM( var ) == 'theta_2m*' ) unit = 'K'
479	IF ( TRIM( var ) == 'wspeed_10m*' ) unit = 'm/s'
480
481	CASE DEFAULT
482	unit = 'illegal'
483
484	END SELECT
485
486
487	END SUBROUTINE doq_check_data_output
488
489	!------------------------------------------------------------------------------!
490	!
491	! Description:
492	! ------------
493	!> Subroutine defining appropriate grid for netcdf variables.
494	!------------------------------------------------------------------------------!
495	SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
496
497	IMPLICIT NONE
498
499	CHARACTER (LEN=*), INTENT(IN) :: variable !<
500	LOGICAL, INTENT(OUT) :: found !<
501	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
502	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
503	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
504
505	found = .TRUE.
506
507	SELECT CASE ( TRIM( variable ) )
508	!
509	!-- s grid
510	CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz' )
511
512	grid_x = 'x'
513	grid_y = 'y'
514	grid_z = 'zu'
515	!
516	!-- s grid surface variables
517	CASE ( 'theta_2m_xy', 'wspeed_10m' )
518
519	grid_x = 'x'
520	grid_y = 'y'
521	grid_z = 'zu'
522	!
523	!-- u grid
524	CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
525
526	grid_x = 'xu'
527	grid_y = 'y'
528	grid_z = 'zu'
529	!
530	!-- v grid
531	CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz' )
532
533	grid_x = 'x'
534	grid_y = 'yv'
535	grid_z = 'zu'
536	!
537	!-- w grid
538	CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz' )
539
540	grid_x = 'x'
541	grid_y = 'y'
542	grid_z = 'zw'
543
544	CASE DEFAULT
545	found = .FALSE.
546	grid_x = 'none'
547	grid_y = 'none'
548	grid_z = 'none'
549
550	END SELECT
551
552
553	END SUBROUTINE doq_define_netcdf_grid
554
555	!------------------------------------------------------------------------------!
556	!
557	! Description:
558	! ------------
559	!> Subroutine defining 2D output variables
560	!------------------------------------------------------------------------------!
561	SUBROUTINE doq_output_2d( av, variable, found, grid, &
562	mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
563
564
565	IMPLICIT NONE
566
567	CHARACTER (LEN=*) :: grid !<
568	CHARACTER (LEN=*) :: mode !<
569	CHARACTER (LEN=*) :: variable !<
570
571	INTEGER(iwp) :: av !< value indicating averaged or non-averaged output
572	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
573	INTEGER(iwp) :: i !< grid index x-direction
574	INTEGER(iwp) :: j !< grid index y-direction
575	INTEGER(iwp) :: k !< grid index z-direction
576	INTEGER(iwp) :: nzb_do !<
577	INTEGER(iwp) :: nzt_do !<
578
579	LOGICAL :: found !< true if variable is in list
580	LOGICAL :: resorted !< true if array is resorted
581	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
582
583	REAL(wp) :: fill_value !< value for the _FillValue attribute
584
585	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
586	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
587
588	flag_nr = 0
589	found = .TRUE.
590	resorted = .FALSE.
591	two_d = .FALSE.
592
593	SELECT CASE ( TRIM( variable ) )
594
595	CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
596	IF ( av == 0 ) THEN
597	to_be_resorted => ti
598	ELSE
599	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
600	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
601	ti_av = REAL( fill_value, KIND = wp )
602	ENDIF
603	to_be_resorted => ti_av
604	ENDIF
605	flag_nr = 0
606
607	IF ( mode == 'xy' ) grid = 'zu'
608
609	CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
610	IF ( av == 0 ) THEN
611	to_be_resorted => uu
612	ELSE
613	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
614	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
615	uu_av = REAL( fill_value, KIND = wp )
616	ENDIF
617	to_be_resorted => uu_av
618	ENDIF
619	flag_nr = 1
620
621	IF ( mode == 'xy' ) grid = 'zu'
622
623	CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
624	IF ( av == 0 ) THEN
625	to_be_resorted => vv
626	ELSE
627	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
628	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
629	vv_av = REAL( fill_value, KIND = wp )
630	ENDIF
631	to_be_resorted => vv_av
632	ENDIF
633	flag_nr = 2
634
635	IF ( mode == 'xy' ) grid = 'zu'
636
637	CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
638	IF ( av == 0 ) THEN
639	to_be_resorted => ww
640	ELSE
641	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
642	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
643	ww_av = REAL( fill_value, KIND = wp )
644	ENDIF
645	to_be_resorted => ww_av
646	ENDIF
647	flag_nr = 3
648
649	IF ( mode == 'xy' ) grid = 'zw'
650
651	CASE ( 'theta_2m*_xy' ) ! 2d-array
652	IF ( av == 0 ) THEN
653	DO i = nxl, nxr
654	DO j = nys, nyn
655	local_pf(i,j,nzb+1) = pt_2m(j,i)
656	ENDDO
657	ENDDO
658	ELSE
659	IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
660	ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
661	pt_2m_av = REAL( fill_value, KIND = wp )
662	ENDIF
663	DO i = nxl, nxr
664	DO j = nys, nyn
665	local_pf(i,j,nzb+1) = pt_2m_av(j,i)
666	ENDDO
667	ENDDO
668	ENDIF
669	resorted = .TRUE.
670	two_d = .TRUE.
671	grid = 'zu1'
672
673	CASE ( 'wspeed_10m*_xy' ) ! 2d-array
674	IF ( av == 0 ) THEN
675	DO i = nxl, nxr
676	DO j = nys, nyn
677	local_pf(i,j,nzb+1) = uv_10m(j,i)
678	ENDDO
679	ENDDO
680	ELSE
681	IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
682	ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
683	uv_10m_av = REAL( fill_value, KIND = wp )
684	ENDIF
685	DO i = nxl, nxr
686	DO j = nys, nyn
687	local_pf(i,j,nzb+1) = uv_10m_av(j,i)
688	ENDDO
689	ENDDO
690	ENDIF
691	resorted = .TRUE.
692	two_d = .TRUE.
693	grid = 'zu1'
694
695	CASE DEFAULT
696	found = .FALSE.
697	grid = 'none'
698
699	END SELECT
700
701	IF ( found .AND. .NOT. resorted ) THEN
702	DO i = nxl, nxr
703	DO j = nys, nyn
704	DO k = nzb_do, nzt_do
705	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
706	REAL( fill_value, KIND = wp ), &
707	BTEST( wall_flags_static_0(k,j,i), flag_nr ) )
708	ENDDO
709	ENDDO
710	ENDDO
711	ENDIF
712
713	END SUBROUTINE doq_output_2d
714
715
716	!------------------------------------------------------------------------------!
717	!
718	! Description:
719	! ------------
720	!> Subroutine defining 3D output variables
721	!------------------------------------------------------------------------------!
722	SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, &
723	nzt_do )
724
725	IMPLICIT NONE
726
727	CHARACTER (LEN=*) :: variable !<
728
729	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
730	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
731	INTEGER(iwp) :: i !< index variable along x-direction
732	INTEGER(iwp) :: j !< index variable along y-direction
733	INTEGER(iwp) :: k !< index variable along z-direction
734	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
735	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
736
737	LOGICAL :: found !< true if variable is in list
738	LOGICAL :: resorted !< true if array is resorted
739
740	REAL(wp) :: fill_value !< value for the _FillValue attribute
741
742	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
743	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
744
745	flag_nr = 0
746	found = .TRUE.
747	resorted = .FALSE.
748
749	SELECT CASE ( TRIM( variable ) )
750
751	CASE ( 'ti' )
752	IF ( av == 0 ) THEN
753	to_be_resorted => ti
754	ELSE
755	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
756	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
757	ti_av = REAL( fill_value, KIND = wp )
758	ENDIF
759	to_be_resorted => ti_av
760	ENDIF
761	flag_nr = 0
762
763	CASE ( 'uu' )
764	IF ( av == 0 ) THEN
765	to_be_resorted => uu
766	ELSE
767	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
768	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
769	uu_av = REAL( fill_value, KIND = wp )
770	ENDIF
771	to_be_resorted => uu_av
772	ENDIF
773	flag_nr = 1
774
775	CASE ( 'vv' )
776	IF ( av == 0 ) THEN
777	to_be_resorted => vv
778	ELSE
779	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
780	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
781	vv_av = REAL( fill_value, KIND = wp )
782	ENDIF
783	to_be_resorted => vv_av
784	ENDIF
785	flag_nr = 2
786
787	CASE ( 'ww' )
788	IF ( av == 0 ) THEN
789	to_be_resorted => ww
790	ELSE
791	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
792	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
793	ww_av = REAL( fill_value, KIND = wp )
794	ENDIF
795	to_be_resorted => ww_av
796	ENDIF
797	flag_nr = 3
798
799	CASE DEFAULT
800	found = .FALSE.
801
802	END SELECT
803
804	IF ( found .AND. .NOT. resorted ) THEN
805	DO i = nxl, nxr
806	DO j = nys, nyn
807	DO k = nzb_do, nzt_do
808	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
809	REAL( fill_value, KIND = wp ), &
810	BTEST( wall_flags_static_0(k,j,i), flag_nr ) )
811	ENDDO
812	ENDDO
813	ENDDO
814	ENDIF
815
816	END SUBROUTINE doq_output_3d
817
818	! Description:
819	! ------------
820	!> Resorts the user-defined output quantity with indices (k,j,i) to a
821	!> temporary array with indices (i,j,k) for masked data output.
822	!------------------------------------------------------------------------------!
823	SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
824
825	USE control_parameters
826
827	USE indices
828
829	IMPLICIT NONE
830
831	CHARACTER (LEN=*) :: variable !<
832	CHARACTER (LEN=5) :: grid !< flag to distinquish between staggered grids
833
834	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
835	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
836	INTEGER(iwp) :: i !< index variable along x-direction
837	INTEGER(iwp) :: j !< index variable along y-direction
838	INTEGER(iwp) :: k !< index variable along z-direction
839	INTEGER(iwp) :: im !< loop index for masked variables
840	INTEGER(iwp) :: jm !< loop index for masked variables
841	INTEGER(iwp) :: kk !< masked output index variable along z-direction
842	INTEGER(iwp) :: mid !< masked output running index
843	INTEGER(iwp) :: ktt !< k index of highest horizontal surface
844
845	LOGICAL :: found !< true if variable is in list
846	LOGICAL :: resorted !< true if array is resorted
847
848	REAL(wp), &
849	DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
850	local_pf !<
851	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
852
853	REAL(wp), PARAMETER :: fill_value = -9999.0_wp !< value for the _FillValue attribute
854
855	flag_nr = 0
856	found = .TRUE.
857	resorted = .FALSE.
858	grid = 's'
859
860	SELECT CASE ( TRIM( variable ) )
861
862	CASE ( 'ti' )
863	IF ( av == 0 ) THEN
864	to_be_resorted => ti
865	ELSE
866	to_be_resorted => ti_av
867	ENDIF
868	grid = 's'
869	flag_nr = 0
870
871	CASE ( 'uu' )
872	IF ( av == 0 ) THEN
873	to_be_resorted => uu
874	ELSE
875	to_be_resorted => uu_av
876	ENDIF
877	grid = 'u'
878	flag_nr = 1
879
880	CASE ( 'vv' )
881	IF ( av == 0 ) THEN
882	to_be_resorted => vv
883	ELSE
884	to_be_resorted => vv_av
885	ENDIF
886	grid = 'v'
887	flag_nr = 2
888
889	CASE ( 'ww' )
890	IF ( av == 0 ) THEN
891	to_be_resorted => ww
892	ELSE
893	to_be_resorted => ww_av
894	ENDIF
895	grid = 'w'
896	flag_nr = 3
897
898	CASE DEFAULT
899	found = .FALSE.
900
901	END SELECT
902
903	IF ( found .AND. .NOT. resorted ) THEN
904	IF ( .NOT. mask_surface(mid) ) THEN
905	!
906	!-- Default masked output
907	DO i = 1, mask_size_l(mid,1)
908	DO j = 1, mask_size_l(mid,2)
909	DO k = 1, mask_size_l(mid,3)
910	local_pf(i,j,k) = to_be_resorted(mask_k(mid,k), &
911	mask_j(mid,j), &
912	mask_i(mid,i))
913	ENDDO
914	ENDDO
915	ENDDO
916
917	ELSE
918	!
919	!-- Terrain-following masked output
920	DO i = 1, mask_size_l(mid,1)
921	DO j = 1, mask_size_l(mid,2)
922	!
923	!-- Get k index of the highest terraing surface
924	im = mask_i(mid,i)
925	jm = mask_j(mid,j)
926	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_static_0(:,jm,im), 5 )), &
927	DIM = 1 ) - 1
928	DO k = 1, mask_size_l(mid,3)
929	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
930	!
931	!-- Set value if not in building
932	IF ( BTEST( wall_flags_static_0(kk,jm,im), 6 ) ) THEN
933	local_pf(i,j,k) = fill_value
934	ELSE
935	local_pf(i,j,k) = to_be_resorted(kk,jm,im)
936	ENDIF
937	ENDDO
938	ENDDO
939	ENDDO
940
941	ENDIF
942	ENDIF
943
944	END SUBROUTINE doq_output_mask
945
946	!------------------------------------------------------------------------------!
947	! Description:
948	! ------------
949	!> Allocate required arrays
950	!------------------------------------------------------------------------------!
951	SUBROUTINE doq_init
952
953	IMPLICIT NONE
954
955	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
956
957	!
958	!-- Next line is to avoid compiler warnings about unused variables
959	IF ( timestep_number_at_prev_calc == 0 ) CONTINUE
960	!
961	!-- Preparatory steps and initialization of output arrays
962	IF ( .NOT. prepared_diagnostic_output_quantities ) CALL doq_prepare
963
964	initialized_diagnostic_output_quantities = .FALSE.
965
966	ivar = 1
967
968	DO WHILE ( ivar <= SIZE( do_all ) )
969
970	SELECT CASE ( TRIM( do_all(ivar) ) )
971	!
972	!-- Allocate array for 'turbulence intensity'
973	CASE ( 'ti' )
974	IF ( .NOT. ALLOCATED( ti ) ) THEN
975	ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
976	ti = 0.0_wp
977	ENDIF
978	!
979	!-- Allocate array for uu
980	CASE ( 'uu' )
981	IF ( .NOT. ALLOCATED( uu ) ) THEN
982	ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) )
983	uu = 0.0_wp
984	ENDIF
985	!
986	!-- Allocate array for vv
987	CASE ( 'vv' )
988	IF ( .NOT. ALLOCATED( vv ) ) THEN
989	ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) )
990	vv = 0.0_wp
991	ENDIF
992	!
993	!-- Allocate array for ww
994	CASE ( 'ww' )
995	IF ( .NOT. ALLOCATED( ww ) ) THEN
996	ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
997	ww = 0.0_wp
998	ENDIF
999	!
1000	!-- Allocate array for 2-m potential temperature
1001	CASE ( 'theta_2m*' )
1002	IF ( .NOT. ALLOCATED( pt_2m ) ) THEN
1003	ALLOCATE( pt_2m(nys:nyn,nxl:nxr) )
1004	pt_2m = 0.0_wp
1005	ENDIF
1006	!
1007	!-- Allocate array for 10-m wind speed
1008	CASE ( 'wspeed_10m*' )
1009	IF ( .NOT. ALLOCATED( uv_10m ) ) THEN
1010	ALLOCATE( uv_10m(nys:nyn,nxl:nxr) )
1011	uv_10m = 0.0_wp
1012	ENDIF
1013
1014	END SELECT
1015
1016	ivar = ivar + 1
1017	ENDDO
1018
1019	initialized_diagnostic_output_quantities = .TRUE.
1020
1021	END SUBROUTINE doq_init
1022
1023
1024	!--------------------------------------------------------------------------------------------------!
1025	! Description:
1026	! ------------
1027	!> Calculate diagnostic quantities
1028	!--------------------------------------------------------------------------------------------------!
1029	SUBROUTINE doq_calculate
1030
1031	IMPLICIT NONE
1032
1033	INTEGER(iwp) :: i !< grid index x-dimension
1034	INTEGER(iwp) :: j !< grid index y-dimension
1035	INTEGER(iwp) :: k !< grid index z-dimension
1036	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
1037
1038	TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines
1039
1040
1041	! CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
1042
1043	!
1044	!-- Save timestep number to check in time_integration if doq_calculate
1045	!-- has been called already, since the CALL occurs at two locations, but the calculations need to be
1046	!-- done only once per timestep.
1047	timestep_number_at_prev_calc = current_timestep_number
1048
1049	ivar = 1
1050
1051	DO WHILE ( ivar <= SIZE( do_all ) )
1052
1053	SELECT CASE ( TRIM( do_all(ivar) ) )
1054	!
1055	!-- Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
1056	CASE ( 'ti' )
1057	DO i = nxl, nxr
1058	DO j = nys, nyn
1059	DO k = nzb+1, nzt
1060	ti(k,j,i) = 0.25_wp * SQRT( &
1061	( ( w(k,j+1,i) + w(k-1,j+1,i) &
1062	- w(k,j-1,i) - w(k-1,j-1,i) ) * ddy &
1063	- ( v(k+1,j,i) + v(k+1,j+1,i) &
1064	- v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2 &
1065	+ ( ( u(k+1,j,i) + u(k+1,j,i+1) &
1066	- u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k) &
1067	- ( w(k,j,i+1) + w(k-1,j,i+1) &
1068	- w(k,j,i-1) - w(k-1,j,i-1) ) * ddx )**2 &
1069	+ ( ( v(k,j,i+1) + v(k,j+1,i+1) &
1070	- v(k,j,i-1) - v(k,j+1,i-1) ) * ddx &
1071	- ( u(k,j+1,i) + u(k,j+1,i+1) &
1072	- u(k,j-1,i) - u(k,j-1,i+1) ) * ddy )**2 ) &
1073	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 0) )
1074	ENDDO
1075	ENDDO
1076	ENDDO
1077	!
1078	!-- uu
1079	CASE ( 'uu' )
1080	DO i = nxl, nxr
1081	DO j = nys, nyn
1082	DO k = nzb+1, nzt
1083	uu(k,j,i) = u(k,j,i) * u(k,j,i) &
1084	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 1) )
1085	ENDDO
1086	ENDDO
1087	ENDDO
1088	!
1089	!-- vv
1090	CASE ( 'vv' )
1091	DO i = nxl, nxr
1092	DO j = nys, nyn
1093	DO k = nzb+1, nzt
1094	vv(k,j,i) = v(k,j,i) * v(k,j,i) &
1095	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 2) )
1096	ENDDO
1097	ENDDO
1098	ENDDO
1099	!
1100	!-- ww
1101	CASE ( 'ww' )
1102	DO i = nxl, nxr
1103	DO j = nys, nyn
1104	DO k = nzb+1, nzt-1
1105	ww(k,j,i) = w(k,j,i) * w(k,j,i) &
1106	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 3) )
1107	ENDDO
1108	ENDDO
1109	ENDDO
1110	!
1111	!-- 2-m potential temperature
1112	CASE ( 'theta_2m*' )
1113	!
1114	!-- 2-m potential temperature is caluclated from surface arrays. In
1115	!-- case the 2m level is below the first grid point, MOST is applied,
1116	!-- else, linear interpolation between two vertical grid levels is
1117	!-- applied. To access all surfaces, iterate over all horizontally-
1118	!-- upward facing surface types.
1119	surf => surf_def_h(0)
1120	CALL calc_pt_2m
1121	surf => surf_lsm_h
1122	CALL calc_pt_2m
1123	surf => surf_usm_h
1124	CALL calc_pt_2m
1125	!
1126	!-- 10-m wind speed
1127	CASE ( 'wspeed_10m*' )
1128	!
1129	!-- 10-m wind speed is caluclated from surface arrays. In
1130	!-- case the 10m level is below the first grid point, MOST is applied,
1131	!-- else, linear interpolation between two vertical grid levels is
1132	!-- applied. To access all surfaces, iterate over all horizontally-
1133	!-- upward facing surface types.
1134	surf => surf_def_h(0)
1135	CALL calc_wind_10m
1136	surf => surf_lsm_h
1137	CALL calc_wind_10m
1138	surf => surf_usm_h
1139	CALL calc_wind_10m
1140
1141	END SELECT
1142
1143	ivar = ivar + 1
1144	ENDDO
1145
1146	! CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
1147
1148	!
1149	!-- The following block contains subroutines to calculate diagnostic
1150	!-- surface quantities.
1151	CONTAINS
1152	!------------------------------------------------------------------------------!
1153	! Description:
1154	! ------------
1155	!> Calculation of 2-m potential temperature.
1156	!------------------------------------------------------------------------------!
1157	SUBROUTINE calc_pt_2m
1158
1159	USE surface_layer_fluxes_mod, &
1160	ONLY: psi_h
1161
1162	IMPLICIT NONE
1163
1164	INTEGER(iwp) :: kk !< running index along the z-dimension
1165	INTEGER(iwp) :: m !< running index for surface elements
1166
1167	DO m = 1, surf%ns
1168
1169	i = surf%i(m)
1170	j = surf%j(m)
1171	k = surf%k(m)
1172	!
1173	!-- If 2-m level is below the first grid level, MOST is
1174	!-- used for calculation of 2-m temperature.
1175	IF ( surf%z_mo(m) > 2.0_wp ) THEN
1176	pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa &
1177	* ( LOG( 2.0_wp / surf%z0h(m) ) &
1178	- psi_h( 2.0_wp / surf%ol(m) ) &
1179	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
1180	!
1181	!-- If 2-m level is above the first grid level, 2-m temperature
1182	!-- is linearly interpolated between the two nearest vertical grid
1183	!-- levels. Note, since 2-m temperature is only computed for
1184	!-- horizontal upward-facing surfaces, only a vertical
1185	!-- interpolation is necessary.
1186	ELSE
1187	!
1188	!-- zw(k-1) defines the height of the surface.
1189	kk = k
1190	DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp .AND. kk <= nzt )
1191	kk = kk + 1
1192	ENDDO
1193	!
1194	!-- kk defines the index of the first grid level >= 2m.
1195	pt_2m(j,i) = pt(kk-1,j,i) + &
1196	( zw(k-1) + 2.0_wp - zu(kk-1) ) * &
1197	( pt(kk,j,i) - pt(kk-1,j,i) ) / &
1198	( zu(kk) - zu(kk-1) )
1199	ENDIF
1200
1201	ENDDO
1202
1203	END SUBROUTINE calc_pt_2m
1204
1205	!------------------------------------------------------------------------------!
1206	! Description:
1207	! ------------
1208	!> Calculation of 10-m wind speed.
1209	!------------------------------------------------------------------------------!
1210	SUBROUTINE calc_wind_10m
1211
1212	USE surface_layer_fluxes_mod, &
1213	ONLY: psi_m
1214
1215	IMPLICIT NONE
1216
1217	INTEGER(iwp) :: kk !< running index along the z-dimension
1218	INTEGER(iwp) :: m !< running index for surface elements
1219
1220	REAL(wp) :: uv_l !< wind speed at lower grid point
1221	REAL(wp) :: uv_u !< wind speed at upper grid point
1222
1223	DO m = 1, surf%ns
1224
1225	i = surf%i(m)
1226	j = surf%j(m)
1227	k = surf%k(m)
1228	!
1229	!-- If 10-m level is below the first grid level, MOST is
1230	!-- used for calculation of 10-m temperature.
1231	IF ( surf%z_mo(m) > 10.0_wp ) THEN
1232	uv_10m(j,i) = surf%us(m) / kappa &
1233	* ( LOG( 10.0_wp / surf%z0(m) ) &
1234	- psi_m( 10.0_wp / surf%ol(m) ) &
1235	+ psi_m( surf%z0(m) / surf%ol(m) ) )
1236	!
1237	!-- If 10-m level is above the first grid level, 10-m wind speed
1238	!-- is linearly interpolated between the two nearest vertical grid
1239	!-- levels. Note, since 10-m temperature is only computed for
1240	!-- horizontal upward-facing surfaces, only a vertical
1241	!-- interpolation is necessary.
1242	ELSE
1243	!
1244	!-- zw(k-1) defines the height of the surface.
1245	kk = k
1246	DO WHILE ( zu(kk) - zw(k-1) < 10.0_wp .AND. kk <= nzt )
1247	kk = kk + 1
1248	ENDDO
1249	!
1250	!-- kk defines the index of the first grid level >= 10m.
1251	uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2 &
1252	+ ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
1253
1254	uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i) + u(kk,j,i+1) ) )**2 &
1255	+ ( 0.5_wp * ( v(kk,j,i) + v(kk,j+1,i) ) )**2 )
1256
1257	uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) * &
1258	( uv_u - uv_l ) / &
1259	( zu(kk) - zu(kk-1) )
1260
1261	ENDIF
1262
1263	ENDDO
1264
1265	END SUBROUTINE calc_wind_10m
1266
1267	END SUBROUTINE doq_calculate
1268
1269
1270	!------------------------------------------------------------------------------!
1271	! Description:
1272	! ------------
1273	!> Preparation of the diagnostic output, counting of the module-specific
1274	!> output quantities and gathering of the output names.
1275	!------------------------------------------------------------------------------!
1276	SUBROUTINE doq_prepare
1277
1278
1279	USE control_parameters, &
1280	ONLY: do2d, do3d, domask, masks
1281
1282	IMPLICIT NONE
1283
1284	CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) :: do2d_var = ' ' !<
1285	!< label array for 2d output quantities
1286
1287	INTEGER(iwp) :: av !< index defining type of output, av=0 instantaneous, av=1 averaged
1288	INTEGER(iwp) :: ivar !< loop index
1289	INTEGER(iwp) :: ivar_all !< loop index
1290	INTEGER(iwp) :: l !< index for cutting string
1291	INTEGER(iwp) :: mid !< masked output running index
1292
1293	prepared_diagnostic_output_quantities = .FALSE.
1294
1295	ivar = 1
1296	ivar_all = 1
1297
1298	DO av = 0, 1
1299	!
1300	!-- Remove _xy, _xz, or _yz from string
1301	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1302	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1303	!
1304	!-- Gather 2d output quantity names.
1305	!-- Check for double occurrence of output quantity, e.g. by _xy,
1306	!-- _yz, _xz.
1307	DO WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
1308	IF ( .NOT. ANY( do_all == do2d_var(av,ivar) ) ) THEN
1309	do_all(ivar_all) = do2d_var(av,ivar)
1310	ENDIF
1311	ivar = ivar + 1
1312	ivar_all = ivar_all + 1
1313	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1314	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1315	ENDDO
1316
1317	ivar = 1
1318	!
1319	!-- Gather 3d output quantity names
1320	DO WHILE ( do3d(av,ivar)(1:1) /= ' ' )
1321	do_all(ivar_all) = do3d(av,ivar)
1322	ivar = ivar + 1
1323	ivar_all = ivar_all + 1
1324	ENDDO
1325
1326	ivar = 1
1327	!
1328	!-- Gather masked output quantity names. Also check for double output
1329	!-- e.g. by different masks.
1330	DO mid = 1, masks
1331	DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
1332	IF ( .NOT. ANY( do_all == domask(mid,av,ivar) ) ) THEN
1333	do_all(ivar_all) = domask(mid,av,ivar)
1334	ENDIF
1335
1336	ivar = ivar + 1
1337	ivar_all = ivar_all + 1
1338	ENDDO
1339	ivar = 1
1340	ENDDO
1341
1342	ENDDO
1343
1344	prepared_diagnostic_output_quantities = .TRUE.
1345
1346	END SUBROUTINE doq_prepare
1347
1348	!------------------------------------------------------------------------------!
1349	! Description:
1350	! ------------
1351	!> Subroutine reads local (subdomain) restart data
1352	!> Note: With the current structure reading of non-standard array is not
1353	!> possible
1354	!------------------------------------------------------------------------------!
1355	! SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
1356	! nxr_on_file, nynf, nync, nyn_on_file, nysf, &
1357	! nysc, nys_on_file, tmp_3d_non_standard, found )
1358	!
1359	!
1360	! USE control_parameters
1361	!
1362	! USE indices
1363	!
1364	! USE kinds
1365	!
1366	! USE pegrid
1367	!
1368	!
1369	! IMPLICIT NONE
1370	!
1371	! INTEGER(iwp) :: k !<
1372	! INTEGER(iwp) :: nxlc !<
1373	! INTEGER(iwp) :: nxlf !<
1374	! INTEGER(iwp) :: nxl_on_file !<
1375	! INTEGER(iwp) :: nxrc !<
1376	! INTEGER(iwp) :: nxrf !<
1377	! INTEGER(iwp) :: nxr_on_file !<
1378	! INTEGER(iwp) :: nync !<
1379	! INTEGER(iwp) :: nynf !<
1380	! INTEGER(iwp) :: nyn_on_file !<
1381	! INTEGER(iwp) :: nysc !<
1382	! INTEGER(iwp) :: nysf !<
1383	! INTEGER(iwp) :: nys_on_file !<
1384	!
1385	! LOGICAL, INTENT(OUT) :: found
1386	!
1387	! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions
1388	! !
1389	! !-- If temporary non-standard array for reading is already allocated,
1390	! !-- deallocate it.
1391	! IF ( ALLOCATED( tmp_3d_non_standard ) ) DEALLOCATE( tmp_3d_non_standard )
1392	!
1393	! found = .TRUE.
1394	!
1395	! SELECT CASE ( restart_string(1:length) )
1396	!
1397	! CASE ( 'ti_av' )
1398	! IF ( .NOT. ALLOCATED( ti_av ) ) THEN
1399	! ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1400	! ENDIF
1401	! IF ( k == 1 ) THEN
1402	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
1403	! nxl_on_file:nxr_on_file) )
1404	! READ ( 13 ) tmp_3d_non_standard
1405	! ENDIF
1406	! ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1407	!
1408	! CASE ( 'uu_av' )
1409	! IF ( .NOT. ALLOCATED( uu_av ) ) THEN
1410	! ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1411	! ENDIF
1412	! IF ( k == 1 ) THEN
1413	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
1414	! nxl_on_file:nxr_on_file) )
1415	! READ ( 13 ) tmp_3d_non_standard
1416	! ENDIF
1417	! uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1418	!
1419	! CASE ( 'vv_av' )
1420	! IF ( .NOT. ALLOCATED( vv_av ) ) THEN
1421	! ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1422	! ENDIF
1423	! IF ( k == 1 ) THEN
1424	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
1425	! nxl_on_file:nxr_on_file) )
1426	! READ ( 13 ) tmp_3d_non_standard
1427	! ENDIF
1428	! vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1429	!
1430	! CASE ( 'ww_av' )
1431	! IF ( .NOT. ALLOCATED( ww_av ) ) THEN
1432	! ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1433	! ENDIF
1434	! IF ( k == 1 ) THEN
1435	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
1436	! nxl_on_file:nxr_on_file) )
1437	! READ ( 13 ) tmp_3d_non_standard
1438	! ENDIF
1439	! ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1440	!
1441	!
1442	! CASE DEFAULT
1443	!
1444	! found = .FALSE.
1445	!
1446	! END SELECT
1447	!
1448	! END SUBROUTINE doq_rrd_local
1449
1450	!------------------------------------------------------------------------------!
1451	! Description:
1452	! ------------
1453	!> Subroutine writes local (subdomain) restart data
1454	!------------------------------------------------------------------------------!
1455	SUBROUTINE doq_wrd_local
1456
1457
1458	IMPLICIT NONE
1459
1460	IF ( ALLOCATED( pt_2m_av ) ) THEN
1461	CALL wrd_write_string( 'pt_2m_av' )
1462	WRITE ( 14 ) pt_2m_av
1463	ENDIF
1464
1465	IF ( ALLOCATED( ti_av ) ) THEN
1466	CALL wrd_write_string( 'ti_av' )
1467	WRITE ( 14 ) ti_av
1468	ENDIF
1469
1470	IF ( ALLOCATED( uu_av ) ) THEN
1471	CALL wrd_write_string( 'uu_av' )
1472	WRITE ( 14 ) uu_av
1473	ENDIF
1474
1475	IF ( ALLOCATED( uv_10m_av ) ) THEN
1476	CALL wrd_write_string( 'uv_10m_av' )
1477	WRITE ( 14 ) uv_10m_av
1478	ENDIF
1479
1480	IF ( ALLOCATED( vv_av ) ) THEN
1481	CALL wrd_write_string( 'vv_av' )
1482	WRITE ( 14 ) vv_av
1483	ENDIF
1484
1485	IF ( ALLOCATED( ww_av ) ) THEN
1486	CALL wrd_write_string( 'ww_av' )
1487	WRITE ( 14 ) ww_av
1488	ENDIF
1489
1490
1491	END SUBROUTINE doq_wrd_local
1492
1493
1494
1495	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 |