Home

Context Navigation

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

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