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

Last change on this file since 4457 was 4457, checked in by raasch, 16 months ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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