Home

Context Navigation

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

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