Home

Context Navigation

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

Last change on this file since 4142 was 4132, checked in by suehring, 5 years ago
Bugfix in masked data output for prognostic quantities
Property svn:keywords set to `Id`
File size: 38.1 KB

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