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

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