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

Last change on this file since 4691 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

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