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

Last change on this file since 4346 was 4346, checked in by motisi, 16 months ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

  • Property svn:keywords set to Id
File size: 50.9 KB
Line 
1!> @file diagnostic_output_quantities_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diagnostic_output_quantities_mod.f90 4346 2019-12-18 11:55:56Z motisi $
27! Introduction of wall_flags_total_0, which currently sets bits based on static
28! topography information used in wall_flags_static_0
29!
30! 4331 2019-12-10 18:25:02Z suehring
31! - Modularize 2-m potential temperature output
32! - New output for 10-m wind speed
33!
34! 4329 2019-12-10 15:46:36Z motisi
35! Renamed wall_flags_0 to wall_flags_static_0
36!
37! 4182 2019-08-22 15:20:23Z scharf
38! Corrected "Former revisions" section
39!
40! 4167 2019-08-16 11:01:48Z suehring
41! Changed behaviour of masked output over surface to follow terrain and ignore
42! buildings (J.Resler, T.Gronemeier)
43!
44! 4157 2019-08-14 09:19:12Z suehring
45! Initialization restructured, in order to work also when data output during
46! spin-up is enabled.
47!
48! 4132 2019-08-02 12:34:17Z suehring
49! Bugfix in masked data output
50!
51! 4069 2019-07-01 14:05:51Z Giersch
52! Masked output running index mid has been introduced as a local variable to
53! avoid runtime error (Loop variable has been modified) in time_integration
54!
55! 4039 2019-06-18 10:32:41Z suehring
56! - Add output of uu, vv, ww to enable variance calculation according temporal
57!   EC method
58! - Allocate arrays only when they are required
59! - Formatting adjustment
60! - Rename subroutines
61! - Further modularization
62!
63! 3998 2019-05-23 13:38:11Z suehring
64! Bugfix in gathering all output strings
65!
66! 3995 2019-05-22 18:59:54Z suehring
67! Avoid compiler warnings about unused variable and fix string operation which
68! is not allowed with PGI compiler
69!
70! 3994 2019-05-22 18:08:09Z suehring
71! Initial revision
72!
73! Authors:
74! --------
75! @author Farah Kanani-Suehring
76!
77!
78! Description:
79! ------------
80!> ...
81!------------------------------------------------------------------------------!
82 MODULE diagnostic_output_quantities_mod
83 
84    USE arrays_3d,                                                             &
85        ONLY:  ddzu,                                                           &
86               pt,                                                             &
87               u,                                                              &
88               v,                                                              &
89               w,                                                              &
90               zu,                                                             &
91               zw
92
93    USE basic_constants_and_equations_mod,                                     &
94        ONLY:  kappa
95
96    USE control_parameters,                                                    &
97        ONLY:  current_timestep_number,                                        &
98               data_output,                                                    &
99               message_string,                                                 &
100               varnamelength
101!
102!     USE cpulog,                                                                &
103!         ONLY:  cpu_log, log_point
104
105   USE grid_variables,                                                         &
106        ONLY:  ddx, ddy
107
108    USE indices,                                                               &
109        ONLY:  nbgp,                                                           &
110               nxl,                                                            &
111               nxlg,                                                           & 
112               nxr,                                                            & 
113               nxrg,                                                           &
114               nyn,                                                            &
115               nyng,                                                           &
116               nys,                                                            &
117               nysg,                                                           &
118               nzb,                                                            &
119               nzt,                                                            &
120               wall_flags_total_0
121
122    USE kinds
123
124    USE surface_mod,                                                           &
125        ONLY:  surf_def_h,                                                     &
126               surf_lsm_h,                                                     &
127               surf_type,                                                      &
128               surf_usm_h
129
130
131    IMPLICIT NONE
132
133    CHARACTER(LEN=varnamelength), DIMENSION(500) ::  do_all = ' '
134
135    INTEGER(iwp) ::  timestep_number_at_prev_calc = 0  !< ...at previous diagnostic output calculation
136 
137    LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
138    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.    !< flag indicating whether output is p
139
140    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m     !< 2-m air potential temperature
141    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_2m_av  !< averaged 2-m air potential temperature
142    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m    !< horizontal wind speed at 10m
143    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  uv_10m_av !< averaged horizontal wind speed at 10m
144
145    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti     !< rotation(u,v,w) aka turbulence intensity
146    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity   
147    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu     !< uu
148    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av  !< mean of uu   
149    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv     !< vv
150    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av  !< mean of vv   
151    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww     !< ww
152    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww_av  !< mean of ww
153
154
155    SAVE
156
157    PRIVATE
158
159!
160!-- Public variables
161    PUBLIC do_all,                                                             &
162           initialized_diagnostic_output_quantities,                           &
163           prepared_diagnostic_output_quantities,                              &
164           timestep_number_at_prev_calc,                                       &
165           pt_2m_av,                                                           &
166           ti_av,                                                              &
167           uu_av,                                                              &
168           uv_10m_av,                                                          &
169           vv_av,                                                              &
170           ww_av
171!                                                                             
172!-- Public routines                                                           
173    PUBLIC doq_3d_data_averaging,                                              &
174           doq_calculate,                                                      &
175           doq_check_data_output,                                              &
176           doq_define_netcdf_grid,                                             &
177           doq_output_2d,                                                      &
178           doq_output_3d,                                                      &
179           doq_output_mask,                                                    &
180           doq_init,                                                           &
181           doq_wrd_local
182!          doq_rrd_local,                                                      &
183
184
185    INTERFACE doq_3d_data_averaging
186       MODULE PROCEDURE doq_3d_data_averaging
187    END INTERFACE doq_3d_data_averaging       
188
189    INTERFACE doq_calculate
190       MODULE PROCEDURE doq_calculate
191    END INTERFACE doq_calculate
192
193    INTERFACE doq_check_data_output
194       MODULE PROCEDURE doq_check_data_output
195    END INTERFACE doq_check_data_output
196   
197    INTERFACE doq_define_netcdf_grid
198       MODULE PROCEDURE doq_define_netcdf_grid
199    END INTERFACE doq_define_netcdf_grid
200   
201    INTERFACE doq_output_2d
202       MODULE PROCEDURE doq_output_2d
203    END INTERFACE doq_output_2d
204   
205    INTERFACE doq_output_3d
206       MODULE PROCEDURE doq_output_3d
207    END INTERFACE doq_output_3d
208   
209    INTERFACE doq_output_mask
210       MODULE PROCEDURE doq_output_mask
211    END INTERFACE doq_output_mask
212     
213    INTERFACE doq_init
214       MODULE PROCEDURE doq_init
215    END INTERFACE doq_init
216
217    INTERFACE doq_prepare
218       MODULE PROCEDURE doq_prepare
219    END INTERFACE doq_prepare
220   
221!     INTERFACE doq_rrd_local
222!        MODULE PROCEDURE doq_rrd_local
223!     END INTERFACE doq_rrd_local
224   
225    INTERFACE doq_wrd_local
226       MODULE PROCEDURE doq_wrd_local
227    END INTERFACE doq_wrd_local
228
229
230 CONTAINS
231 
232!------------------------------------------------------------------------------!
233! Description:
234! ------------
235!> Sum up and time-average diagnostic output quantities as well as allocate
236!> the array necessary for storing the average.
237!------------------------------------------------------------------------------!
238 SUBROUTINE doq_3d_data_averaging( mode, variable )
239
240    USE control_parameters,                                                    &
241        ONLY:  average_count_3d
242
243    CHARACTER (LEN=*) ::  mode     !<
244    CHARACTER (LEN=*) ::  variable !<
245
246    INTEGER(iwp) ::  i !<
247    INTEGER(iwp) ::  j !<
248    INTEGER(iwp) ::  k !<
249
250    IF ( mode == 'allocate' )  THEN
251
252       SELECT CASE ( TRIM( variable ) )
253
254          CASE ( 'ti' )
255             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
256                ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
257             ENDIF
258             ti_av = 0.0_wp
259       
260          CASE ( 'uu' )
261             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
262                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
263             ENDIF
264             uu_av = 0.0_wp
265               
266          CASE ( 'vv' )
267             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
268                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
269             ENDIF
270             vv_av = 0.0_wp
271               
272          CASE ( 'ww' )
273             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
274                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
275             ENDIF
276             ww_av = 0.0_wp
277             
278          CASE ( 'theta_2m*' )
279             IF ( .NOT. ALLOCATED( pt_2m_av ) )  THEN
280                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
281             ENDIF
282             pt_2m_av = 0.0_wp
283             
284          CASE ( 'wspeed_10m*' )
285             IF ( .NOT. ALLOCATED( uv_10m_av ) )  THEN
286                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
287             ENDIF
288             uv_10m_av = 0.0_wp
289
290          CASE DEFAULT
291             CONTINUE
292
293       END SELECT
294
295    ELSEIF ( mode == 'sum' )  THEN
296
297       SELECT CASE ( TRIM( variable ) )
298 
299          CASE ( 'ti' )
300             IF ( ALLOCATED( ti_av ) ) THEN
301                DO  i = nxl, nxr
302                   DO  j = nys, nyn
303                      DO  k = nzb, nzt+1
304                         ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
305                      ENDDO
306                   ENDDO
307                ENDDO
308             ENDIF
309             
310          CASE ( 'uu' )
311             IF ( ALLOCATED( uu_av ) ) THEN
312                DO  i = nxl, nxr
313                   DO  j = nys, nyn
314                      DO  k = nzb, nzt+1
315                         uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
316                      ENDDO
317                   ENDDO
318                ENDDO
319             ENDIF
320             
321          CASE ( 'vv' )
322             IF ( ALLOCATED( vv_av ) ) THEN
323                DO  i = nxl, nxr
324                   DO  j = nys, nyn
325                      DO  k = nzb, nzt+1
326                         vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
327                      ENDDO
328                   ENDDO
329                ENDDO
330             ENDIF
331             
332          CASE ( 'ww' )
333             IF ( ALLOCATED( ww_av ) ) THEN
334                DO  i = nxl, nxr
335                   DO  j = nys, nyn
336                      DO  k = nzb, nzt+1
337                         ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
338                      ENDDO
339                   ENDDO
340                ENDDO
341             ENDIF
342             
343          CASE ( 'theta_2m*' )
344             IF ( ALLOCATED( pt_2m_av ) ) THEN
345                DO  i = nxl, nxr
346                   DO  j = nys, nyn
347                      pt_2m_av(j,i) = pt_2m_av(j,i) + pt_2m(j,i)
348                   ENDDO
349                ENDDO
350             ENDIF
351
352          CASE ( 'wspeed_10m*' )
353             IF ( ALLOCATED( uv_10m_av ) ) THEN
354                DO  i = nxl, nxr
355                   DO  j = nys, nyn
356                      uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i)
357                   ENDDO
358                ENDDO
359             ENDIF
360
361          CASE DEFAULT
362             CONTINUE
363
364       END SELECT
365
366    ELSEIF ( mode == 'average' )  THEN
367
368       SELECT CASE ( TRIM( variable ) )
369
370          CASE ( 'ti' )
371             IF ( ALLOCATED( ti_av ) ) THEN
372                DO  i = nxl, nxr
373                   DO  j = nys, nyn
374                      DO  k = nzb, nzt+1
375                         ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
376                      ENDDO
377                   ENDDO
378                ENDDO
379             ENDIF
380       
381          CASE ( 'uu' )
382             IF ( ALLOCATED( uu_av ) ) THEN
383                DO  i = nxl, nxr
384                   DO  j = nys, nyn
385                      DO  k = nzb, nzt+1
386                         uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
387                      ENDDO
388                   ENDDO
389                ENDDO
390             ENDIF
391             
392          CASE ( 'vv' )
393             IF ( ALLOCATED( vv_av ) ) THEN
394                DO  i = nxl, nxr
395                   DO  j = nys, nyn
396                      DO  k = nzb, nzt+1
397                         vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
398                      ENDDO
399                   ENDDO
400                ENDDO
401             ENDIF
402             
403          CASE ( 'ww' )
404             IF ( ALLOCATED( ww_av ) ) THEN
405                DO  i = nxl, nxr
406                   DO  j = nys, nyn
407                      DO  k = nzb, nzt+1
408                         ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
409                      ENDDO
410                   ENDDO
411                ENDDO
412             ENDIF
413
414         CASE ( 'theta_2m*' )
415            IF ( ALLOCATED( pt_2m_av ) ) THEN
416               DO  i = nxlg, nxrg
417                  DO  j = nysg, nyng
418                     pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
419                  ENDDO
420               ENDDO
421               CALL exchange_horiz_2d( pt_2m_av, nbgp )
422            ENDIF
423
424         CASE ( 'wspeed_10m*' )
425            IF ( ALLOCATED( uv_10m_av ) ) THEN
426               DO  i = nxlg, nxrg
427                  DO  j = nysg, nyng
428                     uv_10m_av(j,i) = uv_10m_av(j,i) / REAL( average_count_3d, KIND=wp )
429                  ENDDO
430               ENDDO
431               CALL exchange_horiz_2d( uv_10m_av, nbgp )
432            ENDIF
433
434       END SELECT
435
436    ENDIF
437
438
439 END SUBROUTINE doq_3d_data_averaging
440 
441!------------------------------------------------------------------------------!
442! Description:
443! ------------
444!> Check data output for diagnostic output
445!------------------------------------------------------------------------------!
446 SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
447
448    IMPLICIT NONE
449
450    CHARACTER (LEN=*) ::  unit  !<
451    CHARACTER (LEN=*) ::  var   !<
452
453    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  i     !< Current element of data_output
454    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  ilen  !< Length of current entry in data_output
455    INTEGER(iwp), OPTIONAL, INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
456
457    SELECT CASE ( TRIM( var ) )
458
459       CASE ( 'ti' )
460          unit = '1/s'
461             
462       CASE ( 'uu' )
463          unit = 'm2/s2'
464             
465       CASE ( 'vv' )
466          unit = 'm2/s2'
467             
468       CASE ( 'ww' )
469          unit = 'm2/s2'
470!
471!--    Treat horizotal cross-section output quanatities
472       CASE ( 'theta_2m*', 'wspeed_10m*' )
473!
474!--       Check if output quantity is _xy only.
475          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
476             message_string = 'illegal value for data_output: "' //            &
477                              TRIM( var ) // '" & only 2d-horizontal ' //      &
478                              'cross sections are allowed for this value'
479             CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
480          ENDIF
481
482          IF ( TRIM( var ) == 'theta_2m*'   )  unit = 'K'
483          IF ( TRIM( var ) == 'wspeed_10m*' )  unit = 'm/s'
484
485       CASE DEFAULT
486          unit = 'illegal'
487
488    END SELECT
489
490
491 END SUBROUTINE doq_check_data_output
492 
493!------------------------------------------------------------------------------!
494!
495! Description:
496! ------------
497!> Subroutine defining appropriate grid for netcdf variables.
498!------------------------------------------------------------------------------!
499 SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
500   
501    IMPLICIT NONE
502
503    CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
504    LOGICAL, INTENT(OUT)           ::  found       !<
505    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
506    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
507    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
508
509    found  = .TRUE.
510   
511    SELECT CASE ( TRIM( variable ) )
512!
513!--    s grid
514       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz'  )
515
516          grid_x = 'x'
517          grid_y = 'y'
518          grid_z = 'zu'
519!
520!--    s grid surface variables
521       CASE ( 'theta_2m*_xy', 'wspeed_10m*' )
522
523          grid_x = 'x'
524          grid_y = 'y'
525          grid_z = 'zu'
526!
527!--    u grid
528       CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
529
530          grid_x = 'xu'
531          grid_y = 'y'
532          grid_z = 'zu'
533!
534!--    v grid
535       CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz'  )
536
537          grid_x = 'x'
538          grid_y = 'yv'
539          grid_z = 'zu'
540!
541!--    w grid
542       CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz'  )
543
544          grid_x = 'x'
545          grid_y = 'y'
546          grid_z = 'zw'
547
548       CASE DEFAULT
549          found  = .FALSE.
550          grid_x = 'none'
551          grid_y = 'none'
552          grid_z = 'none'
553
554    END SELECT
555
556
557 END SUBROUTINE doq_define_netcdf_grid
558 
559!------------------------------------------------------------------------------!
560!
561! Description:
562! ------------
563!> Subroutine defining 2D output variables
564!------------------------------------------------------------------------------!
565 SUBROUTINE doq_output_2d( av, variable, found, grid,                          &
566                           mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
567
568
569    IMPLICIT NONE
570
571    CHARACTER (LEN=*) ::  grid     !<
572    CHARACTER (LEN=*) ::  mode     !<
573    CHARACTER (LEN=*) ::  variable !<
574
575    INTEGER(iwp) ::  av       !< value indicating averaged or non-averaged output
576    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
577    INTEGER(iwp) ::  i        !< grid index x-direction
578    INTEGER(iwp) ::  j        !< grid index y-direction
579    INTEGER(iwp) ::  k        !< grid index z-direction
580    INTEGER(iwp) ::  nzb_do   !<
581    INTEGER(iwp) ::  nzt_do   !<
582
583    LOGICAL ::  found             !< true if variable is in list
584    LOGICAL ::  resorted          !< true if array is resorted
585    LOGICAL ::  two_d             !< flag parameter that indicates 2D variables (horizontal cross sections)
586
587    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
588   
589    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
590    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
591   
592    flag_nr  = 0
593    found    = .TRUE.
594    resorted = .FALSE.
595    two_d    = .FALSE.
596
597    SELECT CASE ( TRIM( variable ) )
598
599       CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
600           IF ( av == 0 )  THEN
601              to_be_resorted => ti
602           ELSE
603              IF ( .NOT. ALLOCATED( ti_av ) ) THEN
604                 ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
605                 ti_av = REAL( fill_value, KIND = wp )
606              ENDIF
607              to_be_resorted => ti_av
608           ENDIF
609           flag_nr = 0
610           
611           IF ( mode == 'xy' )  grid = 'zu'
612   
613       CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
614          IF ( av == 0 )  THEN
615             to_be_resorted => uu
616          ELSE
617             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
618                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
619                uu_av = REAL( fill_value, KIND = wp )
620             ENDIF
621             to_be_resorted => uu_av
622          ENDIF
623          flag_nr = 1
624         
625          IF ( mode == 'xy' )  grid = 'zu'
626         
627       CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
628          IF ( av == 0 )  THEN
629             to_be_resorted => vv
630          ELSE
631             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
632                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
633                vv_av = REAL( fill_value, KIND = wp )
634             ENDIF
635             to_be_resorted => vv_av
636          ENDIF
637          flag_nr = 2
638         
639          IF ( mode == 'xy' )  grid = 'zu'
640               
641       CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
642          IF ( av == 0 )  THEN
643             to_be_resorted => ww
644          ELSE
645             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
646                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
647                ww_av = REAL( fill_value, KIND = wp )
648             ENDIF
649             to_be_resorted => ww_av
650          ENDIF
651          flag_nr = 3
652         
653          IF ( mode == 'xy' )  grid = 'zw'
654         
655       CASE ( 'theta_2m*_xy' )        ! 2d-array
656          IF ( av == 0 )  THEN
657             DO  i = nxl, nxr
658                DO  j = nys, nyn
659                   local_pf(i,j,nzb+1) = pt_2m(j,i)
660                ENDDO
661             ENDDO
662          ELSE
663             IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
664                ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
665                pt_2m_av = REAL( fill_value, KIND = wp )
666             ENDIF
667             DO  i = nxl, nxr
668                DO  j = nys, nyn
669                   local_pf(i,j,nzb+1) = pt_2m_av(j,i)
670                ENDDO
671             ENDDO
672          ENDIF
673          resorted = .TRUE.
674          two_d    = .TRUE.
675          grid     = 'zu1'
676         
677       CASE ( 'wspeed_10m*_xy' )        ! 2d-array
678          IF ( av == 0 )  THEN
679             DO  i = nxl, nxr
680                DO  j = nys, nyn
681                   local_pf(i,j,nzb+1) = uv_10m(j,i)
682                ENDDO
683             ENDDO
684          ELSE
685             IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
686                ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
687                uv_10m_av = REAL( fill_value, KIND = wp )
688             ENDIF
689             DO  i = nxl, nxr
690                DO  j = nys, nyn
691                   local_pf(i,j,nzb+1) = uv_10m_av(j,i)
692                ENDDO
693             ENDDO
694          ENDIF
695          resorted = .TRUE.
696          two_d    = .TRUE.
697          grid     = 'zu1'
698
699       CASE DEFAULT
700          found = .FALSE.
701          grid  = 'none'
702
703    END SELECT
704   
705    IF ( found  .AND.  .NOT. resorted )  THEN     
706       DO  i = nxl, nxr
707          DO  j = nys, nyn
708             DO  k = nzb_do, nzt_do
709                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
710                                     REAL( fill_value, KIND = wp ),            &
711                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 
712             ENDDO
713          ENDDO
714       ENDDO
715    ENDIF
716 
717 END SUBROUTINE doq_output_2d
718 
719 
720!------------------------------------------------------------------------------!
721!
722! Description:
723! ------------
724!> Subroutine defining 3D output variables
725!------------------------------------------------------------------------------!
726 SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do,  &
727                           nzt_do )
728 
729    IMPLICIT NONE
730
731    CHARACTER (LEN=*) ::  variable !<
732
733    INTEGER(iwp) ::  av       !< index indicating averaged or instantaneous output
734    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
735    INTEGER(iwp) ::  i        !< index variable along x-direction
736    INTEGER(iwp) ::  j        !< index variable along y-direction
737    INTEGER(iwp) ::  k        !< index variable along z-direction
738    INTEGER(iwp) ::  nzb_do   !< lower limit of the data output (usually 0)
739    INTEGER(iwp) ::  nzt_do   !< vertical upper limit of the data output (usually nz_do3d)
740
741    LOGICAL ::  found             !< true if variable is in list
742    LOGICAL ::  resorted          !< true if array is resorted
743
744    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
745
746    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf        !<
747    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
748
749    flag_nr  = 0
750    found    = .TRUE.
751    resorted = .FALSE.
752   
753    SELECT CASE ( TRIM( variable ) )
754
755       CASE ( 'ti' )
756          IF ( av == 0 )  THEN
757             to_be_resorted => ti
758          ELSE
759             IF ( .NOT. ALLOCATED( ti_av ) ) THEN
760                ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
761                ti_av = REAL( fill_value, KIND = wp )
762             ENDIF
763             to_be_resorted => ti_av
764          ENDIF
765          flag_nr = 0
766   
767       CASE ( 'uu' )
768          IF ( av == 0 )  THEN
769             to_be_resorted => uu
770          ELSE
771             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
772                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
773                uu_av = REAL( fill_value, KIND = wp )
774             ENDIF
775             to_be_resorted => uu_av
776          ENDIF
777          flag_nr = 1
778             
779       CASE ( 'vv' )
780          IF ( av == 0 )  THEN
781             to_be_resorted => vv
782          ELSE
783             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
784                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
785                vv_av = REAL( fill_value, KIND = wp )
786             ENDIF
787             to_be_resorted => vv_av
788          ENDIF
789          flag_nr = 2
790             
791       CASE ( 'ww' )
792          IF ( av == 0 )  THEN
793             to_be_resorted => ww
794          ELSE
795             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
796                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
797                ww_av = REAL( fill_value, KIND = wp )
798             ENDIF
799             to_be_resorted => ww_av
800          ENDIF
801          flag_nr = 3
802
803       CASE DEFAULT
804          found = .FALSE.
805
806    END SELECT
807   
808    IF ( found  .AND.  .NOT. resorted )  THEN     
809       DO  i = nxl, nxr
810          DO  j = nys, nyn
811             DO  k = nzb_do, nzt_do
812                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
813                                     REAL( fill_value, KIND = wp ),            &
814                                     BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) 
815             ENDDO
816          ENDDO
817       ENDDO
818    ENDIF
819
820 END SUBROUTINE doq_output_3d
821 
822! Description:
823! ------------
824!> Resorts the user-defined output quantity with indices (k,j,i) to a
825!> temporary array with indices (i,j,k) for masked data output.
826!------------------------------------------------------------------------------!
827 SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
828 
829    USE control_parameters
830       
831    USE indices
832
833    IMPLICIT NONE
834
835    CHARACTER (LEN=*) ::  variable   !<
836    CHARACTER (LEN=5) ::  grid       !< flag to distinquish between staggered grids
837
838    INTEGER(iwp) ::  av              !< index indicating averaged or instantaneous output
839    INTEGER(iwp) ::  flag_nr         !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
840    INTEGER(iwp) ::  i               !< index variable along x-direction
841    INTEGER(iwp) ::  j               !< index variable along y-direction
842    INTEGER(iwp) ::  k               !< index variable along z-direction
843    INTEGER(iwp) ::  im              !< loop index for masked variables
844    INTEGER(iwp) ::  jm              !< loop index for masked variables
845    INTEGER(iwp) ::  kk              !< masked output index variable along z-direction
846    INTEGER(iwp) ::  mid             !< masked output running index
847    INTEGER(iwp) ::  ktt             !< k index of highest horizontal surface
848
849    LOGICAL      ::  found           !< true if variable is in list
850    LOGICAL      ::  resorted        !< true if array is resorted
851
852    REAL(wp),                                                                  &
853       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
854          local_pf   !<
855    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
856
857    REAL(wp), PARAMETER   ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
858
859    flag_nr  = 0
860    found    = .TRUE.
861    resorted = .FALSE.
862    grid     = 's'
863
864    SELECT CASE ( TRIM( variable ) )
865
866       CASE ( 'ti' )
867          IF ( av == 0 )  THEN
868             to_be_resorted => ti
869          ELSE
870             to_be_resorted => ti_av
871          ENDIF
872          grid = 's'
873          flag_nr = 0
874   
875       CASE ( 'uu' )
876          IF ( av == 0 )  THEN
877             to_be_resorted => uu
878          ELSE
879             to_be_resorted => uu_av
880          ENDIF
881          grid = 'u'
882          flag_nr = 1
883   
884       CASE ( 'vv' )
885          IF ( av == 0 )  THEN
886             to_be_resorted => vv
887          ELSE
888             to_be_resorted => vv_av
889          ENDIF
890          grid = 'v'
891          flag_nr = 2
892   
893       CASE ( 'ww' )
894          IF ( av == 0 )  THEN
895             to_be_resorted => ww
896          ELSE
897             to_be_resorted => ww_av
898          ENDIF
899          grid = 'w'
900          flag_nr = 3
901
902       CASE DEFAULT
903          found = .FALSE.
904
905    END SELECT
906   
907    IF ( found  .AND.  .NOT. resorted )  THEN
908       IF ( .NOT. mask_surface(mid) )  THEN
909!
910!--       Default masked output
911          DO  i = 1, mask_size_l(mid,1)
912             DO  j = 1, mask_size_l(mid,2)
913                DO  k = 1, mask_size_l(mid,3)
914                   local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k),            &
915                                                     mask_j(mid,j),            &
916                                                     mask_i(mid,i))
917                ENDDO
918             ENDDO
919          ENDDO
920
921       ELSE
922!
923!--       Terrain-following masked output
924          DO  i = 1, mask_size_l(mid,1)
925             DO  j = 1, mask_size_l(mid,2)
926!
927!--             Get k index of the highest terraing surface
928                im = mask_i(mid,i)
929                jm = mask_j(mid,j)
930                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
931                              DIM = 1 ) - 1
932                DO  k = 1, mask_size_l(mid,3)
933                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
934!
935!--                Set value if not in building
936                   IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
937                      local_pf(i,j,k) = fill_value
938                   ELSE
939                      local_pf(i,j,k) = to_be_resorted(kk,jm,im)
940                   ENDIF
941                ENDDO
942             ENDDO
943          ENDDO
944
945       ENDIF
946    ENDIF
947   
948 END SUBROUTINE doq_output_mask
949
950!------------------------------------------------------------------------------!
951! Description:
952! ------------
953!> Allocate required arrays
954!------------------------------------------------------------------------------!
955 SUBROUTINE doq_init
956
957    IMPLICIT NONE
958   
959    INTEGER(iwp) ::  ivar   !< loop index over all 2d/3d/mask output quantities
960
961!
962!-- Next line is to avoid compiler warnings about unused variables
963    IF ( timestep_number_at_prev_calc == 0 )  CONTINUE
964!
965!-- Preparatory steps and initialization of output arrays
966    IF ( .NOT.  prepared_diagnostic_output_quantities )  CALL doq_prepare
967
968    initialized_diagnostic_output_quantities = .FALSE.
969   
970    ivar = 1
971   
972    DO  WHILE ( ivar <= SIZE( do_all ) ) 
973
974       SELECT CASE ( TRIM( do_all(ivar) ) )
975!
976!--       Allocate array for 'turbulence intensity'
977          CASE ( 'ti' )
978             IF ( .NOT. ALLOCATED( ti ) )  THEN
979                ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
980                ti = 0.0_wp
981             ENDIF
982!
983!--       Allocate array for uu
984          CASE ( 'uu' )
985             IF ( .NOT. ALLOCATED( uu ) )  THEN
986                ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) )
987                uu = 0.0_wp
988             ENDIF
989!
990!--       Allocate array for vv
991          CASE ( 'vv' )
992             IF ( .NOT. ALLOCATED( vv ) )  THEN
993                ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) )
994                vv = 0.0_wp
995             ENDIF
996!
997!--       Allocate array for ww
998          CASE ( 'ww' )
999             IF ( .NOT. ALLOCATED( ww ) )  THEN
1000                ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
1001                ww = 0.0_wp
1002             ENDIF
1003!
1004!--       Allocate array for 2-m potential temperature
1005          CASE ( 'theta_2m*' )
1006             IF ( .NOT. ALLOCATED( pt_2m ) )  THEN
1007                ALLOCATE( pt_2m(nys:nyn,nxl:nxr) )
1008                pt_2m = 0.0_wp
1009             ENDIF
1010!
1011!--       Allocate array for 10-m wind speed
1012          CASE ( 'wspeed_10m*' )
1013             IF ( .NOT. ALLOCATED( uv_10m ) )  THEN
1014                ALLOCATE( uv_10m(nys:nyn,nxl:nxr) )
1015                uv_10m = 0.0_wp
1016             ENDIF
1017
1018       END SELECT
1019
1020       ivar = ivar + 1
1021    ENDDO
1022
1023    initialized_diagnostic_output_quantities = .TRUE.
1024   
1025 END SUBROUTINE doq_init
1026
1027
1028!--------------------------------------------------------------------------------------------------!
1029! Description:
1030! ------------
1031!> Calculate diagnostic quantities
1032!--------------------------------------------------------------------------------------------------!
1033 SUBROUTINE doq_calculate
1034
1035    IMPLICIT NONE
1036
1037    INTEGER(iwp) ::  i          !< grid index x-dimension
1038    INTEGER(iwp) ::  j          !< grid index y-dimension
1039    INTEGER(iwp) ::  k          !< grid index z-dimension
1040    INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
1041   
1042    TYPE(surf_type), POINTER ::  surf     !< surf-type array, used to generalize subroutines
1043
1044
1045!     CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
1046
1047!
1048!-- Save timestep number to check in time_integration if doq_calculate
1049!-- has been called already, since the CALL occurs at two locations, but the calculations need to be
1050!-- done only once per timestep.
1051    timestep_number_at_prev_calc = current_timestep_number
1052
1053    ivar = 1
1054
1055    DO  WHILE ( ivar <= SIZE( do_all ) ) 
1056
1057       SELECT CASE ( TRIM( do_all(ivar) ) )
1058!
1059!--       Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
1060          CASE ( 'ti' )
1061             DO  i = nxl, nxr
1062                DO  j = nys, nyn
1063                   DO  k = nzb+1, nzt
1064                      ti(k,j,i) = 0.25_wp * SQRT(                              &
1065                        (   (   w(k,j+1,i) + w(k-1,j+1,i)                      &
1066                              - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy              &
1067                          - (   v(k+1,j,i) + v(k+1,j+1,i)                      &
1068                              - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2     &
1069                      + (   (   u(k+1,j,i) + u(k+1,j,i+1)                      &
1070                              - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)          &
1071                          - (   w(k,j,i+1) + w(k-1,j,i+1)                      &
1072                              - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2     &
1073                      + (   (   v(k,j,i+1) + v(k,j+1,i+1)                      &
1074                              - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx              &
1075                          - (   u(k,j+1,i) + u(k,j+1,i+1)                      &
1076                              - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )  &
1077                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0) )
1078                   ENDDO
1079                ENDDO
1080             ENDDO           
1081!
1082!--       uu
1083          CASE ( 'uu' )
1084             DO  i = nxl, nxr
1085                DO  j = nys, nyn
1086                   DO  k = nzb+1, nzt
1087                      uu(k,j,i) = u(k,j,i) * u(k,j,i)                          &
1088                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1) )
1089                   ENDDO
1090                ENDDO
1091             ENDDO
1092!
1093!--       vv
1094          CASE ( 'vv' )
1095             DO  i = nxl, nxr
1096                DO  j = nys, nyn
1097                   DO  k = nzb+1, nzt
1098                      vv(k,j,i) = v(k,j,i) * v(k,j,i)                          &
1099                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2) )
1100                   ENDDO
1101                ENDDO
1102             ENDDO
1103!
1104!--       ww
1105          CASE ( 'ww' )
1106             DO  i = nxl, nxr
1107                DO  j = nys, nyn
1108                   DO  k = nzb+1, nzt-1
1109                      ww(k,j,i) = w(k,j,i) * w(k,j,i)                          &
1110                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3) )
1111                   ENDDO
1112                ENDDO
1113             ENDDO
1114!
1115!--       2-m potential temperature
1116          CASE ( 'theta_2m*' )
1117!
1118!--          2-m potential temperature is caluclated from surface arrays. In
1119!--          case the 2m level is below the first grid point, MOST is applied,
1120!--          else, linear interpolation between two vertical grid levels is
1121!--          applied. To access all surfaces, iterate over all horizontally-
1122!--          upward facing surface types.
1123             surf => surf_def_h(0)
1124             CALL calc_pt_2m
1125             surf => surf_lsm_h
1126             CALL calc_pt_2m
1127             surf => surf_usm_h
1128             CALL calc_pt_2m
1129!
1130!--       10-m wind speed
1131          CASE ( 'wspeed_10m*' )
1132!
1133!--          10-m wind speed is caluclated from surface arrays. In
1134!--          case the 10m level is below the first grid point, MOST is applied,
1135!--          else, linear interpolation between two vertical grid levels is
1136!--          applied. To access all surfaces, iterate over all horizontally-
1137!--          upward facing surface types.
1138             surf => surf_def_h(0)
1139             CALL calc_wind_10m
1140             surf => surf_lsm_h
1141             CALL calc_wind_10m
1142             surf => surf_usm_h
1143             CALL calc_wind_10m
1144
1145       END SELECT
1146
1147       ivar = ivar + 1
1148    ENDDO
1149
1150!     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
1151
1152!
1153!-- The following block contains subroutines to calculate diagnostic
1154!-- surface quantities.
1155    CONTAINS
1156!------------------------------------------------------------------------------!
1157! Description:
1158! ------------
1159!> Calculation of 2-m potential temperature.
1160!------------------------------------------------------------------------------!
1161       SUBROUTINE calc_pt_2m
1162
1163          USE surface_layer_fluxes_mod,                                        &
1164              ONLY:  psi_h
1165
1166          IMPLICIT NONE
1167
1168          INTEGER(iwp) ::  kk     !< running index along the z-dimension
1169          INTEGER(iwp) ::  m      !< running index for surface elements
1170 
1171          DO  m = 1, surf%ns
1172
1173             i = surf%i(m)
1174             j = surf%j(m)
1175             k = surf%k(m)
1176!
1177!--          If 2-m level is below the first grid level, MOST is
1178!--          used for calculation of 2-m temperature.
1179             IF ( surf%z_mo(m) > 2.0_wp )  THEN
1180                pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa           &
1181                                * ( LOG( 2.0_wp /  surf%z0h(m) )               &
1182                                  - psi_h( 2.0_wp / surf%ol(m) )               &
1183                                  + psi_h( surf%z0h(m) / surf%ol(m) ) )
1184!
1185!--          If 2-m level is above the first grid level, 2-m temperature
1186!--          is linearly interpolated between the two nearest vertical grid
1187!--          levels. Note, since 2-m temperature is only computed for
1188!--          horizontal upward-facing surfaces, only a vertical
1189!--          interpolation is necessary.
1190             ELSE
1191!
1192!--             zw(k-1) defines the height of the surface.
1193                kk = k
1194                DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp  .AND.  kk <= nzt )
1195                   kk = kk + 1 
1196                ENDDO
1197!
1198!--             kk defines the index of the first grid level >= 2m.
1199                pt_2m(j,i) = pt(kk-1,j,i) +                                    &
1200                              ( zw(k-1) + 2.0_wp - zu(kk-1)     ) *            &
1201                              ( pt(kk,j,i)       - pt(kk-1,j,i) ) /            &
1202                              ( zu(kk)           - zu(kk-1)     )
1203             ENDIF
1204
1205          ENDDO
1206
1207       END SUBROUTINE calc_pt_2m
1208       
1209!------------------------------------------------------------------------------!
1210! Description:
1211! ------------
1212!> Calculation of 10-m wind speed.
1213!------------------------------------------------------------------------------!
1214       SUBROUTINE calc_wind_10m
1215
1216          USE surface_layer_fluxes_mod,                                        &
1217              ONLY:  psi_m
1218
1219          IMPLICIT NONE
1220
1221          INTEGER(iwp) ::  kk     !< running index along the z-dimension
1222          INTEGER(iwp) ::  m      !< running index for surface elements
1223
1224          REAL(wp) ::  uv_l !< wind speed at lower grid point
1225          REAL(wp) ::  uv_u !< wind speed at upper grid point
1226 
1227          DO  m = 1, surf%ns
1228
1229             i = surf%i(m)
1230             j = surf%j(m)
1231             k = surf%k(m)
1232!
1233!--          If 10-m level is below the first grid level, MOST is
1234!--          used for calculation of 10-m temperature.
1235             IF ( surf%z_mo(m) > 10.0_wp )  THEN
1236                uv_10m(j,i) = surf%us(m) / kappa                               &
1237                          * ( LOG( 10.0_wp /  surf%z0(m) )                     &
1238                              - psi_m( 10.0_wp    / surf%ol(m) )               &
1239                              + psi_m( surf%z0(m) / surf%ol(m) ) )
1240!
1241!--          If 10-m level is above the first grid level, 10-m wind speed
1242!--          is linearly interpolated between the two nearest vertical grid
1243!--          levels. Note, since 10-m temperature is only computed for
1244!--          horizontal upward-facing surfaces, only a vertical
1245!--          interpolation is necessary.
1246             ELSE
1247!
1248!--             zw(k-1) defines the height of the surface.
1249                kk = k
1250                DO WHILE ( zu(kk) - zw(k-1) < 10.0_wp  .AND.  kk <= nzt )
1251                   kk = kk + 1 
1252                ENDDO
1253!
1254!--             kk defines the index of the first grid level >= 10m.
1255                uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2   &
1256                           + ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
1257
1258                uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i)   + u(kk,j,i+1)   ) )**2   &
1259                           + ( 0.5_wp * ( v(kk,j,i)   + v(kk,j+1,i)   ) )**2 )
1260
1261                uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) *        &
1262                                     ( uv_u              - uv_l     ) /        &
1263                                     ( zu(kk)            - zu(kk-1) )
1264
1265             ENDIF
1266
1267          ENDDO
1268
1269       END SUBROUTINE calc_wind_10m
1270
1271 END SUBROUTINE doq_calculate
1272
1273
1274!------------------------------------------------------------------------------!
1275! Description:
1276! ------------
1277!> Preparation of the diagnostic output, counting of the module-specific
1278!> output quantities and gathering of the output names.
1279!------------------------------------------------------------------------------!
1280 SUBROUTINE doq_prepare
1281
1282
1283    USE control_parameters,                                                    &
1284        ONLY:  do2d, do3d, domask, masks
1285
1286    IMPLICIT NONE
1287
1288    CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d_var = ' '  !<
1289                                                          !< label array for 2d output quantities
1290
1291    INTEGER(iwp) ::  av         !< index defining type of output, av=0 instantaneous, av=1 averaged
1292    INTEGER(iwp) ::  ivar       !< loop index
1293    INTEGER(iwp) ::  ivar_all   !< loop index
1294    INTEGER(iwp) ::  l          !< index for cutting string
1295    INTEGER(iwp) ::  mid          !< masked output running index
1296
1297    prepared_diagnostic_output_quantities = .FALSE.
1298
1299    ivar     = 1
1300    ivar_all = 1
1301
1302    DO  av = 0, 1
1303!
1304!--    Remove _xy, _xz, or _yz from string
1305       l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1306       do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1307!
1308!--    Gather 2d output quantity names.
1309!--    Check for double occurrence of output quantity, e.g. by _xy,
1310!--    _yz, _xz.
1311       DO  WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
1312          IF ( .NOT.  ANY( do_all == do2d_var(av,ivar) ) )  THEN
1313             do_all(ivar_all) = do2d_var(av,ivar)
1314          ENDIF
1315          ivar = ivar + 1
1316          ivar_all = ivar_all + 1
1317          l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1318          do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1319       ENDDO
1320
1321       ivar = 1
1322!
1323!--    Gather 3d output quantity names
1324       DO  WHILE ( do3d(av,ivar)(1:1) /= ' ' )
1325          do_all(ivar_all) = do3d(av,ivar)
1326          ivar = ivar + 1
1327          ivar_all = ivar_all + 1
1328       ENDDO
1329
1330       ivar = 1
1331!
1332!--    Gather masked output quantity names. Also check for double output
1333!--    e.g. by different masks.
1334       DO  mid = 1, masks
1335          DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
1336             IF ( .NOT.  ANY( do_all == domask(mid,av,ivar) ) )  THEN
1337                do_all(ivar_all) = domask(mid,av,ivar)
1338             ENDIF
1339
1340             ivar = ivar + 1
1341             ivar_all = ivar_all + 1
1342          ENDDO
1343          ivar = 1
1344       ENDDO
1345
1346    ENDDO
1347
1348    prepared_diagnostic_output_quantities = .TRUE.
1349
1350 END SUBROUTINE doq_prepare
1351 
1352!------------------------------------------------------------------------------!
1353! Description:
1354! ------------
1355!> Subroutine reads local (subdomain) restart data
1356!> Note: With the current structure reading of non-standard array is not
1357!> possible
1358!------------------------------------------------------------------------------!
1359!  SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
1360!                            nxr_on_file, nynf, nync, nyn_on_file, nysf,         &
1361!                            nysc, nys_on_file, tmp_3d_non_standard, found )
1362
1363!
1364!     USE control_parameters
1365!         
1366!     USE indices
1367!     
1368!     USE kinds
1369!     
1370!     USE pegrid
1371!
1372!
1373!     IMPLICIT NONE
1374!
1375!     INTEGER(iwp) ::  k               !<
1376!     INTEGER(iwp) ::  nxlc            !<
1377!     INTEGER(iwp) ::  nxlf            !<
1378!     INTEGER(iwp) ::  nxl_on_file     !<
1379!     INTEGER(iwp) ::  nxrc            !<
1380!     INTEGER(iwp) ::  nxrf            !<
1381!     INTEGER(iwp) ::  nxr_on_file     !<
1382!     INTEGER(iwp) ::  nync            !<
1383!     INTEGER(iwp) ::  nynf            !<
1384!     INTEGER(iwp) ::  nyn_on_file     !<
1385!     INTEGER(iwp) ::  nysc            !<
1386!     INTEGER(iwp) ::  nysf            !<
1387!     INTEGER(iwp) ::  nys_on_file     !<
1388!
1389!     LOGICAL, INTENT(OUT)  :: found
1390!     
1391!     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions
1392! !
1393! !-- If temporary non-standard array for reading is already allocated,
1394! !-- deallocate it.
1395!     IF ( ALLOCATED( tmp_3d_non_standard ) )  DEALLOCATE( tmp_3d_non_standard )
1396!     
1397!     found = .TRUE.
1398!
1399!     SELECT CASE ( restart_string(1:length) )
1400!
1401!        CASE ( 'ti_av' )
1402!           IF ( .NOT. ALLOCATED( ti_av ) )  THEN
1403!              ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1404!           ENDIF
1405!           IF ( k == 1 )  THEN
1406!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
1407!                                            nxl_on_file:nxr_on_file) )
1408!              READ ( 13 )  tmp_3d_non_standard
1409!           ENDIF
1410!           ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1411!     
1412!        CASE ( 'uu_av' )
1413!           IF ( .NOT. ALLOCATED( uu_av ) )  THEN
1414!              ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1415!           ENDIF
1416!           IF ( k == 1 )  THEN
1417!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
1418!                                            nxl_on_file:nxr_on_file) )
1419!              READ ( 13 )  tmp_3d_non_standard
1420!           ENDIF
1421!           uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1422!                   
1423!        CASE ( 'vv_av' )
1424!           IF ( .NOT. ALLOCATED( vv_av ) )  THEN
1425!              ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1426!           ENDIF
1427!           IF ( k == 1 )  THEN
1428!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
1429!                                            nxl_on_file:nxr_on_file) )
1430!              READ ( 13 )  tmp_3d_non_standard
1431!           ENDIF
1432!           vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1433!                   
1434!        CASE ( 'ww_av' )
1435!           IF ( .NOT. ALLOCATED( ww_av ) )  THEN
1436!              ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1437!           ENDIF
1438!           IF ( k == 1 )  THEN
1439!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
1440!                                            nxl_on_file:nxr_on_file) )
1441!              READ ( 13 )  tmp_3d_non_standard
1442!           ENDIF
1443!           ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1444!                         
1445!
1446!        CASE DEFAULT
1447!
1448!           found = .FALSE.
1449!
1450!     END SELECT
1451!
1452!  END SUBROUTINE doq_rrd_local
1453 
1454!------------------------------------------------------------------------------!
1455! Description:
1456! ------------
1457!> Subroutine writes local (subdomain) restart data
1458!------------------------------------------------------------------------------!
1459 SUBROUTINE doq_wrd_local
1460
1461
1462    IMPLICIT NONE
1463
1464    IF ( ALLOCATED( pt_2m_av ) )  THEN
1465       CALL wrd_write_string( 'pt_2m_av' )
1466       WRITE ( 14 )  pt_2m_av
1467    ENDIF
1468
1469    IF ( ALLOCATED( ti_av ) )  THEN
1470       CALL wrd_write_string( 'ti_av' )
1471       WRITE ( 14 )  ti_av
1472    ENDIF
1473
1474    IF ( ALLOCATED( uu_av ) )  THEN
1475       CALL wrd_write_string( 'uu_av' )
1476       WRITE ( 14 )  uu_av
1477    ENDIF
1478
1479    IF ( ALLOCATED( uv_10m_av ) )  THEN
1480       CALL wrd_write_string( 'uv_10m_av' )
1481       WRITE ( 14 )  uv_10m_av
1482    ENDIF
1483
1484    IF ( ALLOCATED( vv_av ) )  THEN
1485       CALL wrd_write_string( 'vv_av' )
1486       WRITE ( 14 )  vv_av
1487    ENDIF
1488
1489    IF ( ALLOCATED( ww_av ) )  THEN
1490       CALL wrd_write_string( 'ww_av' )
1491       WRITE ( 14 )  ww_av
1492    ENDIF
1493
1494
1495 END SUBROUTINE doq_wrd_local
1496 
1497 
1498
1499 END MODULE diagnostic_output_quantities_mod
Note: See TracBrowser for help on using the repository browser.