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

Last change on this file since 4670 was 4583, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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