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

Last change on this file since 4784 was 4768, checked in by suehring, 3 years ago

Enable 3D data output also with 64-bit precision

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