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

Last change on this file since 4331 was 4331, checked in by suehring, 2 years ago

New diagnostic output for 10-m wind speed; Diagnostic output of 2-m potential temperature moved to diagnostic output

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