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

Last change on this file since 4069 was 4069, checked in by Giersch, 2 years ago

Bugfix for masked output, compiler warning removed, test case for wind turbine model revised

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