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

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

Bugfix in calculation of vertical momentum and scalar fluxes

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