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

Last change on this file since 4172 was 4167, checked in by suehring, 5 years ago

Merge from branch resler: Changed behaviour of masked output over surface to follow terrain and ignore buildings

  • Property svn:keywords set to Id
File size: 38.4 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 4167 2019-08-16 11:01:48Z oliver.maas $
27! Changed behaviour of masked output over surface to follow terrain and ignore
28! buildings (J.Resler, T.Gronemeier)
29!
30! 4157 2019-08-14 09:19:12Z suehring
31! Initialization restructured, in order to work also when data output during
32! spin-up is enabled.
33!
34! 4132 2019-08-02 12:34:17Z suehring
35! Bugfix in masked data output
36!
37! 4069 2019-07-01 14:05:51Z Giersch
38! Masked output running index mid has been introduced as a local variable to
39! avoid runtime error (Loop variable has been modified) in time_integration
40!
41! 4039 2019-06-18 10:32:41Z suehring
42! - Add output of uu, vv, ww to enable variance calculation according temporal
43!   EC method
44! - Allocate arrays only when they are required
45! - Formatting adjustment
46! - Rename subroutines
47! - Further modularization
48!
49! 3998 2019-05-23 13:38:11Z suehring
50! Bugfix in gathering all output strings
51!
52! 3995 2019-05-22 18:59:54Z suehring
53! Avoid compiler warnings about unused variable and fix string operation which
54! is not allowed with PGI compiler
55!
56! 3994 2019-05-22 18:08:09Z suehring
57! Initial revision
58!
59!
60! @author Farah Kanani-Suehring
61!
62! Description:
63! ------------
64!> ...
65!------------------------------------------------------------------------------!
66 MODULE diagnostic_output_quantities_mod
67 
68    USE arrays_3d,                                                             &
69        ONLY:  ddzu, u, v, w
70!
71!     USE averaging
72!
73!     USE basic_constants_and_equations_mod,                                     &
74!         ONLY:  c_p, lv_d_cp, l_v
75!
76!     USE bulk_cloud_model_mod,                                                  &
77!         ONLY:  bulk_cloud_model
78!
79    USE control_parameters,                                                    &
80        ONLY:  current_timestep_number, varnamelength
81!
82!     USE cpulog,                                                                &
83!         ONLY:  cpu_log, log_point
84
85   USE grid_variables,                                                         &
86        ONLY:  ddx, ddy
87!
88    USE indices,                                                               &
89        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
90!
91    USE kinds
92!
93!     USE land_surface_model_mod,                                                &
94!         ONLY:  zs
95!
96!     USE module_interface,                                                      &
97!         ONLY:  module_interface_data_output_2d
98!
99! #if defined( __netcdf )
100!     USE NETCDF
101! #endif
102!
103!     USE netcdf_interface,                                                      &
104!         ONLY:  fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d,       &
105!                id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat,        &
106!                netcdf_data_format, netcdf_handle_error
107!
108!     USE particle_attributes,                                                   &
109!         ONLY:  grid_particles, number_of_particles, particle_advection_start,  &
110!                particles, prt_count
111!     
112!     USE pegrid
113!
114!     USE surface_mod,                                                           &
115!         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
116!                surf_lsm_h, surf_usm_h
117!
118!     USE turbulence_closure_mod,                                                &
119!         ONLY:  tcm_data_output_2d
120
121
122    IMPLICIT NONE
123
124    CHARACTER(LEN=varnamelength), DIMENSION(500) ::  do_all = ' '
125
126    INTEGER(iwp) ::  timestep_number_at_prev_calc = 0  !< ...at previous diagnostic output calculation
127 
128    LOGICAL ::  initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
129    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.    !< flag indicating whether output is p
130
131    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti     !< rotation(u,v,w) aka turbulence intensity
132    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ti_av  !< avg. rotation(u,v,w) aka turbulence intensity
133   
134    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu     !< uu
135    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  uu_av  !< mean of uu
136   
137    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv     !< vv
138    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  vv_av  !< mean of vv
139   
140    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww     !< ww
141    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ww_av  !< mean of ww
142
143
144    SAVE
145
146    PRIVATE
147
148!
149!-- Public variables
150    PUBLIC do_all,                                                             &
151           initialized_diagnostic_output_quantities,                           &
152           prepared_diagnostic_output_quantities,                              &
153           timestep_number_at_prev_calc,                                       &
154           ti_av,                                                              &
155           uu_av,                                                              &
156           vv_av,                                                              &
157           ww_av
158!                                                                             
159!-- Public routines                                                           
160    PUBLIC doq_3d_data_averaging,                                              &
161           doq_calculate,                                                      &
162           doq_check_data_output,                                              &
163           doq_define_netcdf_grid,                                             &
164           doq_output_2d,                                                      &
165           doq_output_3d,                                                      &
166           doq_output_mask,                                                    &
167           doq_init,                                                           &
168           doq_wrd_local
169!          doq_rrd_local,                                                      &
170
171
172    INTERFACE doq_3d_data_averaging
173       MODULE PROCEDURE doq_3d_data_averaging
174    END INTERFACE doq_3d_data_averaging       
175
176    INTERFACE doq_calculate
177       MODULE PROCEDURE doq_calculate
178    END INTERFACE doq_calculate
179
180    INTERFACE doq_check_data_output
181       MODULE PROCEDURE doq_check_data_output
182    END INTERFACE doq_check_data_output
183   
184    INTERFACE doq_define_netcdf_grid
185       MODULE PROCEDURE doq_define_netcdf_grid
186    END INTERFACE doq_define_netcdf_grid
187   
188    INTERFACE doq_output_2d
189       MODULE PROCEDURE doq_output_2d
190    END INTERFACE doq_output_2d
191   
192    INTERFACE doq_output_3d
193       MODULE PROCEDURE doq_output_3d
194    END INTERFACE doq_output_3d
195   
196    INTERFACE doq_output_mask
197       MODULE PROCEDURE doq_output_mask
198    END INTERFACE doq_output_mask
199     
200    INTERFACE doq_init
201       MODULE PROCEDURE doq_init
202    END INTERFACE doq_init
203
204    INTERFACE doq_prepare
205       MODULE PROCEDURE doq_prepare
206    END INTERFACE doq_prepare
207   
208!     INTERFACE doq_rrd_local
209!        MODULE PROCEDURE doq_rrd_local
210!     END INTERFACE doq_rrd_local
211   
212    INTERFACE doq_wrd_local
213       MODULE PROCEDURE doq_wrd_local
214    END INTERFACE doq_wrd_local
215
216
217 CONTAINS
218 
219!------------------------------------------------------------------------------!
220! Description:
221! ------------
222!> Sum up and time-average diagnostic output quantities as well as allocate
223!> the array necessary for storing the average.
224!------------------------------------------------------------------------------!
225 SUBROUTINE doq_3d_data_averaging( mode, variable )
226
227    USE control_parameters,                                                    &
228        ONLY:  average_count_3d
229
230    CHARACTER (LEN=*) ::  mode     !<
231    CHARACTER (LEN=*) ::  variable !<
232
233    INTEGER(iwp) ::  i !<
234    INTEGER(iwp) ::  j !<
235    INTEGER(iwp) ::  k !<
236
237    IF ( mode == 'allocate' )  THEN
238
239       SELECT CASE ( TRIM( variable ) )
240
241          CASE ( 'ti' )
242             IF ( .NOT. ALLOCATED( ti_av ) )  THEN
243                ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
244             ENDIF
245             ti_av = 0.0_wp
246       
247          CASE ( 'uu' )
248             IF ( .NOT. ALLOCATED( uu_av ) )  THEN
249                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
250             ENDIF
251             uu_av = 0.0_wp
252               
253          CASE ( 'vv' )
254             IF ( .NOT. ALLOCATED( vv_av ) )  THEN
255                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
256             ENDIF
257             vv_av = 0.0_wp
258               
259          CASE ( 'ww' )
260             IF ( .NOT. ALLOCATED( ww_av ) )  THEN
261                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
262             ENDIF
263             ww_av = 0.0_wp
264
265          CASE DEFAULT
266             CONTINUE
267
268       END SELECT
269
270    ELSEIF ( mode == 'sum' )  THEN
271
272       SELECT CASE ( TRIM( variable ) )
273 
274          CASE ( 'ti' )
275             IF ( ALLOCATED( ti_av ) ) THEN
276                DO  i = nxl, nxr
277                   DO  j = nys, nyn
278                      DO  k = nzb, nzt+1
279                         ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
280                      ENDDO
281                   ENDDO
282                ENDDO
283             ENDIF
284             
285          CASE ( 'uu' )
286             IF ( ALLOCATED( uu_av ) ) THEN
287                DO  i = nxl, nxr
288                   DO  j = nys, nyn
289                      DO  k = nzb, nzt+1
290                         uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
291                      ENDDO
292                   ENDDO
293                ENDDO
294             ENDIF
295             
296          CASE ( 'vv' )
297             IF ( ALLOCATED( vv_av ) ) THEN
298                DO  i = nxl, nxr
299                   DO  j = nys, nyn
300                      DO  k = nzb, nzt+1
301                         vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
302                      ENDDO
303                   ENDDO
304                ENDDO
305             ENDIF
306             
307          CASE ( 'ww' )
308             IF ( ALLOCATED( ww_av ) ) THEN
309                DO  i = nxl, nxr
310                   DO  j = nys, nyn
311                      DO  k = nzb, nzt+1
312                         ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
313                      ENDDO
314                   ENDDO
315                ENDDO
316             ENDIF
317
318          CASE DEFAULT
319             CONTINUE
320
321       END SELECT
322
323    ELSEIF ( mode == 'average' )  THEN
324
325       SELECT CASE ( TRIM( variable ) )
326
327          CASE ( 'ti' )
328             IF ( ALLOCATED( ti_av ) ) THEN
329                DO  i = nxl, nxr
330                   DO  j = nys, nyn
331                      DO  k = nzb, nzt+1
332                         ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
333                      ENDDO
334                   ENDDO
335                ENDDO
336             ENDIF
337       
338          CASE ( 'uu' )
339             IF ( ALLOCATED( uu_av ) ) THEN
340                DO  i = nxl, nxr
341                   DO  j = nys, nyn
342                      DO  k = nzb, nzt+1
343                         uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
344                      ENDDO
345                   ENDDO
346                ENDDO
347             ENDIF
348             
349          CASE ( 'vv' )
350             IF ( ALLOCATED( vv_av ) ) THEN
351                DO  i = nxl, nxr
352                   DO  j = nys, nyn
353                      DO  k = nzb, nzt+1
354                         vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
355                      ENDDO
356                   ENDDO
357                ENDDO
358             ENDIF
359             
360          CASE ( 'ww' )
361             IF ( ALLOCATED( ww_av ) ) THEN
362                DO  i = nxl, nxr
363                   DO  j = nys, nyn
364                      DO  k = nzb, nzt+1
365                         ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
366                      ENDDO
367                   ENDDO
368                ENDDO
369             ENDIF
370
371       END SELECT
372
373    ENDIF
374
375
376 END SUBROUTINE doq_3d_data_averaging 
377 
378!------------------------------------------------------------------------------!
379! Description:
380! ------------
381!> Check data output for diagnostic output
382!------------------------------------------------------------------------------!
383 SUBROUTINE doq_check_data_output( var, unit )
384
385    IMPLICIT NONE
386
387    CHARACTER (LEN=*) ::  unit  !<
388    CHARACTER (LEN=*) ::  var   !<
389
390    SELECT CASE ( TRIM( var ) )
391
392       CASE ( 'ti' )
393          unit = '1/s'
394             
395       CASE ( 'uu' )
396          unit = 'm2/s2'
397             
398       CASE ( 'vv' )
399          unit = 'm2/s2'
400             
401       CASE ( 'ww' )
402          unit = 'm2/s2'
403             
404       CASE DEFAULT
405          unit = 'illegal'
406
407    END SELECT
408
409
410 END SUBROUTINE doq_check_data_output
411 
412!------------------------------------------------------------------------------!
413!
414! Description:
415! ------------
416!> Subroutine defining appropriate grid for netcdf variables.
417!------------------------------------------------------------------------------!
418 SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
419   
420    IMPLICIT NONE
421
422    CHARACTER (LEN=*), INTENT(IN)  ::  variable    !<
423    LOGICAL, INTENT(OUT)           ::  found       !<
424    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
425    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
426    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
427
428    found  = .TRUE.
429   
430    SELECT CASE ( TRIM( variable ) )
431!
432!--    s grid
433       CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz'  )
434
435          grid_x = 'x'
436          grid_y = 'y'
437          grid_z = 'zu'   
438!
439!--    u grid
440       CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
441
442          grid_x = 'xu'
443          grid_y = 'y'
444          grid_z = 'zu'
445!
446!--    v grid
447       CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz'  )
448
449          grid_x = 'x'
450          grid_y = 'yv'
451          grid_z = 'zu'
452!
453!--    w grid
454       CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz'  )
455
456          grid_x = 'x'
457          grid_y = 'y'
458          grid_z = 'zw'
459
460       CASE DEFAULT
461          found  = .FALSE.
462          grid_x = 'none'
463          grid_y = 'none'
464          grid_z = 'none'
465
466    END SELECT
467
468
469 END SUBROUTINE doq_define_netcdf_grid
470 
471!------------------------------------------------------------------------------!
472!
473! Description:
474! ------------
475!> Subroutine defining 2D output variables
476!------------------------------------------------------------------------------!
477 SUBROUTINE doq_output_2d( av, variable, found, grid,                          &
478                           mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
479
480
481    IMPLICIT NONE
482
483    CHARACTER (LEN=*) ::  grid     !<
484    CHARACTER (LEN=*) ::  mode     !<
485    CHARACTER (LEN=*) ::  variable !<
486
487    INTEGER(iwp) ::  av       !< value indicating averaged or non-averaged output
488    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
489    INTEGER(iwp) ::  i        !< grid index x-direction
490    INTEGER(iwp) ::  j        !< grid index y-direction
491    INTEGER(iwp) ::  k        !< grid index z-direction
492    INTEGER(iwp) ::  nzb_do   !<
493    INTEGER(iwp) ::  nzt_do   !<
494
495    LOGICAL ::  found             !< true if variable is in list
496    LOGICAL ::  resorted          !< true if array is resorted
497    LOGICAL ::  two_d             !< flag parameter that indicates 2D variables (horizontal cross sections)
498
499    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
500   
501    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
502    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
503   
504    flag_nr  = 0
505    found    = .TRUE.
506    resorted = .FALSE.
507    two_d    = .FALSE.
508
509    SELECT CASE ( TRIM( variable ) )
510
511       CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
512           IF ( av == 0 )  THEN
513              to_be_resorted => ti
514           ELSE
515              IF ( .NOT. ALLOCATED( ti_av ) ) THEN
516                 ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
517                 ti_av = REAL( fill_value, KIND = wp )
518              ENDIF
519              to_be_resorted => ti_av
520           ENDIF
521           flag_nr = 0
522           
523           IF ( mode == 'xy' )  grid = 'zu'
524   
525       CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
526          IF ( av == 0 )  THEN
527             to_be_resorted => uu
528          ELSE
529             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
530                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
531                uu_av = REAL( fill_value, KIND = wp )
532             ENDIF
533             to_be_resorted => uu_av
534          ENDIF
535          flag_nr = 1
536         
537          IF ( mode == 'xy' )  grid = 'zu'
538         
539       CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
540          IF ( av == 0 )  THEN
541             to_be_resorted => vv
542          ELSE
543             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
544                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
545                vv_av = REAL( fill_value, KIND = wp )
546             ENDIF
547             to_be_resorted => vv_av
548          ENDIF
549          flag_nr = 2
550         
551          IF ( mode == 'xy' )  grid = 'zu'
552               
553       CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
554          IF ( av == 0 )  THEN
555             to_be_resorted => ww
556          ELSE
557             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
558                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
559                ww_av = REAL( fill_value, KIND = wp )
560             ENDIF
561             to_be_resorted => ww_av
562          ENDIF
563          flag_nr = 3
564         
565          IF ( mode == 'xy' )  grid = 'zw'
566
567       CASE DEFAULT
568          found = .FALSE.
569          grid  = 'none'
570
571    END SELECT
572   
573    IF ( found  .AND.  .NOT. resorted )  THEN     
574       DO  i = nxl, nxr
575          DO  j = nys, nyn
576             DO  k = nzb_do, nzt_do
577                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
578                                         REAL( fill_value, KIND = wp ),        &
579                                         BTEST( wall_flags_0(k,j,i), flag_nr ) ) 
580             ENDDO
581          ENDDO
582       ENDDO
583    ENDIF
584 
585 END SUBROUTINE doq_output_2d
586 
587 
588!------------------------------------------------------------------------------!
589!
590! Description:
591! ------------
592!> Subroutine defining 3D output variables
593!------------------------------------------------------------------------------!
594 SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do,  &
595                           nzt_do )
596 
597    IMPLICIT NONE
598
599    CHARACTER (LEN=*) ::  variable !<
600
601    INTEGER(iwp) ::  av       !< index indicating averaged or instantaneous output
602    INTEGER(iwp) ::  flag_nr  !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
603    INTEGER(iwp) ::  i        !< index variable along x-direction
604    INTEGER(iwp) ::  j        !< index variable along y-direction
605    INTEGER(iwp) ::  k        !< index variable along z-direction
606    INTEGER(iwp) ::  nzb_do   !< lower limit of the data output (usually 0)
607    INTEGER(iwp) ::  nzt_do   !< vertical upper limit of the data output (usually nz_do3d)
608
609    LOGICAL ::  found             !< true if variable is in list
610    LOGICAL ::  resorted          !< true if array is resorted
611
612    REAL(wp) ::  fill_value       !< value for the _FillValue attribute
613
614    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf        !<
615    REAL(wp), DIMENSION(:,:,:), POINTER ::                 to_be_resorted  !< points to array which needs to be resorted for output
616
617    flag_nr  = 0
618    found    = .TRUE.
619    resorted = .FALSE.
620   
621    SELECT CASE ( TRIM( variable ) )
622
623       CASE ( 'ti' )
624          IF ( av == 0 )  THEN
625             to_be_resorted => ti
626          ELSE
627             IF ( .NOT. ALLOCATED( ti_av ) ) THEN
628                ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
629                ti_av = REAL( fill_value, KIND = wp )
630             ENDIF
631             to_be_resorted => ti_av
632          ENDIF
633          flag_nr = 0
634   
635       CASE ( 'uu' )
636          IF ( av == 0 )  THEN
637             to_be_resorted => uu
638          ELSE
639             IF ( .NOT. ALLOCATED( uu_av ) ) THEN
640                ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
641                uu_av = REAL( fill_value, KIND = wp )
642             ENDIF
643             to_be_resorted => uu_av
644          ENDIF
645          flag_nr = 1
646             
647       CASE ( 'vv' )
648          IF ( av == 0 )  THEN
649             to_be_resorted => vv
650          ELSE
651             IF ( .NOT. ALLOCATED( vv_av ) ) THEN
652                ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
653                vv_av = REAL( fill_value, KIND = wp )
654             ENDIF
655             to_be_resorted => vv_av
656          ENDIF
657          flag_nr = 2
658             
659       CASE ( 'ww' )
660          IF ( av == 0 )  THEN
661             to_be_resorted => ww
662          ELSE
663             IF ( .NOT. ALLOCATED( ww_av ) ) THEN
664                ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
665                ww_av = REAL( fill_value, KIND = wp )
666             ENDIF
667             to_be_resorted => ww_av
668          ENDIF
669          flag_nr = 3
670
671       CASE DEFAULT
672          found = .FALSE.
673
674    END SELECT
675   
676    IF ( found  .AND.  .NOT. resorted )  THEN     
677       DO  i = nxl, nxr
678          DO  j = nys, nyn
679             DO  k = nzb_do, nzt_do
680                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
681                                         REAL( fill_value, KIND = wp ),        &
682                                         BTEST( wall_flags_0(k,j,i), flag_nr ) ) 
683             ENDDO
684          ENDDO
685       ENDDO
686    ENDIF
687
688 END SUBROUTINE doq_output_3d
689 
690! Description:
691! ------------
692!> Resorts the user-defined output quantity with indices (k,j,i) to a
693!> temporary array with indices (i,j,k) for masked data output.
694!------------------------------------------------------------------------------!
695 SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
696 
697    USE control_parameters
698       
699    USE indices
700
701    IMPLICIT NONE
702
703    CHARACTER (LEN=*) ::  variable   !<
704    CHARACTER (LEN=5) ::  grid       !< flag to distinquish between staggered grids
705
706    INTEGER(iwp) ::  av              !< index indicating averaged or instantaneous output
707    INTEGER(iwp) ::  flag_nr         !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
708    INTEGER(iwp) ::  i               !< index variable along x-direction
709    INTEGER(iwp) ::  j               !< index variable along y-direction
710    INTEGER(iwp) ::  k               !< index variable along z-direction
711    INTEGER(iwp) ::  im              !< loop index for masked variables
712    INTEGER(iwp) ::  jm              !< loop index for masked variables
713    INTEGER(iwp) ::  kk              !< masked output index variable along z-direction
714    INTEGER(iwp) ::  mid             !< masked output running index
715    INTEGER(iwp) ::  ktt             !< k index of highest horizontal surface
716
717    LOGICAL      ::  found           !< true if variable is in list
718    LOGICAL      ::  resorted        !< true if array is resorted
719
720    REAL(wp),                                                                  &
721       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
722          local_pf   !<
723    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
724
725    REAL(wp), PARAMETER   ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
726
727    flag_nr  = 0
728    found    = .TRUE.
729    resorted = .FALSE.
730    grid     = 's'
731
732    SELECT CASE ( TRIM( variable ) )
733
734       CASE ( 'ti' )
735          IF ( av == 0 )  THEN
736             to_be_resorted => ti
737          ELSE
738             to_be_resorted => ti_av
739          ENDIF
740          grid = 's'
741          flag_nr = 0
742   
743       CASE ( 'uu' )
744          IF ( av == 0 )  THEN
745             to_be_resorted => uu
746          ELSE
747             to_be_resorted => uu_av
748          ENDIF
749          grid = 'u'
750          flag_nr = 1
751   
752       CASE ( 'vv' )
753          IF ( av == 0 )  THEN
754             to_be_resorted => vv
755          ELSE
756             to_be_resorted => vv_av
757          ENDIF
758          grid = 'v'
759          flag_nr = 2
760   
761       CASE ( 'ww' )
762          IF ( av == 0 )  THEN
763             to_be_resorted => ww
764          ELSE
765             to_be_resorted => ww_av
766          ENDIF
767          grid = 'w'
768          flag_nr = 3
769
770
771       CASE DEFAULT
772          found = .FALSE.
773
774    END SELECT
775   
776    IF ( found  .AND.  .NOT. resorted )  THEN
777       IF ( .NOT. mask_surface(mid) )  THEN
778!
779!--       Default masked output
780          DO  i = 1, mask_size_l(mid,1)
781             DO  j = 1, mask_size_l(mid,2)
782                DO  k = 1, mask_size_l(mid,3)
783                   local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k),            &
784                                                     mask_j(mid,j),            &
785                                                     mask_i(mid,i))
786                ENDDO
787             ENDDO
788          ENDDO
789
790       ELSE
791!
792!--       Terrain-following masked output
793          DO  i = 1, mask_size_l(mid,1)
794             DO  j = 1, mask_size_l(mid,2)
795!
796!--             Get k index of the highest terraing surface
797                im = mask_i(mid,i)
798                jm = mask_j(mid,j)
799                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), &
800                              DIM = 1 ) - 1
801                DO  k = 1, mask_size_l(mid,3)
802                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
803!
804!--                Set value if not in building
805                   IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
806                      local_pf(i,j,k) = fill_value
807                   ELSE
808                      local_pf(i,j,k) = to_be_resorted(kk,jm,im)
809                   ENDIF
810                ENDDO
811             ENDDO
812          ENDDO
813
814       ENDIF
815    ENDIF
816   
817 END SUBROUTINE doq_output_mask
818
819!------------------------------------------------------------------------------!
820! Description:
821! ------------
822!> Allocate required arrays
823!------------------------------------------------------------------------------!
824 SUBROUTINE doq_init
825
826    IMPLICIT NONE
827   
828    INTEGER(iwp) ::  ivar   !< loop index over all 2d/3d/mask output quantities
829
830!
831!-- Next line is to avoid compiler warnings about unused variables
832    IF ( timestep_number_at_prev_calc == 0 )  CONTINUE
833!
834!-- Preparatory steps and initialization of output arrays
835    IF ( .NOT.  prepared_diagnostic_output_quantities )  CALL doq_prepare
836
837    initialized_diagnostic_output_quantities = .FALSE.
838   
839    ivar = 1
840   
841    DO  WHILE ( ivar <= SIZE( do_all ) ) 
842
843       SELECT CASE ( TRIM( do_all(ivar) ) )
844!
845!--       Allocate array for 'turbulence intensity'
846          CASE ( 'ti' )
847             IF ( .NOT. ALLOCATED( ti ) )  THEN
848                ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
849                ti = 0.0_wp
850             ENDIF
851!
852!--       Allocate array for uu
853          CASE ( 'uu' )
854             IF ( .NOT. ALLOCATED( uu ) )  THEN
855                ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) )
856                uu = 0.0_wp
857             ENDIF
858
859!
860!--       Allocate array for vv
861          CASE ( 'vv' )
862             IF ( .NOT. ALLOCATED( vv ) )  THEN
863                ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) )
864                vv = 0.0_wp
865             ENDIF
866
867!
868!--       Allocate array for ww
869          CASE ( 'ww' )
870             IF ( .NOT. ALLOCATED( ww ) )  THEN
871                ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
872                ww = 0.0_wp
873             ENDIF
874
875       END SELECT
876
877       ivar = ivar + 1
878    ENDDO
879
880    initialized_diagnostic_output_quantities = .TRUE.
881   
882 END SUBROUTINE doq_init
883
884
885!--------------------------------------------------------------------------------------------------!
886! Description:
887! ------------
888!> Calculate diagnostic quantities
889!--------------------------------------------------------------------------------------------------!
890 SUBROUTINE doq_calculate
891
892    IMPLICIT NONE
893
894    INTEGER(iwp) ::  i, j, k    !< grid loop index in x-, y-, z-direction
895    INTEGER(iwp) ::  ivar       !< loop index over all 2d/3d/mask output quantities
896
897
898!     CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
899
900!
901!-- Save timestep number to check in time_integration if doq_calculate
902!-- has been called already, since the CALL occurs at two locations, but the calculations need to be
903!-- done only once per timestep.
904    timestep_number_at_prev_calc = current_timestep_number
905
906    ivar = 1
907
908    DO  WHILE ( ivar <= SIZE( do_all ) ) 
909
910       SELECT CASE ( TRIM( do_all(ivar) ) )
911!
912!--       Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
913          CASE ( 'ti' )
914             DO  i = nxl, nxr
915                DO  j = nys, nyn
916                   DO  k = nzb+1, nzt
917                      ti(k,j,i) = 0.25_wp * SQRT(                              &
918                        (   (   w(k,j+1,i) + w(k-1,j+1,i)                      &
919                              - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy              &
920                          - (   v(k+1,j,i) + v(k+1,j+1,i)                      &
921                              - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2     &
922                      + (   (   u(k+1,j,i) + u(k+1,j,i+1)                      &
923                              - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k)          &
924                          - (   w(k,j,i+1) + w(k-1,j,i+1)                      &
925                              - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx     )**2     &
926                      + (   (   v(k,j,i+1) + v(k,j+1,i+1)                      &
927                              - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx              &
928                          - (   u(k,j+1,i) + u(k,j+1,i+1)                      &
929                              - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy     )**2  )  &
930                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) )
931                   ENDDO
932                ENDDO
933             ENDDO           
934!
935!--       uu
936          CASE ( 'uu' )
937             DO  i = nxl, nxr
938                DO  j = nys, nyn
939                   DO  k = nzb+1, nzt
940                      uu(k,j,i) = u(k,j,i) * u(k,j,i)                          &
941                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1) )
942                   ENDDO
943                ENDDO
944             ENDDO
945!
946!--       vv
947          CASE ( 'vv' )
948             DO  i = nxl, nxr
949                DO  j = nys, nyn
950                   DO  k = nzb+1, nzt
951                      vv(k,j,i) = v(k,j,i) * v(k,j,i)                          &
952                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2) )
953                   ENDDO
954                ENDDO
955             ENDDO
956!
957!--       ww
958          CASE ( 'ww' )
959             DO  i = nxl, nxr
960                DO  j = nys, nyn
961                   DO  k = nzb+1, nzt-1
962                      ww(k,j,i) = w(k,j,i) * w(k,j,i)                          &
963                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3) )
964                   ENDDO
965                ENDDO
966             ENDDO
967
968       END SELECT
969
970       ivar = ivar + 1
971    ENDDO
972
973!     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
974
975 END SUBROUTINE doq_calculate
976
977
978!------------------------------------------------------------------------------!
979! Description:
980! ------------
981!> ...
982!------------------------------------------------------------------------------!
983 SUBROUTINE doq_prepare
984
985
986    USE control_parameters,                                                    &
987        ONLY:  do2d, do3d, domask, masks
988
989    IMPLICIT NONE
990
991    CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) ::  do2d_var = ' '  !<
992                                                          !< label array for 2d output quantities
993
994    INTEGER(iwp) ::  av         !< index defining type of output, av=0 instantaneous, av=1 averaged
995    INTEGER(iwp) ::  ivar       !< loop index
996    INTEGER(iwp) ::  ivar_all   !< loop index
997    INTEGER(iwp) ::  l          !< index for cutting string
998    INTEGER(iwp) ::  mid          !< masked output running index
999
1000    prepared_diagnostic_output_quantities = .FALSE.
1001
1002    ivar     = 1
1003    ivar_all = 1
1004
1005    DO  av = 0, 1
1006!
1007!--    Remove _xy, _xz, or _yz from string
1008       l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1009       do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1010!
1011!--    Gather 2d output quantity names.
1012!--    Check for double occurrence of output quantity, e.g. by _xy,
1013!--    _yz, _xz.
1014       DO  WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
1015          IF ( .NOT.  ANY( do_all == do2d_var(av,ivar) ) )  THEN
1016             do_all(ivar_all) = do2d_var(av,ivar)
1017          ENDIF
1018          ivar = ivar + 1
1019          ivar_all = ivar_all + 1
1020          l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
1021          do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
1022       ENDDO
1023
1024       ivar = 1
1025!
1026!--    Gather 3d output quantity names
1027       DO  WHILE ( do3d(av,ivar)(1:1) /= ' ' )
1028          do_all(ivar_all) = do3d(av,ivar)
1029          ivar = ivar + 1
1030          ivar_all = ivar_all + 1
1031       ENDDO
1032
1033       ivar = 1
1034!
1035!--    Gather masked output quantity names. Also check for double output
1036!--    e.g. by different masks.
1037       DO  mid = 1, masks
1038          DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
1039             IF ( .NOT.  ANY( do_all == domask(mid,av,ivar) ) )  THEN
1040                do_all(ivar_all) = domask(mid,av,ivar)
1041             ENDIF
1042
1043             ivar = ivar + 1
1044             ivar_all = ivar_all + 1
1045          ENDDO
1046          ivar = 1
1047       ENDDO
1048
1049    ENDDO
1050
1051    prepared_diagnostic_output_quantities = .TRUE.
1052
1053 END SUBROUTINE doq_prepare
1054 
1055!------------------------------------------------------------------------------!
1056! Description:
1057! ------------
1058!> Subroutine reads local (subdomain) restart data
1059!> Note: With the current structure reading of non-standard array is not
1060!> possible
1061!------------------------------------------------------------------------------!
1062!  SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,             &
1063!                            nxr_on_file, nynf, nync, nyn_on_file, nysf,         &
1064!                            nysc, nys_on_file, tmp_3d_non_standard, found )
1065
1066!
1067!     USE control_parameters
1068!         
1069!     USE indices
1070!     
1071!     USE kinds
1072!     
1073!     USE pegrid
1074!
1075!
1076!     IMPLICIT NONE
1077!
1078!     INTEGER(iwp) ::  k               !<
1079!     INTEGER(iwp) ::  nxlc            !<
1080!     INTEGER(iwp) ::  nxlf            !<
1081!     INTEGER(iwp) ::  nxl_on_file     !<
1082!     INTEGER(iwp) ::  nxrc            !<
1083!     INTEGER(iwp) ::  nxrf            !<
1084!     INTEGER(iwp) ::  nxr_on_file     !<
1085!     INTEGER(iwp) ::  nync            !<
1086!     INTEGER(iwp) ::  nynf            !<
1087!     INTEGER(iwp) ::  nyn_on_file     !<
1088!     INTEGER(iwp) ::  nysc            !<
1089!     INTEGER(iwp) ::  nysf            !<
1090!     INTEGER(iwp) ::  nys_on_file     !<
1091!
1092!     LOGICAL, INTENT(OUT)  :: found
1093!     
1094!     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions
1095! !
1096! !-- If temporary non-standard array for reading is already allocated,
1097! !-- deallocate it.
1098!     IF ( ALLOCATED( tmp_3d_non_standard ) )  DEALLOCATE( tmp_3d_non_standard )
1099!     
1100!     found = .TRUE.
1101!
1102!     SELECT CASE ( restart_string(1:length) )
1103!
1104!        CASE ( 'ti_av' )
1105!           IF ( .NOT. ALLOCATED( ti_av ) )  THEN
1106!              ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1107!           ENDIF
1108!           IF ( k == 1 )  THEN
1109!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
1110!                                            nxl_on_file:nxr_on_file) )
1111!              READ ( 13 )  tmp_3d_non_standard
1112!           ENDIF
1113!           ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1114!     
1115!        CASE ( 'uu_av' )
1116!           IF ( .NOT. ALLOCATED( uu_av ) )  THEN
1117!              ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1118!           ENDIF
1119!           IF ( k == 1 )  THEN
1120!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
1121!                                            nxl_on_file:nxr_on_file) )
1122!              READ ( 13 )  tmp_3d_non_standard
1123!           ENDIF
1124!           uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1125!                   
1126!        CASE ( 'vv_av' )
1127!           IF ( .NOT. ALLOCATED( vv_av ) )  THEN
1128!              ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1129!           ENDIF
1130!           IF ( k == 1 )  THEN
1131!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
1132!                                            nxl_on_file:nxr_on_file) )
1133!              READ ( 13 )  tmp_3d_non_standard
1134!           ENDIF
1135!           vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1136!                   
1137!        CASE ( 'ww_av' )
1138!           IF ( .NOT. ALLOCATED( ww_av ) )  THEN
1139!              ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
1140!           ENDIF
1141!           IF ( k == 1 )  THEN
1142!              ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file,  &
1143!                                            nxl_on_file:nxr_on_file) )
1144!              READ ( 13 )  tmp_3d_non_standard
1145!           ENDIF
1146!           ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
1147!                         
1148!
1149!        CASE DEFAULT
1150!
1151!           found = .FALSE.
1152!
1153!     END SELECT
1154!
1155!  END SUBROUTINE doq_rrd_local
1156 
1157!------------------------------------------------------------------------------!
1158! Description:
1159! ------------
1160!> Subroutine writes local (subdomain) restart data
1161!------------------------------------------------------------------------------!
1162 SUBROUTINE doq_wrd_local
1163
1164
1165    IMPLICIT NONE
1166
1167    IF ( ALLOCATED( ti_av ) )  THEN
1168       CALL wrd_write_string( 'ti_av' ) 
1169       WRITE ( 14 )  ti_av
1170    ENDIF
1171   
1172    IF ( ALLOCATED( uu_av ) )  THEN
1173       CALL wrd_write_string( 'uu_av' ) 
1174       WRITE ( 14 )  uu_av
1175    ENDIF
1176       
1177    IF ( ALLOCATED( vv_av ) )  THEN
1178       CALL wrd_write_string( 'vv_av' ) 
1179       WRITE ( 14 )  vv_av
1180    ENDIF
1181       
1182    IF ( ALLOCATED( ww_av ) )  THEN
1183       CALL wrd_write_string( 'ww_av' ) 
1184       WRITE ( 14 )  ww_av
1185    ENDIF
1186
1187
1188 END SUBROUTINE doq_wrd_local
1189 
1190 
1191
1192 END MODULE diagnostic_output_quantities_mod
Note: See TracBrowser for help on using the repository browser.