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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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