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

Last change on this file since 4351 was 4351, checked in by oliver.maas, 16 months ago

added output of wu, wv, wtheta and wq to enable covariance calculation according to temporal EC method

  • 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diagnostic_output_quantities_mod.f90 4351 2019-12-19 10:43:19Z oliver.maas $
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.