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

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

chem_emissions: Replace global arrays also in mode_emis branch; diagnostic output: restructure initialization in order to work also when data output during spin-up is enabled; radiation: give informative message on raytracing distance only by core zero not by all cores

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