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

Last change on this file since 4631 was 4583, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

  • Property svn:keywords set to Id
File size: 81.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 terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: diagnostic_output_quantities_mod.f90 4583 2020-06-29 12:36:47Z suehring $
26! file re-formatted to follow the PALM coding standard
27!
28! 4560 2020-06-11 12:19:47Z suehring
29! - Bugfix in calculation of vertical momentum and scalar fluxes
30! - remove averaged output variables from PUBLIC list
31!
32! 4535 2020-05-15 12:07:23Z raasch
33! bugfix for restart data format query
34!
35! 4518 2020-05-04 15:44:28Z suehring
36! * Define arrays over ghost points in order to allow for standard mpi-io treatment. By this
37!   modularization of restart-data input is possible with the module interface.
38! * Bugfix: add missing restart input of wtheta_av, wq_av, wu_av, and wv_av.
39!
40! 4517 2020-05-03 14:29:30Z raasch
41! use statement for exchange horiz added,
42! bugfix for call of exchange horiz 2d
43!
44! 4431 2020-02-27 23:23:01Z gronemeier
45! added wspeed and wdir output; bugfix: set fill_value in case of masked output
46!
47! 4360 2020-01-07 11:25:50Z suehring
48! added output of wu, wv, wtheta and wq to enable covariance calculation according to temporal EC
49! method
50!
51! 4346 2019-12-18 11:55:56Z motisi
52! Introduction of wall_flags_total_0, which currently sets bits based on static topography
53! information used in wall_flags_static_0
54!
55! 4331 2019-12-10 18:25:02Z suehring
56! - Modularize 2-m potential temperature output
57! - New output for 10-m wind speed
58!
59! 4329 2019-12-10 15:46:36Z motisi
60! Renamed wall_flags_0 to wall_flags_static_0
61!
62! 4182 2019-08-22 15:20:23Z scharf
63! Corrected "Former revisions" section
64!
65! 4167 2019-08-16 11:01:48Z suehring
66! Changed behaviour of masked output over surface to follow terrain and ignore buildings
67! (J.Resler, T.Gronemeier)
68!
69! 4157 2019-08-14 09:19:12Z suehring
70! Initialization restructured, in order to work also when data output during spin-up is enabled.
71!
72! 4132 2019-08-02 12:34:17Z suehring
73! Bugfix in masked data output
74!
75! 4069 2019-07-01 14:05:51Z Giersch
76! Masked output running index mid has been introduced as a local variable to avoid runtime error
77! (Loop variable has been modified) in time_integration
78!
79! 4039 2019-06-18 10:32:41Z suehring
80! - Add output of uu, vv, ww to enable variance calculation according temporal EC method
81! - Allocate arrays only when they are required
82! - Formatting adjustment
83! - Rename subroutines
84! - Further modularization
85!
86! 3998 2019-05-23 13:38:11Z suehring
87! Bugfix in gathering all output strings
88!
89! 3995 2019-05-22 18:59:54Z suehring
90! Avoid compiler warnings about unused variable and fix string operation which is not allowed with
91! PGI compiler
92!
93! 3994 2019-05-22 18:08:09Z suehring
94! Initial revision
95!
96! Authors:
97! --------
98! @author Farah Kanani-Suehring
99!
100!
101! Description:
102! ------------
103!> ...
104!--------------------------------------------------------------------------------------------------!
105 MODULE diagnostic_output_quantities_mod
106
107    USE arrays_3d,                                                                                 &
108        ONLY:  ddzu,                                                                               &
109               pt,                                                                                 &
110               q,                                                                                  &
111               u,                                                                                  &
112               v,                                                                                  &
113               w,                                                                                  &
114               zu,                                                                                 &
115               zw
116
117    USE basic_constants_and_equations_mod,                                                         &
118        ONLY:  kappa, pi
119
120    USE control_parameters,                                                                        &
121        ONLY:  current_timestep_number,                                                            &
122               data_output,                                                                        &
123               message_string,                                                                     &
124               restart_data_format_output,                                                         &
125               varnamelength
126!
127!     USE cpulog,                                                                                   &
128!         ONLY:  cpu_log, log_point
129
130    USE exchange_horiz_mod,                                                                        &
131        ONLY:  exchange_horiz_2d
132
133   USE grid_variables,                                                                             &
134        ONLY:  ddx, ddy
135
136    USE indices,                                                                                   &
137        ONLY:  nbgp,                                                                               &
138               nxl,                                                                                &
139               nxlg,                                                                               &
140               nxr,                                                                                &
141               nxrg,                                                                               &
142               nyn,                                                                                &
143               nyng,                                                                               &
144               nys,                                                                                &
145               nysg,                                                                               &
146               nzb,                                                                                &
147               nzt,                                                                                &
148               wall_flags_total_0
149
150    USE kinds
151
152    USE restart_data_mpi_io_mod,                                                                   &
153        ONLY:  rd_mpi_io_check_array,                                                              &
154               rrd_mpi_io,                                                                         &
155               wrd_mpi_io
156
157    USE surface_mod,                                                                               &
158        ONLY:  surf_def_h,                                                                         &
159               surf_lsm_h,                                                                         &
160               surf_type,                                                                          &
161               surf_usm_h
162
163
164    IMPLICIT NONE
165
166    CHARACTER(LEN=varnamelength), DIMENSION(500) ::  do_all = ' '
167
168    INTEGER(iwp) ::  timestep_number_at_prev_calc = 0  !< ...at previous diagnostic output calculation
169
170    LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
171    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.    !< flag indicating whether output is p
172
173    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m     !< 2-m air potential temperature
174    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m_av  !< averaged 2-m air potential temperature
175    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m    !< horizontal wind speed at 10m
176    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m_av !< averaged horizontal wind speed at 10m
177
178    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti     !< rotation(u,v,w) aka turbulence intensity
179    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity
180    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  u_center     !< u at center of grid box
181    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  u_center_av  !< mean of u_center
182    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu           !< uu
183    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av        !< mean of uu
184    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wspeed       !< horizontal wind speed
185    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wspeed_av    !< mean of horizotal wind speed
186    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  v_center     !< v at center of grid box
187    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  v_center_av  !< mean of v_center
188    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv           !< vv
189    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av        !< mean of vv
190    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wdir         !< wind direction
191    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wdir_av      !< mean wind direction
192    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww           !< ww
193    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww_av        !< mean of ww
194    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wu           !< wu
195    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wu_av        !< mean of wu
196    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wv           !< wv
197    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wv_av        !< mean of wv
198    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wtheta       !< wtheta
199    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wtheta_av    !< mean of wtheta
200    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wq           !< wq
201    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wq_av        !< mean of wq
202
203
204    SAVE
205
206    PRIVATE
207
208!
209!-- Public variables
210    PUBLIC do_all,                                                                                 &
211           initialized_diagnostic_output_quantities,                                               &
212           prepared_diagnostic_output_quantities,                                                  &
213           timestep_number_at_prev_calc
214!
215!-- Public routines
216    PUBLIC doq_3d_data_averaging,                                                                  &
217           doq_calculate,                                                                          &
218           doq_check_data_output,                                                                  &
219           doq_define_netcdf_grid,                                                                 &
220           doq_init,                                                                               &
221           doq_output_2d,                                                                          &
222           doq_output_3d,                                                                          &
223           doq_output_mask,                                                                        &
224           doq_rrd_local,                                                                          &
225           doq_wrd_local
226
227
228    INTERFACE doq_3d_data_averaging
229       MODULE PROCEDURE doq_3d_data_averaging
230    END INTERFACE doq_3d_data_averaging
231
232    INTERFACE doq_calculate
233       MODULE PROCEDURE doq_calculate
234    END INTERFACE doq_calculate
235
236    INTERFACE doq_check_data_output
237       MODULE PROCEDURE doq_check_data_output
238    END INTERFACE doq_check_data_output
239
240    INTERFACE doq_define_netcdf_grid
241       MODULE PROCEDURE doq_define_netcdf_grid
242    END INTERFACE doq_define_netcdf_grid
243
244    INTERFACE doq_output_2d
245       MODULE PROCEDURE doq_output_2d
246    END INTERFACE doq_output_2d
247
248    INTERFACE doq_output_3d
249       MODULE PROCEDURE doq_output_3d
250    END INTERFACE doq_output_3d
251
252    INTERFACE doq_output_mask
253       MODULE PROCEDURE doq_output_mask
254    END INTERFACE doq_output_mask
255
256    INTERFACE doq_init
257       MODULE PROCEDURE doq_init
258    END INTERFACE doq_init
259
260    INTERFACE doq_prepare
261       MODULE PROCEDURE doq_prepare
262    END INTERFACE doq_prepare
263
264   INTERFACE doq_rrd_local
265       MODULE PROCEDURE doq_rrd_local_ftn
266       MODULE PROCEDURE doq_rrd_local_mpi
267    END INTERFACE doq_rrd_local
268
269    INTERFACE doq_wrd_local
270       MODULE PROCEDURE doq_wrd_local
271    END INTERFACE doq_wrd_local
272
273
274 CONTAINS
275
276!--------------------------------------------------------------------------------------------------!
277! Description:
278! ------------
279!> Sum up and time-average diagnostic output quantities as well as allocate the array necessary for
280!> storing the average.
281!--------------------------------------------------------------------------------------------------!
282 SUBROUTINE doq_3d_data_averaging( mode, variable )
283
284    USE control_parameters,                                                                        &
285        ONLY:  average_count_3d
286
287    CHARACTER (LEN=*) ::  mode     !<
288    CHARACTER (LEN=*) ::  variable !<
289
290    INTEGER(iwp) ::  i !<
291    INTEGER(iwp) ::  j !<
292    INTEGER(iwp) ::  k !<
293
294
295    IF ( mode == 'allocate' )  THEN
296
297       SELECT CASE ( TRIM( variable ) )
298
299          CASE ( 'ti' )
300             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
301                ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
302             ENDIF
303             ti_av = 0.0_wp
304
305          CASE ( 'uu' )
306             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
307                ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
308             ENDIF
309             uu_av = 0.0_wp
310
311          CASE ( 'vv' )
312             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
313                ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
314             ENDIF
315             vv_av = 0.0_wp
316
317          CASE ( 'ww' )
318             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
319                ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
320             ENDIF
321             ww_av = 0.0_wp
322
323           CASE ( 'wu' )
324             IF ( .NOT. ALLOCATED( wu_av ) )  THEN
325                ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
326             ENDIF
327             wu_av = 0.0_wp
328
329           CASE ( 'wv' )
330             IF ( .NOT. ALLOCATED( wv_av ) )  THEN
331                ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
332             ENDIF
333             wv_av = 0.0_wp
334
335           CASE ( 'wtheta' )
336             IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
337                ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
338             ENDIF
339             wtheta_av = 0.0_wp
340
341           CASE ( 'wq' )
342             IF ( .NOT. ALLOCATED( wq_av ) )  THEN
343                ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
344             ENDIF
345             wq_av = 0.0_wp
346
347          CASE ( 'theta_2m*' )
348             IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
349                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
350             ENDIF
351             pt_2m_av = 0.0_wp
352
353          CASE ( 'wspeed_10m*' )
354             IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
355                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
356             ENDIF
357             uv_10m_av = 0.0_wp
358
359          CASE ( 'wspeed' )
360             IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
361                ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
362             ENDIF
363             wspeed_av = 0.0_wp
364
365          CASE ( 'wdir' )
366             IF ( .NOT. ALLOCATED( u_center_av ) )  THEN
367                ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
368             ENDIF
369             IF ( .NOT. ALLOCATED( v_center_av ) )  THEN
370                ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
371             ENDIF
372             u_center_av = 0.0_wp
373             v_center_av = 0.0_wp
374
375          CASE DEFAULT
376             CONTINUE
377
378       END SELECT
379
380    ELSEIF ( mode == 'sum' )  THEN
381
382       SELECT CASE ( TRIM( variable ) )
383
384          CASE ( 'ti' )
385             IF ( ALLOCATED( ti_av ) )  THEN
386                DO  i = nxl, nxr
387                   DO  j = nys, nyn
388                      DO  k = nzb, nzt+1
389                         ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
390                      ENDDO
391                   ENDDO
392                ENDDO
393             ENDIF
394
395          CASE ( 'uu' )
396             IF ( ALLOCATED( uu_av ) )  THEN
397                DO  i = nxl, nxr
398                   DO  j = nys, nyn
399                      DO  k = nzb, nzt+1
400                         uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
401                      ENDDO
402                   ENDDO
403                ENDDO
404             ENDIF
405
406          CASE ( 'vv' )
407             IF ( ALLOCATED( vv_av ) )  THEN
408                DO  i = nxl, nxr
409                   DO  j = nys, nyn
410                      DO  k = nzb, nzt+1
411                         vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
412                      ENDDO
413                   ENDDO
414                ENDDO
415             ENDIF
416
417          CASE ( 'ww' )
418             IF ( ALLOCATED( ww_av ) )  THEN
419                DO  i = nxl, nxr
420                   DO  j = nys, nyn
421                      DO  k = nzb, nzt+1
422                         ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
423                      ENDDO
424                   ENDDO
425                ENDDO
426             ENDIF
427
428          CASE ( 'wu' )
429             IF ( ALLOCATED( wu_av ) )  THEN
430                DO  i = nxl, nxr
431                   DO  j = nys, nyn
432                      DO  k = nzb, nzt+1
433                         wu_av(k,j,i) = wu_av(k,j,i) + wu(k,j,i)
434                      ENDDO
435                   ENDDO
436                ENDDO
437             ENDIF
438
439          CASE ( 'wv' )
440             IF ( ALLOCATED( wv_av ) )  THEN
441                DO  i = nxl, nxr
442                   DO  j = nys, nyn
443                      DO  k = nzb, nzt+1
444                         wv_av(k,j,i) = wv_av(k,j,i) + wv(k,j,i)
445                      ENDDO
446                   ENDDO
447                ENDDO
448             ENDIF
449
450          CASE ( 'wtheta' )
451             IF ( ALLOCATED( wtheta_av ) )  THEN
452                DO  i = nxl, nxr
453                   DO  j = nys, nyn
454                      DO  k = nzb, nzt+1
455                         wtheta_av(k,j,i) = wtheta_av(k,j,i) + wtheta(k,j,i)
456                      ENDDO
457                   ENDDO
458                ENDDO
459             ENDIF
460
461          CASE ( 'wq' )
462             IF ( ALLOCATED( wq_av ) )  THEN
463                DO  i = nxl, nxr
464                   DO  j = nys, nyn
465                      DO  k = nzb, nzt+1
466                         wq_av(k,j,i) = wq_av(k,j,i) + wq(k,j,i)
467                      ENDDO
468                   ENDDO
469                ENDDO
470             ENDIF
471
472          CASE ( 'theta_2m*' )
473             IF ( ALLOCATED( pt_2m_av ) )  THEN
474                DO  i = nxl, nxr
475                   DO  j = nys, nyn
476                      pt_2m_av(j,i) = pt_2m_av(j,i) + pt_2m(j,i)
477                   ENDDO
478                ENDDO
479             ENDIF
480
481          CASE ( 'wspeed_10m*' )
482             IF ( ALLOCATED( uv_10m_av ) )  THEN
483                DO  i = nxl, nxr
484                   DO  j = nys, nyn
485                      uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i)
486                   ENDDO
487                ENDDO
488             ENDIF
489
490          CASE ( 'wspeed' )
491            IF ( ALLOCATED( wspeed_av ) )  THEN
492               DO  i = nxl, nxr
493                  DO  j = nys, nyn
494                     DO  k = nzb, nzt+1
495                         wspeed_av(k,j,i) = wspeed_av(k,j,i) + wspeed(k,j,i)
496                     ENDDO
497                  ENDDO
498               ENDDO
499            ENDIF
500
501          CASE ( 'wdir' )
502             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) )  THEN
503                DO  i = nxl, nxr
504                   DO  j = nys, nyn
505                      DO  k = nzb, nzt+1
506                        u_center_av(k,j,i) = u_center_av(k,j,i) + u_center(k,j,i)
507                        v_center_av(k,j,i) = v_center_av(k,j,i) + v_center(k,j,i)
508                      ENDDO
509                   ENDDO
510                ENDDO
511             ENDIF
512
513          CASE DEFAULT
514             CONTINUE
515
516       END SELECT
517
518    ELSEIF ( mode == 'average' )  THEN
519
520       SELECT CASE ( TRIM( variable ) )
521
522          CASE ( 'ti' )
523             IF ( ALLOCATED( ti_av ) )  THEN
524                DO  i = nxl, nxr
525                   DO  j = nys, nyn
526                      DO  k = nzb, nzt+1
527                         ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
528                      ENDDO
529                   ENDDO
530                ENDDO
531             ENDIF
532
533          CASE ( 'uu' )
534             IF ( ALLOCATED( uu_av ) )  THEN
535                DO  i = nxl, nxr
536                   DO  j = nys, nyn
537                      DO  k = nzb, nzt+1
538                         uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
539                      ENDDO
540                   ENDDO
541                ENDDO
542             ENDIF
543
544          CASE ( 'vv' )
545             IF ( ALLOCATED( vv_av ) )  THEN
546                DO  i = nxl, nxr
547                   DO  j = nys, nyn
548                      DO  k = nzb, nzt+1
549                         vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
550                      ENDDO
551                   ENDDO
552                ENDDO
553             ENDIF
554
555          CASE ( 'ww' )
556             IF ( ALLOCATED( ww_av ) )  THEN
557                DO  i = nxl, nxr
558                   DO  j = nys, nyn
559                      DO  k = nzb, nzt+1
560                         ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
561                      ENDDO
562                   ENDDO
563                ENDDO
564             ENDIF
565
566          CASE ( 'wu' )
567             IF ( ALLOCATED( wu_av ) )  THEN
568                DO  i = nxl, nxr
569                   DO  j = nys, nyn
570                      DO  k = nzb, nzt+1
571                         wu_av(k,j,i) = wu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
572                      ENDDO
573                   ENDDO
574                ENDDO
575             ENDIF
576
577          CASE ( 'wv' )
578             IF ( ALLOCATED( wv_av ) )  THEN
579                DO  i = nxl, nxr
580                   DO  j = nys, nyn
581                      DO  k = nzb, nzt+1
582                         wv_av(k,j,i) = wv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
583                      ENDDO
584                   ENDDO
585                ENDDO
586             ENDIF
587
588          CASE ( 'wtheta' )
589             IF ( ALLOCATED( wtheta_av ) )  THEN
590                DO  i = nxl, nxr
591                   DO  j = nys, nyn
592                      DO  k = nzb, nzt+1
593                         wtheta_av(k,j,i) = wtheta_av(k,j,i) / REAL( average_count_3d, KIND=wp )
594                      ENDDO
595                   ENDDO
596                ENDDO
597             ENDIF
598
599          CASE ( 'wq' )
600             IF ( ALLOCATED( wq_av ) )  THEN
601                DO  i = nxl, nxr
602                   DO  j = nys, nyn
603                      DO  k = nzb, nzt+1
604                         wq_av(k,j,i) = wq_av(k,j,i) / REAL( average_count_3d, KIND=wp )
605                      ENDDO
606                   ENDDO
607                ENDDO
608             ENDIF
609
610         CASE ( 'theta_2m*' )
611            IF ( ALLOCATED( pt_2m_av ) )  THEN
612               DO  i = nxlg, nxrg
613                  DO  j = nysg, nyng
614                     pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
615                  ENDDO
616               ENDDO
617               CALL exchange_horiz_2d( pt_2m_av )
618            ENDIF
619
620         CASE ( 'wspeed_10m*' )
621            IF ( ALLOCATED( uv_10m_av ) )  THEN
622               DO  i = nxlg, nxrg
623                  DO  j = nysg, nyng
624                     uv_10m_av(j,i) = uv_10m_av(j,i) / REAL( average_count_3d, KIND=wp )
625                  ENDDO
626               ENDDO
627               CALL exchange_horiz_2d( uv_10m_av )
628            ENDIF
629
630         CASE ( 'wspeed' )
631             IF ( ALLOCATED( wspeed_av ) )  THEN
632                DO  i = nxl, nxr
633                   DO  j = nys, nyn
634                      DO  k = nzb, nzt+1
635                         wspeed_av(k,j,i) = wspeed_av(k,j,i) / REAL( average_count_3d, KIND=wp )
636                      ENDDO
637                   ENDDO
638                ENDDO
639             ENDIF
640
641          CASE ( 'wdir' )
642             IF ( ALLOCATED( u_center_av )  .AND.  ALLOCATED( v_center_av ) )  THEN
643
644                IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
645                   ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
646                ENDIF
647                wdir_av = 0.0_wp
648
649                DO  i = nxl, nxr
650                   DO  j = nys, nyn
651                      DO  k = nzb, nzt+1
652                         u_center_av(k,j,i) = u_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
653                         v_center_av(k,j,i) = v_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
654                         wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) )          &
655                                          / pi * 180.0_wp + 180.0_wp
656                      ENDDO
657                   ENDDO
658                ENDDO
659             ENDIF
660
661       END SELECT
662
663    ENDIF
664
665
666 END SUBROUTINE doq_3d_data_averaging
667
668!--------------------------------------------------------------------------------------------------!
669! Description:
670! ------------
671!> Check data output for diagnostic output
672!--------------------------------------------------------------------------------------------------!
673 SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
674
675    IMPLICIT NONE
676
677    CHARACTER (LEN=*) ::  unit  !<
678    CHARACTER (LEN=*) ::  var   !<
679
680    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  i     !< Current element of data_output
681    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  ilen  !< Length of current entry in data_output
682    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
683
684
685    SELECT CASE ( TRIM( var ) )
686
687       CASE ( 'ti' )
688          unit = '1/s'
689
690       CASE ( 'uu' )
691          unit = 'm2/s2'
692
693       CASE ( 'vv' )
694          unit = 'm2/s2'
695
696       CASE ( 'ww' )
697          unit = 'm2/s2'
698
699       CASE ( 'wu' )
700          unit = 'm2/s2'
701
702       CASE ( 'wv' )
703          unit = 'm2/s2'
704
705       CASE ( 'wtheta' )
706          unit = 'Km/s'
707
708       CASE ( 'wq' )
709          unit = 'm/s'
710
711       CASE ( 'wspeed' )
712          unit = 'm/s'
713
714       CASE ( 'wdir' )
715          unit = 'degree'
716!
717!--    Treat horizotal cross-section output quanatities
718       CASE ( 'theta_2m*', 'wspeed_10m*' )
719!
720!--       Check if output quantity is _xy only.
721          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
722             message_string = 'illegal value for data_output: "' //                                &
723                              TRIM( var ) // '" & only 2d-horizontal ' //                          &
724                              'cross sections are allowed for this value'
725             CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
726          ENDIF
727
728          IF ( TRIM( var ) == 'theta_2m*'   )  unit = 'K'
729          IF ( TRIM( var ) == 'wspeed_10m*' )  unit = 'm/s'
730
731       CASE DEFAULT
732          unit = 'illegal'
733
734    END SELECT
735
736
737 END SUBROUTINE doq_check_data_output
738
739!--------------------------------------------------------------------------------------------------!
740!
741! Description:
742! ------------
743!> Subroutine defining appropriate grid for netcdf variables.
744!--------------------------------------------------------------------------------------------------!
745 SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
746
747    IMPLICIT NONE
748
749    CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
750    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
751    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
752    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
753
754    LOGICAL, INTENT(OUT)           ::  found       !<
755
756
757    found  = .TRUE.
758
759    SELECT CASE ( TRIM( variable ) )
760!
761!--    s grid
762       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz',                                                     &
763              'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz',                                     &
764              'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' )
765
766          grid_x = 'x'
767          grid_y = 'y'
768          grid_z = 'zu'
769!
770!--    s grid surface variables
771       CASE ( 'theta_2m*_xy', 'wspeed_10m*' )
772
773          grid_x = 'x'
774          grid_y = 'y'
775          grid_z = 'zu'
776!
777!--    u grid
778       CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
779
780          grid_x = 'xu'
781          grid_y = 'y'
782          grid_z = 'zu'
783!
784!--    v grid
785       CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz'  )
786
787          grid_x = 'x'
788          grid_y = 'yv'
789          grid_z = 'zu'
790!
791!--    w grid
792       CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz',                                                     &
793              'wu', 'wu_xy', 'wu_xz', 'wu_yz',                                                     &
794              'wv', 'wv_xy', 'wv_xz', 'wv_yz',                                                     &
795              'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz',                                     &
796              'wq', 'wq_xy', 'wq_xz', 'wq_yz' )
797
798          grid_x = 'x'
799          grid_y = 'y'
800          grid_z = 'zw'
801
802       CASE DEFAULT
803          found  = .FALSE.
804          grid_x = 'none'
805          grid_y = 'none'
806          grid_z = 'none'
807
808    END SELECT
809
810
811 END SUBROUTINE doq_define_netcdf_grid
812
813!--------------------------------------------------------------------------------------------------!
814!
815! Description:
816! ------------
817!> Subroutine defining 2D output variables
818!--------------------------------------------------------------------------------------------------!
819 SUBROUTINE doq_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do,       &
820                           fill_value )
821
822
823    IMPLICIT NONE
824
825    CHARACTER (LEN=*) ::  grid     !<
826    CHARACTER (LEN=*) ::  mode     !<
827    CHARACTER (LEN=*) ::  variable !<
828
829    INTEGER(iwp) ::  av       !< value indicating averaged or non-averaged output
830    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
831    INTEGER(iwp) ::  i        !< grid index x-direction
832    INTEGER(iwp) ::  j        !< grid index y-direction
833    INTEGER(iwp) ::  k        !< grid index z-direction
834    INTEGER(iwp) ::  nzb_do   !<
835    INTEGER(iwp) ::  nzt_do   !<
836
837    LOGICAL ::  found             !< true if variable is in list
838    LOGICAL ::  resorted          !< true if array is resorted
839    LOGICAL ::  two_d             !< flag parameter that indicates 2D variables (horizontal cross sections)
840
841    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
842
843    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
844    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
845
846
847    flag_nr  = 0
848    found    = .TRUE.
849    resorted = .FALSE.
850    two_d    = .FALSE.
851
852    SELECT CASE ( TRIM( variable ) )
853
854       CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
855           IF ( av == 0 )  THEN
856              to_be_resorted => ti
857           ELSE
858              IF ( .NOT. ALLOCATED( ti_av ) )  THEN
859                 ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
860                 ti_av = REAL( fill_value, KIND = wp )
861              ENDIF
862              to_be_resorted => ti_av
863           ENDIF
864           flag_nr = 0
865
866           IF ( mode == 'xy' )  grid = 'zu'
867
868       CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
869          IF ( av == 0 )  THEN
870             to_be_resorted => uu
871          ELSE
872             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
873                ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
874                uu_av = REAL( fill_value, KIND = wp )
875             ENDIF
876             to_be_resorted => uu_av
877          ENDIF
878          flag_nr = 1
879
880          IF ( mode == 'xy' )  grid = 'zu'
881
882       CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
883          IF ( av == 0 )  THEN
884             to_be_resorted => vv
885          ELSE
886             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
887                ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
888                vv_av = REAL( fill_value, KIND = wp )
889             ENDIF
890             to_be_resorted => vv_av
891          ENDIF
892          flag_nr = 2
893
894          IF ( mode == 'xy' )  grid = 'zu'
895
896       CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
897          IF ( av == 0 )  THEN
898             to_be_resorted => ww
899          ELSE
900             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
901                ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
902                ww_av = REAL( fill_value, KIND = wp )
903             ENDIF
904             to_be_resorted => ww_av
905          ENDIF
906          flag_nr = 3
907
908          IF ( mode == 'xy' )  grid = 'zw'
909
910       CASE ( 'wu_xy', 'wu_xz', 'wu_yz' )
911          IF ( av == 0 )  THEN
912             to_be_resorted => wu
913          ELSE
914             IF ( .NOT. ALLOCATED( wu_av ) )  THEN
915                ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
916                wu_av = REAL( fill_value, KIND = wp )
917             ENDIF
918             to_be_resorted => wu_av
919          ENDIF
920          flag_nr = 0
921
922          IF ( mode == 'xy' )  grid = 'zw'
923
924       CASE ( 'wv_xy', 'wv_xz', 'wv_yz' )
925          IF ( av == 0 )  THEN
926             to_be_resorted => wv
927          ELSE
928             IF ( .NOT. ALLOCATED( wv_av ) )  THEN
929                ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
930                wv_av = REAL( fill_value, KIND = wp )
931             ENDIF
932             to_be_resorted => wv_av
933          ENDIF
934          flag_nr = 0
935
936          IF ( mode == 'xy' )  grid = 'zw'
937
938       CASE ( 'wtheta_xy', 'wtheta_xz', 'wtheta_yz' )
939          IF ( av == 0 )  THEN
940             to_be_resorted => wtheta
941          ELSE
942             IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
943                ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
944                wtheta_av = REAL( fill_value, KIND = wp )
945             ENDIF
946             to_be_resorted => wtheta_av
947          ENDIF
948          flag_nr = 0
949
950          IF ( mode == 'xy' )  grid = 'zw'
951
952       CASE ( 'wq_xy', 'wq_xz', 'wq_yz' )
953          IF ( av == 0 )  THEN
954             to_be_resorted => wq
955          ELSE
956             IF ( .NOT. ALLOCATED( wq_av ) )  THEN
957                ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
958                wq_av = REAL( fill_value, KIND = wp )
959             ENDIF
960             to_be_resorted => wq_av
961          ENDIF
962          flag_nr = 0
963
964          IF ( mode == 'xy' )  grid = 'zw'
965
966       CASE ( 'theta_2m*_xy' )        ! 2d-array
967          IF ( av == 0 )  THEN
968             DO  i = nxl, nxr
969                DO  j = nys, nyn
970                   local_pf(i,j,nzb+1) = pt_2m(j,i)
971                ENDDO
972             ENDDO
973          ELSE
974             IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
975                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
976                pt_2m_av = REAL( fill_value, KIND = wp )
977             ENDIF
978             DO  i = nxl, nxr
979                DO  j = nys, nyn
980                   local_pf(i,j,nzb+1) = pt_2m_av(j,i)
981                ENDDO
982             ENDDO
983          ENDIF
984          resorted = .TRUE.
985          two_d    = .TRUE.
986          grid     = 'zu1'
987
988       CASE ( 'wspeed_10m*_xy' )        ! 2d-array
989          IF ( av == 0 )  THEN
990             DO  i = nxl, nxr
991                DO  j = nys, nyn
992                   local_pf(i,j,nzb+1) = uv_10m(j,i)
993                ENDDO
994             ENDDO
995          ELSE
996             IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
997                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
998                uv_10m_av = REAL( fill_value, KIND = wp )
999             ENDIF
1000             DO  i = nxl, nxr
1001                DO  j = nys, nyn
1002                   local_pf(i,j,nzb+1) = uv_10m_av(j,i)
1003                ENDDO
1004             ENDDO
1005          ENDIF
1006          resorted = .TRUE.
1007          two_d    = .TRUE.
1008          grid     = 'zu1'
1009
1010       CASE ( 'wspeed_xy', 'wspeed_xz', 'wspeed_yz' )
1011          IF ( av == 0 )  THEN
1012             to_be_resorted => wspeed
1013          ELSE
1014             IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
1015                ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1016                wspeed_av = REAL( fill_value, KIND = wp )
1017             ENDIF
1018             to_be_resorted => wspeed_av
1019          ENDIF
1020          flag_nr = 0
1021
1022          IF ( mode == 'xy' )  grid = 'zu'
1023
1024       CASE ( 'wdir_xy', 'wdir_xz', 'wdir_yz' )
1025          IF ( av == 0 )  THEN
1026             to_be_resorted => wdir
1027          ELSE
1028             IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
1029                ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1030                wdir_av = REAL( fill_value, KIND = wp )
1031             ENDIF
1032             to_be_resorted => wdir_av
1033          ENDIF
1034          flag_nr = 0
1035
1036          IF ( mode == 'xy' )  grid = 'zu'
1037
1038       CASE DEFAULT
1039          found = .FALSE.
1040          grid  = 'none'
1041
1042    END SELECT
1043
1044    IF ( found  .AND.  .NOT. resorted )  THEN
1045       DO  i = nxl, nxr
1046          DO  j = nys, nyn
1047             DO  k = nzb_do, nzt_do
1048                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                                    &
1049                                         REAL( fill_value, KIND = wp ),                            &
1050                                         BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
1051             ENDDO
1052          ENDDO
1053       ENDDO
1054    ENDIF
1055
1056 END SUBROUTINE doq_output_2d
1057
1058
1059!--------------------------------------------------------------------------------------------------!
1060!
1061! Description:
1062! ------------
1063!> Subroutine defining 3D output variables
1064!--------------------------------------------------------------------------------------------------!
1065 SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1066
1067    IMPLICIT NONE
1068
1069    CHARACTER (LEN=*) ::  variable !<
1070
1071    INTEGER(iwp) ::  av       !< index indicating averaged or instantaneous output
1072    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
1073    INTEGER(iwp) ::  i        !< index variable along x-direction
1074    INTEGER(iwp) ::  j        !< index variable along y-direction
1075    INTEGER(iwp) ::  k        !< index variable along z-direction
1076    INTEGER(iwp) ::  nzb_do   !< lower limit of the data output (usually 0)
1077    INTEGER(iwp) ::  nzt_do   !< vertical upper limit of the data output (usually nz_do3d)
1078
1079    LOGICAL ::  found             !< true if variable is in list
1080    LOGICAL ::  resorted          !< true if array is resorted
1081
1082    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
1083
1084    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf        !<
1085    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
1086
1087
1088    flag_nr  = 0
1089    found    = .TRUE.
1090    resorted = .FALSE.
1091
1092    SELECT CASE ( TRIM( variable ) )
1093
1094       CASE ( 'ti' )
1095          IF ( av == 0 )  THEN
1096             to_be_resorted => ti
1097          ELSE
1098             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
1099                ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1100                ti_av = REAL( fill_value, KIND = wp )
1101             ENDIF
1102             to_be_resorted => ti_av
1103          ENDIF
1104          flag_nr = 0
1105
1106       CASE ( 'uu' )
1107          IF ( av == 0 )  THEN
1108             to_be_resorted => uu
1109          ELSE
1110             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
1111                ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1112                uu_av = REAL( fill_value, KIND = wp )
1113             ENDIF
1114             to_be_resorted => uu_av
1115          ENDIF
1116          flag_nr = 1
1117
1118       CASE ( 'vv' )
1119          IF ( av == 0 )  THEN
1120             to_be_resorted => vv
1121          ELSE
1122             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
1123                ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1124                vv_av = REAL( fill_value, KIND = wp )
1125             ENDIF
1126             to_be_resorted => vv_av
1127          ENDIF
1128          flag_nr = 2
1129
1130       CASE ( 'ww' )
1131          IF ( av == 0 )  THEN
1132             to_be_resorted => ww
1133          ELSE
1134             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
1135                ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1136                ww_av = REAL( fill_value, KIND = wp )
1137             ENDIF
1138             to_be_resorted => ww_av
1139          ENDIF
1140          flag_nr = 3
1141
1142       CASE ( 'wu' )
1143          IF ( av == 0 )  THEN
1144             to_be_resorted => wu
1145          ELSE
1146             IF ( .NOT. ALLOCATED( wu_av ) )  THEN
1147                ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1148                wu_av = REAL( fill_value, KIND = wp )
1149             ENDIF
1150             to_be_resorted => wu_av
1151          ENDIF
1152          flag_nr = 0
1153
1154       CASE ( 'wv' )
1155          IF ( av == 0 )  THEN
1156             to_be_resorted => wv
1157          ELSE
1158             IF ( .NOT. ALLOCATED( wv_av ) )  THEN
1159                ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1160                wv_av = REAL( fill_value, KIND = wp )
1161             ENDIF
1162             to_be_resorted => wv_av
1163          ENDIF
1164          flag_nr = 0
1165
1166       CASE ( 'wtheta' )
1167          IF ( av == 0 )  THEN
1168             to_be_resorted => wtheta
1169          ELSE
1170             IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
1171                ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1172                wtheta_av = REAL( fill_value, KIND = wp )
1173             ENDIF
1174             to_be_resorted => wtheta_av
1175          ENDIF
1176          flag_nr = 0
1177
1178       CASE ( 'wq' )
1179          IF ( av == 0 )  THEN
1180             to_be_resorted => wq
1181          ELSE
1182             IF ( .NOT. ALLOCATED( wq_av ) )  THEN
1183                ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1184                wq_av = REAL( fill_value, KIND = wp )
1185             ENDIF
1186             to_be_resorted => wq_av
1187          ENDIF
1188          flag_nr = 0
1189
1190       CASE ( 'wspeed' )
1191          IF ( av == 0 )  THEN
1192             to_be_resorted => wspeed
1193          ELSE
1194             IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
1195                ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1196                wspeed_av = REAL( fill_value, KIND = wp )
1197             ENDIF
1198             to_be_resorted => wspeed_av
1199          ENDIF
1200          flag_nr = 0
1201
1202       CASE ( 'wdir' )
1203          IF ( av == 0 )  THEN
1204             to_be_resorted => wdir
1205          ELSE
1206             IF ( .NOT. ALLOCATED( wdir_av ) )  THEN
1207                ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1208                wdir_av = REAL( fill_value, KIND = wp )
1209             ENDIF
1210             to_be_resorted => wdir_av
1211          ENDIF
1212          flag_nr = 0
1213
1214       CASE DEFAULT
1215          found = .FALSE.
1216
1217    END SELECT
1218
1219    IF ( found  .AND.  .NOT. resorted )  THEN
1220       DO  i = nxl, nxr
1221          DO  j = nys, nyn
1222             DO  k = nzb_do, nzt_do
1223                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                                    &
1224                                         REAL( fill_value, KIND = wp ),                            &
1225                                         BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
1226             ENDDO
1227          ENDDO
1228       ENDDO
1229    ENDIF
1230
1231 END SUBROUTINE doq_output_3d
1232
1233
1234!--------------------------------------------------------------------------------------------------!
1235!
1236! Description:
1237! ------------
1238!> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with indices
1239!> (i,j,k) for masked data output.
1240!--------------------------------------------------------------------------------------------------!
1241 SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
1242
1243    USE control_parameters
1244
1245    USE indices
1246
1247    IMPLICIT NONE
1248
1249    REAL(wp), PARAMETER   ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
1250
1251    CHARACTER (LEN=*) ::  variable   !<
1252    CHARACTER (LEN=5) ::  grid       !< flag to distinquish between staggered grids
1253
1254    INTEGER(iwp) ::  av              !< index indicating averaged or instantaneous output
1255    INTEGER(iwp) ::  flag_nr         !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
1256    INTEGER(iwp) ::  i               !< index variable along x-direction
1257    INTEGER(iwp) ::  j               !< index variable along y-direction
1258    INTEGER(iwp) ::  k               !< index variable along z-direction
1259    INTEGER(iwp) ::  im              !< loop index for masked variables
1260    INTEGER(iwp) ::  jm              !< loop index for masked variables
1261    INTEGER(iwp) ::  kk              !< masked output index variable along z-direction
1262    INTEGER(iwp) ::  mid             !< masked output running index
1263    INTEGER(iwp) ::  ktt             !< k index of highest horizontal surface
1264
1265    LOGICAL      ::  found           !< true if variable is in list
1266    LOGICAL      ::  resorted        !< true if array is resorted
1267
1268    REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  local_pf   !<
1269
1270    REAL(wp), DIMENSION(:,:,:), POINTER  ::  to_be_resorted  !< points to array which needs to be resorted for output
1271
1272
1273    flag_nr  = 0
1274    found    = .TRUE.
1275    resorted = .FALSE.
1276    grid     = 's'
1277
1278    SELECT CASE ( TRIM( variable ) )
1279
1280       CASE ( 'ti' )
1281          IF ( av == 0 )  THEN
1282             to_be_resorted => ti
1283          ELSE
1284             to_be_resorted => ti_av
1285          ENDIF
1286          grid = 's'
1287          flag_nr = 0
1288
1289       CASE ( 'uu' )
1290          IF ( av == 0 )  THEN
1291             to_be_resorted => uu
1292          ELSE
1293             to_be_resorted => uu_av
1294          ENDIF
1295          grid = 'u'
1296          flag_nr = 1
1297
1298       CASE ( 'vv' )
1299          IF ( av == 0 )  THEN
1300             to_be_resorted => vv
1301          ELSE
1302             to_be_resorted => vv_av
1303          ENDIF
1304          grid = 'v'
1305          flag_nr = 2
1306
1307       CASE ( 'ww' )
1308          IF ( av == 0 )  THEN
1309             to_be_resorted => ww
1310          ELSE
1311             to_be_resorted => ww_av
1312          ENDIF
1313          grid = 'w'
1314          flag_nr = 3
1315
1316       CASE ( 'wu' )
1317          IF ( av == 0 )  THEN
1318             to_be_resorted => wu
1319          ELSE
1320             to_be_resorted => wu_av
1321          ENDIF
1322          grid = 's'
1323          flag_nr = 0
1324
1325       CASE ( 'wv' )
1326          IF ( av == 0 )  THEN
1327             to_be_resorted => wv
1328          ELSE
1329             to_be_resorted => wv_av
1330          ENDIF
1331          grid = 's'
1332          flag_nr = 0
1333
1334       CASE ( 'wtheta' )
1335          IF ( av == 0 )  THEN
1336             to_be_resorted => wtheta
1337          ELSE
1338             to_be_resorted => wtheta_av
1339          ENDIF
1340          grid = 's'
1341          flag_nr = 0
1342
1343       CASE ( 'wq' )
1344          IF ( av == 0 )  THEN
1345             to_be_resorted => wq
1346          ELSE
1347             to_be_resorted => wq_av
1348          ENDIF
1349          grid = 's'
1350          flag_nr = 0
1351
1352       CASE ( 'wspeed' )
1353          IF ( av == 0 )  THEN
1354             to_be_resorted => wspeed
1355          ELSE
1356             to_be_resorted => wspeed_av
1357          ENDIF
1358          grid = 's'
1359          flag_nr = 0
1360
1361       CASE ( 'wdir' )
1362          IF ( av == 0 )  THEN
1363             to_be_resorted => wdir
1364          ELSE
1365             to_be_resorted => wdir_av
1366          ENDIF
1367          grid = 's'
1368          flag_nr = 0
1369
1370       CASE DEFAULT
1371          found = .FALSE.
1372
1373    END SELECT
1374
1375    IF ( found  .AND.  .NOT. resorted )  THEN
1376       IF ( .NOT. mask_surface(mid) )  THEN
1377!
1378!--       Default masked output
1379          DO  i = 1, mask_size_l(mid,1)
1380             DO  j = 1, mask_size_l(mid,2)
1381                DO  k = 1, mask_size_l(mid,3)
1382                   local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k),                          &
1383                                                           mask_j(mid,j),                          &
1384                                                           mask_i(mid,i)),                         &
1385                                            REAL( fill_value, KIND = wp ),                         &
1386                                            BTEST( wall_flags_total_0(mask_k(mid,k),               &
1387                                                                      mask_j(mid,j),               &
1388                                                                      mask_i(mid,i)),              &
1389                                                   flag_nr ) )
1390                ENDDO
1391             ENDDO
1392          ENDDO
1393
1394       ELSE
1395!
1396!--       Terrain-following masked output
1397          DO  i = 1, mask_size_l(mid,1)
1398             DO  j = 1, mask_size_l(mid,2)
1399!
1400!--             Get k index of the highest terraing surface
1401                im = mask_i(mid,i)
1402                jm = mask_j(mid,j)
1403                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ), DIM=1 ) - 1
1404                DO  k = 1, mask_size_l(mid,3)
1405                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
1406!
1407!--                Set value if not in building
1408                   IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
1409                      local_pf(i,j,k) = fill_value
1410                   ELSE
1411                      local_pf(i,j,k) = to_be_resorted(kk,jm,im)
1412                   ENDIF
1413                ENDDO
1414             ENDDO
1415          ENDDO
1416
1417       ENDIF
1418    ENDIF
1419
1420 END SUBROUTINE doq_output_mask
1421
1422!--------------------------------------------------------------------------------------------------!
1423! Description:
1424! ------------
1425!> Allocate required arrays
1426!--------------------------------------------------------------------------------------------------!
1427 SUBROUTINE doq_init
1428
1429    IMPLICIT NONE
1430
1431    INTEGER(iwp) ::  ivar   !< loop index over all 2d/3d/mask output quantities
1432
1433
1434!
1435!-- Next line is to avoid compiler warnings about unused variables
1436    IF ( timestep_number_at_prev_calc == 0 )  CONTINUE
1437!
1438!-- Preparatory steps and initialization of output arrays
1439    IF ( .NOT.  prepared_diagnostic_output_quantities )  CALL doq_prepare
1440
1441    initialized_diagnostic_output_quantities = .FALSE.
1442
1443    ivar = 1
1444
1445    DO  WHILE ( ivar <= SIZE( do_all ) )
1446
1447       SELECT CASE ( TRIM( do_all(ivar) ) )
1448!
1449!--       Allocate array for 'turbulence intensity'
1450          CASE ( 'ti' )
1451             IF ( .NOT. ALLOCATED( ti ) )  THEN
1452                ALLOCATE( ti(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1453                ti = 0.0_wp
1454             ENDIF
1455!
1456!--       Allocate array for uu
1457          CASE ( 'uu' )
1458             IF ( .NOT. ALLOCATED( uu ) )  THEN
1459                ALLOCATE( uu(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1460                uu = 0.0_wp
1461             ENDIF
1462!
1463!--       Allocate array for vv
1464          CASE ( 'vv' )
1465             IF ( .NOT. ALLOCATED( vv ) )  THEN
1466                ALLOCATE( vv(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1467                vv = 0.0_wp
1468             ENDIF
1469!
1470!--       Allocate array for ww
1471          CASE ( 'ww' )
1472             IF ( .NOT. ALLOCATED( ww ) )  THEN
1473                ALLOCATE( ww(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1474                ww = 0.0_wp
1475             ENDIF
1476!
1477!--       Allocate array for wu
1478          CASE ( 'wu' )
1479             IF ( .NOT. ALLOCATED( wu ) )  THEN
1480                ALLOCATE( wu(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1481                wu = 0.0_wp
1482             ENDIF
1483!
1484!--       Allocate array for wv
1485          CASE ( 'wv' )
1486             IF ( .NOT. ALLOCATED( wv ) )  THEN
1487                ALLOCATE( wv(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1488                wv = 0.0_wp
1489             ENDIF
1490!
1491!--       Allocate array for wtheta
1492          CASE ( 'wtheta' )
1493             IF ( .NOT. ALLOCATED( wtheta ) )  THEN
1494                ALLOCATE( wtheta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1495                wtheta = 0.0_wp
1496             ENDIF
1497!
1498!--       Allocate array for wq
1499          CASE ( 'wq' )
1500             IF ( .NOT. ALLOCATED( wq ) )  THEN
1501                ALLOCATE( wq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1502                wq = 0.0_wp
1503             ENDIF
1504!
1505!--       Allocate array for 2-m potential temperature
1506          CASE ( 'theta_2m*' )
1507             IF ( .NOT. ALLOCATED( pt_2m ) )  THEN
1508                ALLOCATE( pt_2m(nysg:nyng,nxlg:nxrg) )
1509                pt_2m = 0.0_wp
1510             ENDIF
1511!
1512!--       Allocate array for 10-m wind speed
1513          CASE ( 'wspeed_10m*' )
1514             IF ( .NOT. ALLOCATED( uv_10m ) )  THEN
1515                ALLOCATE( uv_10m(nysg:nyng,nxlg:nxrg) )
1516                uv_10m = 0.0_wp
1517             ENDIF
1518!
1519!--       Allocate array for wspeed
1520          CASE ( 'wspeed' )
1521             IF ( .NOT. ALLOCATED( wspeed ) )  THEN
1522                ALLOCATE( wspeed(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1523                wspeed = 0.0_wp
1524             ENDIF
1525
1526!
1527!--       Allocate array for wdir
1528          CASE ( 'wdir' )
1529             IF ( .NOT. ALLOCATED( u_center ) )  THEN
1530                ALLOCATE( u_center(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1531                u_center = 0.0_wp
1532             ENDIF
1533             IF ( .NOT. ALLOCATED( v_center ) )  THEN
1534                ALLOCATE( v_center(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1535                v_center = 0.0_wp
1536             ENDIF
1537             IF ( .NOT. ALLOCATED( wdir ) )  THEN
1538                ALLOCATE( wdir(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1539                wdir = 0.0_wp
1540             ENDIF
1541
1542       END SELECT
1543
1544       ivar = ivar + 1
1545    ENDDO
1546
1547    initialized_diagnostic_output_quantities = .TRUE.
1548
1549 END SUBROUTINE doq_init
1550
1551
1552!--------------------------------------------------------------------------------------------------!
1553! Description:
1554! ------------
1555!> Calculate diagnostic quantities
1556!--------------------------------------------------------------------------------------------------!
1557 SUBROUTINE doq_calculate
1558
1559    IMPLICIT NONE
1560
1561    INTEGER(iwp) ::  i          !< grid index x-dimension
1562    INTEGER(iwp) ::  j          !< grid index y-dimension
1563    INTEGER(iwp) ::  k          !< grid index z-dimension
1564    INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
1565
1566    TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
1567
1568
1569!     CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
1570
1571!
1572!-- Save timestep number to check in time_integration if doq_calculate has been called already,
1573!-- since the CALL occurs at two locations, but the calculations need to be done only once per
1574!-- timestep.
1575    timestep_number_at_prev_calc = current_timestep_number
1576
1577    ivar = 1
1578
1579    DO  WHILE ( ivar <= SIZE( do_all ) )
1580
1581       SELECT CASE ( TRIM( do_all(ivar) ) )
1582!
1583!--       Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
1584          CASE ( 'ti' )
1585             DO  i = nxl, nxr
1586                DO  j = nys, nyn
1587                   DO  k = nzb+1, nzt
1588                      ti(k,j,i) = 0.25_wp * SQRT(                                                  &
1589                                       (   (   w(k,j+1,i) + w(k-1,j+1,i)                           &
1590                                             - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy                   &
1591                                         - (   v(k+1,j,i) + v(k+1,j+1,i)                           &
1592                                             - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2          &
1593                                     + (   (   u(k+1,j,i) + u(k+1,j,i+1)                           &
1594                                             - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)               &
1595                                         - (   w(k,j,i+1) + w(k-1,j,i+1)                           &
1596                                             - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2          &
1597                                     + (   (   v(k,j,i+1) + v(k,j+1,i+1)                           &
1598                                             - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx                   &
1599                                         - (   u(k,j+1,i) + u(k,j+1,i+1)                           &
1600                                             - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )       &
1601                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
1602                   ENDDO
1603                ENDDO
1604             ENDDO
1605!
1606!--       uu
1607          CASE ( 'uu' )
1608             DO  i = nxl, nxr
1609                DO  j = nys, nyn
1610                   DO  k = nzb+1, nzt
1611                      uu(k,j,i) = u(k,j,i) * u(k,j,i)                                              &
1612                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
1613                   ENDDO
1614                ENDDO
1615             ENDDO
1616!
1617!--       vv
1618          CASE ( 'vv' )
1619             DO  i = nxl, nxr
1620                DO  j = nys, nyn
1621                   DO  k = nzb+1, nzt
1622                      vv(k,j,i) = v(k,j,i) * v(k,j,i)                                              &
1623                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
1624                   ENDDO
1625                ENDDO
1626             ENDDO
1627!
1628!--       ww
1629          CASE ( 'ww' )
1630             DO  i = nxl, nxr
1631                DO  j = nys, nyn
1632                   DO  k = nzb+1, nzt-1
1633                      ww(k,j,i) = w(k,j,i) * w(k,j,i)                                              &
1634                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
1635                   ENDDO
1636                ENDDO
1637             ENDDO
1638!
1639!--       wu
1640          CASE ( 'wu' )
1641             DO  i = nxl, nxr
1642                DO  j = nys, nyn
1643                   DO  k = nzb+1, nzt-1
1644                      wu(k,j,i) = w(k,j,i)                                                         &
1645                                  * 0.25_wp * ( u(k,j,i) + u(k,j,i+1) + u(k+1,j,i) + u(k+1,j,i+1) )&
1646                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
1647                   ENDDO
1648                ENDDO
1649             ENDDO
1650!
1651!--       wv
1652          CASE ( 'wv' )
1653             DO  i = nxl, nxr
1654                DO  j = nys, nyn
1655                   DO  k = nzb+1, nzt-1
1656                      wv(k,j,i) = w(k,j,i)                                                         &
1657                                  * 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k+1,j,i) + v(k+1,j+1,i) )&
1658                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
1659                   ENDDO
1660                ENDDO
1661             ENDDO
1662!
1663!--       wtheta
1664          CASE ( 'wtheta' )
1665             DO  i = nxl, nxr
1666                DO  j = nys, nyn
1667                   DO  k = nzb+1, nzt-1
1668                      wtheta(k,j,i) = w(k,j,i)                                                     &
1669                                     *  0.5_wp  * ( pt(k,j,i) + pt(k+1,j,i) )                      &
1670                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ))
1671                   ENDDO
1672                ENDDO
1673             ENDDO
1674!
1675!--       wq
1676          CASE ( 'wq' )
1677             DO  i = nxl, nxr
1678                DO  j = nys, nyn
1679                   DO  k = nzb+1, nzt-1
1680                      wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) )                    &
1681                                  * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
1682                   ENDDO
1683                ENDDO
1684             ENDDO
1685!
1686!--       2-m potential temperature
1687          CASE ( 'theta_2m*' )
1688!
1689!--          2-m potential temperature is caluclated from surface arrays. In case the 2m level is
1690!--          below the first grid point, MOST is applied, else, linear interpolation between two
1691!--          vertical grid levels is applied. To access all surfaces, iterate over all horizontally-
1692!--          upward facing surface types.
1693             surf => surf_def_h(0)
1694             CALL calc_pt_2m
1695             surf => surf_lsm_h
1696             CALL calc_pt_2m
1697             surf => surf_usm_h
1698             CALL calc_pt_2m
1699!
1700!--       10-m wind speed
1701          CASE ( 'wspeed_10m*' )
1702!
1703!--          10-m wind speed is caluclated from surface arrays. In case the 10m level is below the
1704!--          first grid point, MOST is applied, else, linear interpolation between two vertical grid
1705!--          levels is applied. To access all surfaces, iterate over all horizontally-upward facing
1706!--          surface types.
1707             surf => surf_def_h(0)
1708             CALL calc_wind_10m
1709             surf => surf_lsm_h
1710             CALL calc_wind_10m
1711             surf => surf_usm_h
1712             CALL calc_wind_10m
1713!
1714!--       Horizontal wind speed
1715          CASE ( 'wspeed' )
1716             DO  i = nxl, nxr
1717                DO  j = nys, nyn
1718                   DO  k = nzb, nzt+1
1719                      wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2              &
1720                                          + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 )            &
1721                                     * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ))
1722                   ENDDO
1723                ENDDO
1724             ENDDO
1725
1726!
1727!--       Horizontal wind direction
1728          CASE ( 'wdir' )
1729             DO  i = nxl, nxr
1730                DO  j = nys, nyn
1731                   DO  k = nzb, nzt+1
1732                      u_center(k,j,i) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
1733                      v_center(k,j,i) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
1734
1735                      wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) )                      &
1736                                    / pi * 180.0_wp + 180.0_wp
1737                   ENDDO
1738                ENDDO
1739             ENDDO
1740
1741       END SELECT
1742
1743       ivar = ivar + 1
1744    ENDDO
1745
1746!     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
1747
1748!
1749!-- The following block contains subroutines to calculate diagnostic surface quantities.
1750    CONTAINS
1751!--------------------------------------------------------------------------------------------------!
1752! Description:
1753! ------------
1754!> Calculation of 2-m potential temperature.
1755!--------------------------------------------------------------------------------------------------!
1756       SUBROUTINE calc_pt_2m
1757
1758          USE surface_layer_fluxes_mod,                                                            &
1759              ONLY:  psi_h
1760
1761          IMPLICIT NONE
1762
1763          INTEGER(iwp) ::  kk     !< running index along the z-dimension
1764          INTEGER(iwp) ::  m      !< running index for surface elements
1765
1766          DO  m = 1, surf%ns
1767
1768             i = surf%i(m)
1769             j = surf%j(m)
1770             k = surf%k(m)
1771!
1772!--          If 2-m level is below the first grid level, MOST is used for calculation of
1773!--          2-m temperature.
1774             IF ( surf%z_mo(m) > 2.0_wp )  THEN
1775                pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa                               &
1776                                * ( LOG( 2.0_wp / surf%z0h(m) )                                    &
1777                                    - psi_h( 2.0_wp      / surf%ol(m) )                            &
1778                                    + psi_h( surf%z0h(m) / surf%ol(m) ) )
1779!
1780!--          If 2-m level is above the first grid level, 2-m temperature is linearly interpolated
1781!--          between the two nearest vertical grid levels. Note, since 2-m temperature is only
1782!--          computed for horizontal upward-facing surfaces, only a vertical interpolation is
1783!--          necessary.
1784             ELSE
1785!
1786!--             zw(k-1) defines the height of the surface.
1787                kk = k
1788                DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp  .AND.  kk <= nzt )
1789                   kk = kk + 1
1790                ENDDO
1791!
1792!--             kk defines the index of the first grid level >= 2m.
1793                pt_2m(j,i) = pt(kk-1,j,i) +                                                        &
1794                              ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *                                &
1795                              ( pt(kk,j,i)       - pt(kk-1,j,i) ) /                                &
1796                              ( zu(kk)           - zu(kk-1)     )
1797             ENDIF
1798
1799          ENDDO
1800
1801       END SUBROUTINE calc_pt_2m
1802
1803!--------------------------------------------------------------------------------------------------!
1804! Description:
1805! ------------
1806!> Calculation of 10-m wind speed.
1807!--------------------------------------------------------------------------------------------------!
1808       SUBROUTINE calc_wind_10m
1809
1810          USE surface_layer_fluxes_mod,                                                            &
1811              ONLY:  psi_m
1812
1813          IMPLICIT NONE
1814
1815          INTEGER(iwp) ::  kk     !< running index along the z-dimension
1816          INTEGER(iwp) ::  m      !< running index for surface elements
1817
1818          REAL(wp) ::  uv_l !< wind speed at lower grid point
1819          REAL(wp) ::  uv_u !< wind speed at upper grid point
1820
1821
1822          DO  m = 1, surf%ns
1823
1824             i = surf%i(m)
1825             j = surf%j(m)
1826             k = surf%k(m)
1827!
1828!--          If 10-m level is below the first grid level, MOST is used for calculation of 10-m
1829!--          temperature.
1830             IF ( surf%z_mo(m) > 10.0_wp )  THEN
1831                uv_10m(j,i) = surf%us(m) / kappa                                                   &
1832                              * ( LOG( 10.0_wp /  surf%z0(m) )                                     &
1833                                  - psi_m( 10.0_wp    / surf%ol(m) )                               &
1834                                  + psi_m( surf%z0(m) / surf%ol(m) ) )
1835!
1836!--          If 10-m level is above the first grid level, 10-m wind speed is linearly interpolated
1837!--          between the two nearest vertical grid levels. Note, since 10-m temperature is only
1838!--          computed for horizontal upward-facing surfaces, only a vertical interpolation is
1839!--          necessary.
1840             ELSE
1841!
1842!--             zw(k-1) defines the height of the surface.
1843                kk = k
1844                DO WHILE ( zu(kk) - zw(k-1) < 10.0_wp  .AND.  kk <= nzt )
1845                   kk = kk + 1
1846                ENDDO
1847!
1848!--             kk defines the index of the first grid level >= 10m.
1849                uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2                       &
1850                           + ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
1851
1852                uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i)   + u(kk,j,i+1)   ) )**2                       &
1853                           + ( 0.5_wp * ( v(kk,j,i)   + v(kk,j+1,i)   ) )**2 )
1854
1855                uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) *                            &
1856                                     ( uv_u              - uv_l     ) /                            &
1857                                     ( zu(kk)            - zu(kk-1) )
1858
1859             ENDIF
1860
1861          ENDDO
1862
1863       END SUBROUTINE calc_wind_10m
1864
1865 END SUBROUTINE doq_calculate
1866
1867
1868!--------------------------------------------------------------------------------------------------!
1869! Description:
1870! ------------
1871!> Preparation of the diagnostic output, counting of the module-specific output quantities and
1872!> gathering of the output names.
1873!--------------------------------------------------------------------------------------------------!
1874 SUBROUTINE doq_prepare
1875
1876    USE control_parameters,                                                                        &
1877        ONLY:  do2d, do3d, domask, masks
1878
1879    IMPLICIT NONE
1880
1881    CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d_var = ' '  !< label array for 2d output quantities
1882
1883    INTEGER(iwp) ::  av         !< index defining type of output, av=0 instantaneous, av=1 averaged
1884    INTEGER(iwp) ::  ivar       !< loop index
1885    INTEGER(iwp) ::  ivar_all   !< loop index
1886    INTEGER(iwp) ::  l          !< index for cutting string
1887    INTEGER(iwp) ::  mid          !< masked output running index
1888
1889
1890    prepared_diagnostic_output_quantities = .FALSE.
1891
1892    ivar     = 1
1893    ivar_all = 1
1894
1895    DO  av = 0, 1
1896!
1897!--    Remove _xy, _xz, or _yz from string
1898       l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1899       do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1900!
1901!--    Gather 2d output quantity names.
1902!--    Check for double occurrence of output quantity, e.g. by _xy, _yz, _xz.
1903       DO  WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
1904          IF ( .NOT.  ANY( do_all == do2d_var(av,ivar) ) )  THEN
1905             do_all(ivar_all) = do2d_var(av,ivar)
1906          ENDIF
1907          ivar = ivar + 1
1908          ivar_all = ivar_all + 1
1909          l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1910          do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1911       ENDDO
1912
1913       ivar = 1
1914!
1915!--    Gather 3d output quantity names
1916       DO  WHILE ( do3d(av,ivar)(1:1) /= ' ' )
1917          do_all(ivar_all) = do3d(av,ivar)
1918          ivar = ivar + 1
1919          ivar_all = ivar_all + 1
1920       ENDDO
1921
1922       ivar = 1
1923!
1924!--    Gather masked output quantity names. Also check for double output e.g. by different masks.
1925       DO  mid = 1, masks
1926          DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
1927             IF ( .NOT.  ANY( do_all == domask(mid,av,ivar) ) )  THEN
1928                do_all(ivar_all) = domask(mid,av,ivar)
1929             ENDIF
1930
1931             ivar = ivar + 1
1932             ivar_all = ivar_all + 1
1933          ENDDO
1934          ivar = 1
1935       ENDDO
1936
1937    ENDDO
1938
1939    prepared_diagnostic_output_quantities = .TRUE.
1940
1941 END SUBROUTINE doq_prepare
1942
1943!--------------------------------------------------------------------------------------------------!
1944! Description:
1945! ------------
1946!> Subroutine reads local (subdomain) restart data
1947!--------------------------------------------------------------------------------------------------!
1948 SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,    &
1949                               nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, found )
1950
1951    USE control_parameters
1952
1953    USE indices
1954
1955    USE kinds
1956
1957    USE pegrid
1958
1959
1960    IMPLICIT NONE
1961
1962    INTEGER(iwp) ::  k               !<
1963    INTEGER(iwp) ::  nxlc            !<
1964    INTEGER(iwp) ::  nxlf            !<
1965    INTEGER(iwp) ::  nxl_on_file     !<
1966    INTEGER(iwp) ::  nxrc            !<
1967    INTEGER(iwp) ::  nxrf            !<
1968    INTEGER(iwp) ::  nxr_on_file     !<
1969    INTEGER(iwp) ::  nync            !<
1970    INTEGER(iwp) ::  nynf            !<
1971    INTEGER(iwp) ::  nyn_on_file     !<
1972    INTEGER(iwp) ::  nysc            !<
1973    INTEGER(iwp) ::  nysf            !<
1974    INTEGER(iwp) ::  nys_on_file     !<
1975
1976    LOGICAL, INTENT(OUT)  :: found
1977
1978    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
1979
1980    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
1981
1982
1983    found = .TRUE.
1984
1985    SELECT CASE ( restart_string(1:length) )
1986
1987       CASE ( 'pt_2m_av' )
1988          IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
1989             ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
1990          ENDIF
1991          IF ( k == 1 )  READ ( 13 )  tmp_2d
1992          pt_2m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =                                     &
1993                      tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
1994
1995       CASE ( 'ti_av' )
1996          IF ( .NOT. ALLOCATED( ti_av ) )  THEN
1997             ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1998          ENDIF
1999          IF ( k == 1 )  READ ( 13 )  tmp_3d
2000          ti_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
2001                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2002
2003       CASE ( 'u_center_av' )
2004          IF ( .NOT. ALLOCATED( u_center_av ) )  THEN
2005             ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2006          ENDIF
2007          IF ( k == 1 )  READ ( 13 )  tmp_3d
2008          u_center_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                 &
2009                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2010
2011       CASE ( 'uu_av' )
2012          IF ( .NOT. ALLOCATED( uu_av ) )  THEN
2013             ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2014          ENDIF
2015          IF ( k == 1 )  READ ( 13 )  tmp_3d
2016          uu_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
2017                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2018
2019       CASE ( 'uv_10m_av' )
2020          IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
2021             ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
2022          ENDIF
2023          IF ( k == 1 )  READ ( 13 )  tmp_2d
2024          uv_10m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =                                    &
2025                       tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2026
2027       CASE ( 'v_center_av' )
2028          IF ( .NOT. ALLOCATED( v_center_av ) )  THEN
2029             ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2030          ENDIF
2031          IF ( k == 1 )  READ ( 13 )  tmp_3d
2032          v_center_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                 &
2033                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2034
2035       CASE ( 'vv_av' )
2036          IF ( .NOT. ALLOCATED( vv_av ) )  THEN
2037             ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2038          ENDIF
2039          IF ( k == 1 )  READ ( 13 )  tmp_3d
2040          vv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
2041                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2042
2043       CASE ( 'wtheta_av' )
2044          IF ( .NOT. ALLOCATED( wtheta_av ) )  THEN
2045             ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2046          ENDIF
2047          IF ( k == 1 )  READ ( 13 )  tmp_3d
2048          wtheta_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                   &
2049                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2050
2051       CASE ( 'wq_av' )
2052          IF ( .NOT. ALLOCATED( wq_av ) )  THEN
2053             ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2054          ENDIF
2055          IF ( k == 1 )  READ ( 13 )  tmp_3d
2056          wq_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
2057                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2058
2059       CASE ( 'wspeed_av' )
2060          IF ( .NOT. ALLOCATED( wspeed_av ) )  THEN
2061             ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2062          ENDIF
2063          IF ( k == 1 )  READ ( 13 )  tmp_3d
2064          wspeed_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                   &
2065                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2066
2067       CASE ( 'wu_av' )
2068          IF ( .NOT. ALLOCATED( wu_av ) )  THEN
2069             ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2070          ENDIF
2071          IF ( k == 1 )  READ ( 13 )  tmp_3d
2072          wu_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
2073                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2074
2075       CASE ( 'wv_av' )
2076          IF ( .NOT. ALLOCATED( wv_av ) )  THEN
2077             ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2078          ENDIF
2079          IF ( k == 1 )  READ ( 13 )  tmp_3d
2080          wv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
2081                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2082
2083       CASE ( 'ww_av' )
2084          IF ( .NOT. ALLOCATED( ww_av ) )  THEN
2085             ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2086          ENDIF
2087          IF ( k == 1 )  READ ( 13 )  tmp_3d
2088          ww_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                       &
2089                      tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2090
2091
2092       CASE DEFAULT
2093
2094          found = .FALSE.
2095
2096    END SELECT
2097
2098 END SUBROUTINE doq_rrd_local_ftn
2099
2100!--------------------------------------------------------------------------------------------------!
2101! Description:
2102! ------------
2103!> Read module-specific local restart data arrays (MPI-IO).
2104!--------------------------------------------------------------------------------------------------!
2105 SUBROUTINE doq_rrd_local_mpi
2106
2107    LOGICAL ::  array_found  !< control flad indicating whether the array is found in the file or not
2108
2109
2110    CALL rd_mpi_io_check_array( 'pt_2m_av' , found = array_found )
2111    IF ( array_found )  THEN
2112       IF ( .NOT. ALLOCATED( pt_2m_av ) )  ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
2113       CALL rrd_mpi_io( 'pt_2m_av', pt_2m_av )
2114    ENDIF
2115
2116    CALL rd_mpi_io_check_array( 'ti_av' , found = array_found )
2117    IF ( array_found )  THEN
2118       IF ( .NOT. ALLOCATED( ti_av ) )  ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2119       CALL rrd_mpi_io( 'ti_av', ti_av )
2120    ENDIF
2121
2122    CALL rd_mpi_io_check_array( 'u_center_av' , found = array_found )
2123    IF ( array_found )  THEN
2124       IF ( .NOT. ALLOCATED( u_center_av ) )  ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2125       CALL rrd_mpi_io( 'u_center_av', u_center_av )
2126    ENDIF
2127
2128    CALL rd_mpi_io_check_array( 'uu_av' , found = array_found )
2129    IF ( array_found )  THEN
2130       IF ( .NOT. ALLOCATED( uu_av ) )  ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2131       CALL rrd_mpi_io( 'uu_av', uu_av )
2132    ENDIF
2133
2134    CALL rd_mpi_io_check_array( 'uv_10m_av' , found = array_found )
2135    IF ( array_found )  THEN
2136       IF ( .NOT. ALLOCATED( uv_10m_av ) )  ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
2137       CALL rrd_mpi_io( 'uv_10m_av', uv_10m_av )
2138    ENDIF
2139
2140    CALL rd_mpi_io_check_array( 'v_center_av' , found = array_found )
2141    IF ( array_found )  THEN
2142       IF ( .NOT. ALLOCATED( v_center_av ) )  ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2143       CALL rrd_mpi_io( 'v_center_av', v_center_av )
2144    ENDIF
2145
2146    CALL rd_mpi_io_check_array( 'vv_av' , found = array_found )
2147    IF ( array_found )  THEN
2148       IF ( .NOT. ALLOCATED( vv_av ) )  ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2149       CALL rrd_mpi_io( 'vv_av', vv_av )
2150    ENDIF
2151
2152    CALL rd_mpi_io_check_array( 'ww_av' , found = array_found )
2153    IF ( array_found )  THEN
2154       IF ( .NOT. ALLOCATED( ww_av ) )  ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2155       CALL rrd_mpi_io( 'ww_av', ww_av )
2156    ENDIF
2157
2158    CALL rd_mpi_io_check_array( 'wu_av' , found = array_found )
2159    IF ( array_found )  THEN
2160       IF ( .NOT. ALLOCATED( wu_av ) )  ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2161       CALL rrd_mpi_io( 'wu_av', wu_av )
2162    ENDIF
2163
2164    CALL rd_mpi_io_check_array( 'wv_av' , found = array_found )
2165    IF ( array_found )  THEN
2166       IF ( .NOT. ALLOCATED( wv_av ) )  ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2167       CALL rrd_mpi_io( 'wv_av', wv_av )
2168    ENDIF
2169
2170    CALL rd_mpi_io_check_array( 'wtheta_av' , found = array_found )
2171    IF ( array_found )  THEN
2172       IF ( .NOT. ALLOCATED( wtheta_av ) )  ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2173       CALL rrd_mpi_io( 'wtheta_av', wtheta_av )
2174    ENDIF
2175
2176    CALL rd_mpi_io_check_array( 'wq_av' , found = array_found )
2177    IF ( array_found )  THEN
2178       IF ( .NOT. ALLOCATED( wq_av ) )  ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2179       CALL rrd_mpi_io( 'wq_av', wq_av )
2180    ENDIF
2181
2182    CALL rd_mpi_io_check_array( 'wspeed_av' , found = array_found )
2183    IF ( array_found )  THEN
2184       IF ( .NOT. ALLOCATED( wspeed_av ) )  ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2185       CALL rrd_mpi_io( 'wspeed_av', wspeed_av )
2186    ENDIF
2187
2188 END SUBROUTINE doq_rrd_local_mpi
2189
2190!--------------------------------------------------------------------------------------------------!
2191! Description:
2192! ------------
2193!> Subroutine writes local (subdomain) restart data
2194!--------------------------------------------------------------------------------------------------!
2195 SUBROUTINE doq_wrd_local
2196
2197    IMPLICIT NONE
2198
2199
2200    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
2201
2202       IF ( ALLOCATED( pt_2m_av ) )  THEN
2203          CALL wrd_write_string( 'pt_2m_av' )
2204          WRITE ( 14 )  pt_2m_av
2205       ENDIF
2206
2207       IF ( ALLOCATED( ti_av ) )  THEN
2208          CALL wrd_write_string( 'ti_av' )
2209          WRITE ( 14 )  ti_av
2210       ENDIF
2211
2212       IF ( ALLOCATED( u_center_av ) )  THEN
2213          CALL wrd_write_string( 'u_center_av' )
2214          WRITE ( 14 )  u_center_av
2215       ENDIF
2216
2217       IF ( ALLOCATED( uu_av ) )  THEN
2218          CALL wrd_write_string( 'uu_av' )
2219          WRITE ( 14 )  uu_av
2220       ENDIF
2221
2222       IF ( ALLOCATED( uv_10m_av ) )  THEN
2223          CALL wrd_write_string( 'uv_10m_av' )
2224          WRITE ( 14 )  uv_10m_av
2225       ENDIF
2226
2227       IF ( ALLOCATED( v_center_av ) )  THEN
2228          CALL wrd_write_string( 'v_center_av' )
2229          WRITE ( 14 )  v_center_av
2230       ENDIF
2231
2232       IF ( ALLOCATED( vv_av ) )  THEN
2233          CALL wrd_write_string( 'vv_av' )
2234          WRITE ( 14 )  vv_av
2235       ENDIF
2236
2237       IF ( ALLOCATED( ww_av ) )  THEN
2238          CALL wrd_write_string( 'ww_av' )
2239          WRITE ( 14 )  ww_av
2240       ENDIF
2241
2242       IF ( ALLOCATED( wu_av ) )  THEN
2243          CALL wrd_write_string( 'wu_av' )
2244          WRITE ( 14 )  wu_av
2245       ENDIF
2246
2247       IF ( ALLOCATED( wv_av ) )  THEN
2248          CALL wrd_write_string( 'wv_av' )
2249          WRITE ( 14 )  wv_av
2250       ENDIF
2251
2252       IF ( ALLOCATED( wtheta_av ) )  THEN
2253          CALL wrd_write_string( 'wtheta_av' )
2254          WRITE ( 14 )  wtheta_av
2255       ENDIF
2256
2257       IF ( ALLOCATED( wq_av ) )  THEN
2258          CALL wrd_write_string( 'wq_av' )
2259          WRITE ( 14 )  wq_av
2260       ENDIF
2261
2262       IF ( ALLOCATED( wspeed_av ) )  THEN
2263          CALL wrd_write_string( 'wspeed_av' )
2264          WRITE ( 14 )  wspeed_av
2265       ENDIF
2266
2267    ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
2268
2269       IF ( ALLOCATED( pt_2m_av ) )     CALL wrd_mpi_io( 'pt_2m_av', pt_2m_av )
2270       IF ( ALLOCATED( ti_av ) )        CALL wrd_mpi_io( 'ti_av', ti_av )
2271       IF ( ALLOCATED( u_center_av ) )  CALL wrd_mpi_io( 'u_center_av', u_center_av )
2272       IF ( ALLOCATED( uu_av ) )        CALL wrd_mpi_io( 'uu_av', uu_av )
2273       IF ( ALLOCATED( uv_10m_av ) )    CALL wrd_mpi_io( 'uv_10m_av', uv_10m_av )
2274       IF ( ALLOCATED( vv_av ) )        CALL wrd_mpi_io( 'vv_av', vv_av )
2275       IF ( ALLOCATED( v_center_av ) )  CALL wrd_mpi_io( 'v_center_av', v_center_av )
2276       IF ( ALLOCATED( wtheta_av ) )    CALL wrd_mpi_io( 'wtheta_av', wtheta_av )
2277       IF ( ALLOCATED( wq_av ) )        CALL wrd_mpi_io( 'wq_av', wq_av )
2278       IF ( ALLOCATED( wspeed_av ) )    CALL wrd_mpi_io( 'wspeed_av', wspeed_av )
2279       IF ( ALLOCATED( wu_av ) )        CALL wrd_mpi_io( 'wu_av', wu_av )
2280       IF ( ALLOCATED( wv_av ) )        CALL wrd_mpi_io( 'wv_av', wv_av )
2281       IF ( ALLOCATED( ww_av ) )        CALL wrd_mpi_io( 'ww_av', ww_av )
2282
2283    ENDIF
2284
2285 END SUBROUTINE doq_wrd_local
2286
2287
2288 END MODULE diagnostic_output_quantities_mod
Note: See TracBrowser for help on using the repository browser.