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

Last change on this file since 4518 was 4518, checked in by suehring, 4 years ago

Diagnostic output: Define arrays over ghost points in order to allow for standard mpi-io treatment. By this modularization of restart-data input is possible with the module interface. Move input of restart data to doq_rrd_local. Enable mpi-io for restart data. Bugfix: add missing restart input of wtheta_av, wq_av, wu_av, and wv_av.

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