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

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

diagnostic output: Modularize diagnostic output, rename subroutines; formatting adjustments; allocate arrays only when required; add output of uu, vv, ww to enable variance calculation via temporal EC method; radiation: bugfix in masked data output; flow_statistics: Correct conversion to kinematic vertical scalar fluxes in case of pw-scheme and statistic regions

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