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

Last change on this file since 4147 was 4132, checked in by suehring, 5 years ago

Bugfix in masked data output for prognostic quantities

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