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

Last change on this file since 4455 was 4431, checked in by gronemeier, 4 years ago

diagnostic_output_quantities: added wspeed and wdir output; bugfix: set fill_value in case of masked output

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