source: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90

Last change on this file was 4895, checked in by suehring, 5 months ago

Remove offset in terrain-following masked output and allow only mask_k_over_surface >= 1

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