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

Last change on this file since 4381 was 4360, checked in by suehring, 4 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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