!> @file diagnostic_output_quantities_mod.f90 !------------------------------------------------------------------------------! ! This file is part of the PALM model system. ! ! PALM is free software: you can redistribute it and/or modify it under the ! terms of the GNU General Public License as published by the Free Software ! Foundation, either version 3 of the License, or (at your option) any later ! version. ! ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License along with ! PALM. If not, see . ! ! Copyright 1997-2019 Leibniz Universitaet Hannover !------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: diagnostic_output_quantities_mod.f90 4039 2019-06-18 10:32:41Z raasch $ ! - Add output of uu, vv, ww to enable variance calculation according temporal ! EC method ! - Allocate arrays only when they are required ! - Formatting adjustment ! - Rename subroutines ! - Further modularization ! ! 3998 2019-05-23 13:38:11Z suehring ! Bugfix in gathering all output strings ! ! 3995 2019-05-22 18:59:54Z suehring ! Avoid compiler warnings about unused variable and fix string operation which ! is not allowed with PGI compiler ! ! 3994 2019-05-22 18:08:09Z suehring ! Initial revision ! ! ! @author Farah Kanani-Suehring ! ! Description: ! ------------ !> ... !------------------------------------------------------------------------------! MODULE diagnostic_output_quantities_mod USE arrays_3d, & ONLY: ddzu, u, v, w ! ! USE averaging ! ! USE basic_constants_and_equations_mod, & ! ONLY: c_p, lv_d_cp, l_v ! ! USE bulk_cloud_model_mod, & ! ONLY: bulk_cloud_model ! USE control_parameters, & ONLY: current_timestep_number, varnamelength ! ! USE cpulog, & ! ONLY: cpu_log, log_point USE grid_variables, & ONLY: ddx, ddy ! USE indices, & ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0 ! USE kinds ! ! USE land_surface_model_mod, & ! ONLY: zs ! ! USE module_interface, & ! ONLY: module_interface_data_output_2d ! ! #if defined( __netcdf ) ! USE NETCDF ! #endif ! ! USE netcdf_interface, & ! ONLY: fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d, & ! id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat, & ! netcdf_data_format, netcdf_handle_error ! ! USE particle_attributes, & ! ONLY: grid_particles, number_of_particles, particle_advection_start, & ! particles, prt_count ! ! USE pegrid ! ! USE surface_mod, & ! ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h, & ! surf_lsm_h, surf_usm_h ! ! USE turbulence_closure_mod, & ! ONLY: tcm_data_output_2d IMPLICIT NONE CHARACTER(LEN=varnamelength), DIMENSION(500) :: do_all = ' ' INTEGER(iwp) :: timestep_number_at_prev_calc = 0 !< ...at previous diagnostic output calculation LOGICAL :: initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is p REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu !< uu REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv !< vv REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww !< ww REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww_av !< mean of ww SAVE PRIVATE ! !-- Public variables PUBLIC do_all, & initialized_diagnostic_output_quantities, & prepared_diagnostic_output_quantities, & timestep_number_at_prev_calc, & ti_av, & uu_av, & vv_av, & ww_av ! !-- Public routines PUBLIC doq_3d_data_averaging, & doq_calculate, & doq_check_data_output, & doq_define_netcdf_grid, & doq_output_2d, & doq_output_3d, & doq_output_mask, & doq_init, & doq_prepare, & doq_wrd_local ! doq_rrd_local, & INTERFACE doq_3d_data_averaging MODULE PROCEDURE doq_3d_data_averaging END INTERFACE doq_3d_data_averaging INTERFACE doq_calculate MODULE PROCEDURE doq_calculate END INTERFACE doq_calculate INTERFACE doq_check_data_output MODULE PROCEDURE doq_check_data_output END INTERFACE doq_check_data_output INTERFACE doq_define_netcdf_grid MODULE PROCEDURE doq_define_netcdf_grid END INTERFACE doq_define_netcdf_grid INTERFACE doq_output_2d MODULE PROCEDURE doq_output_2d END INTERFACE doq_output_2d INTERFACE doq_output_3d MODULE PROCEDURE doq_output_3d END INTERFACE doq_output_3d INTERFACE doq_output_mask MODULE PROCEDURE doq_output_mask END INTERFACE doq_output_mask INTERFACE doq_init MODULE PROCEDURE doq_init END INTERFACE doq_init INTERFACE doq_prepare MODULE PROCEDURE doq_prepare END INTERFACE doq_prepare ! INTERFACE doq_rrd_local ! MODULE PROCEDURE doq_rrd_local ! END INTERFACE doq_rrd_local INTERFACE doq_wrd_local MODULE PROCEDURE doq_wrd_local END INTERFACE doq_wrd_local CONTAINS !------------------------------------------------------------------------------! ! Description: ! ------------ !> Sum up and time-average diagnostic output quantities as well as allocate !> the array necessary for storing the average. !------------------------------------------------------------------------------! SUBROUTINE doq_3d_data_averaging( mode, variable ) USE control_parameters, & ONLY: average_count_3d CHARACTER (LEN=*) :: mode !< CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: i !< INTEGER(iwp) :: j !< INTEGER(iwp) :: k !< IF ( mode == 'allocate' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'ti' ) IF ( .NOT. ALLOCATED( ti_av ) ) THEN ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ENDIF ti_av = 0.0_wp CASE ( 'uu' ) IF ( .NOT. ALLOCATED( uu_av ) ) THEN ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ENDIF uu_av = 0.0_wp CASE ( 'vv' ) IF ( .NOT. ALLOCATED( vv_av ) ) THEN ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ENDIF vv_av = 0.0_wp CASE ( 'ww' ) IF ( .NOT. ALLOCATED( ww_av ) ) THEN ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ENDIF ww_av = 0.0_wp CASE DEFAULT CONTINUE END SELECT ELSEIF ( mode == 'sum' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'ti' ) IF ( ALLOCATED( ti_av ) ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'uu' ) IF ( ALLOCATED( uu_av ) ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'vv' ) IF ( ALLOCATED( vv_av ) ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE ( 'ww' ) IF ( ALLOCATED( ww_av ) ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i) ENDDO ENDDO ENDDO ENDIF CASE DEFAULT CONTINUE END SELECT ELSEIF ( mode == 'average' ) THEN SELECT CASE ( TRIM( variable ) ) CASE ( 'ti' ) IF ( ALLOCATED( ti_av ) ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE ( 'uu' ) IF ( ALLOCATED( uu_av ) ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE ( 'vv' ) IF ( ALLOCATED( vv_av ) ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF CASE ( 'ww' ) IF ( ALLOCATED( ww_av ) ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb, nzt+1 ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp ) ENDDO ENDDO ENDDO ENDIF END SELECT ENDIF END SUBROUTINE doq_3d_data_averaging !------------------------------------------------------------------------------! ! Description: ! ------------ !> Check data output for diagnostic output !------------------------------------------------------------------------------! SUBROUTINE doq_check_data_output( var, unit ) IMPLICIT NONE CHARACTER (LEN=*) :: unit !< CHARACTER (LEN=*) :: var !< SELECT CASE ( TRIM( var ) ) CASE ( 'ti' ) unit = '1/s' CASE ( 'uu' ) unit = 'm2/s2' CASE ( 'vv' ) unit = 'm2/s2' CASE ( 'ww' ) unit = 'm2/s2' CASE DEFAULT unit = 'illegal' END SELECT END SUBROUTINE doq_check_data_output !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining appropriate grid for netcdf variables. !------------------------------------------------------------------------------! SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z ) IMPLICIT NONE CHARACTER (LEN=*), INTENT(IN) :: variable !< LOGICAL, INTENT(OUT) :: found !< CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< found = .TRUE. SELECT CASE ( TRIM( variable ) ) ! !-- s grid CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz' ) grid_x = 'x' grid_y = 'y' grid_z = 'zu' ! !-- u grid CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' ) grid_x = 'xu' grid_y = 'y' grid_z = 'zu' ! !-- v grid CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz' ) grid_x = 'x' grid_y = 'yv' grid_z = 'zu' ! !-- w grid CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz' ) grid_x = 'x' grid_y = 'y' grid_z = 'zw' CASE DEFAULT found = .FALSE. grid_x = 'none' grid_y = 'none' grid_z = 'none' END SELECT END SUBROUTINE doq_define_netcdf_grid !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining 2D output variables !------------------------------------------------------------------------------! SUBROUTINE doq_output_2d( av, variable, found, grid, & mode, local_pf, two_d, nzb_do, nzt_do, fill_value ) IMPLICIT NONE CHARACTER (LEN=*) :: grid !< CHARACTER (LEN=*) :: mode !< CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: av !< value indicating averaged or non-averaged output INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w) INTEGER(iwp) :: i !< grid index x-direction INTEGER(iwp) :: j !< grid index y-direction INTEGER(iwp) :: k !< grid index z-direction INTEGER(iwp) :: nzb_do !< INTEGER(iwp) :: nzt_do !< LOGICAL :: found !< true if variable is in list LOGICAL :: resorted !< true if array is resorted LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections) REAL(wp) :: fill_value !< value for the _FillValue attribute REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output flag_nr = 0 found = .TRUE. resorted = .FALSE. two_d = .FALSE. SELECT CASE ( TRIM( variable ) ) CASE ( 'ti_xy', 'ti_xz', 'ti_yz' ) IF ( av == 0 ) THEN to_be_resorted => ti ELSE IF ( .NOT. ALLOCATED( ti_av ) ) THEN ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ti_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => ti_av ENDIF flag_nr = 0 IF ( mode == 'xy' ) grid = 'zu' CASE ( 'uu_xy', 'uu_xz', 'uu_yz' ) IF ( av == 0 ) THEN to_be_resorted => uu ELSE IF ( .NOT. ALLOCATED( uu_av ) ) THEN ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) uu_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => uu_av ENDIF flag_nr = 1 IF ( mode == 'xy' ) grid = 'zu' CASE ( 'vv_xy', 'vv_xz', 'vv_yz' ) IF ( av == 0 ) THEN to_be_resorted => vv ELSE IF ( .NOT. ALLOCATED( vv_av ) ) THEN ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) vv_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => vv_av ENDIF flag_nr = 2 IF ( mode == 'xy' ) grid = 'zu' CASE ( 'ww_xy', 'ww_xz', 'ww_yz' ) IF ( av == 0 ) THEN to_be_resorted => ww ELSE IF ( .NOT. ALLOCATED( ww_av ) ) THEN ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ww_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => ww_av ENDIF flag_nr = 3 IF ( mode == 'xy' ) grid = 'zw' CASE DEFAULT found = .FALSE. grid = 'none' END SELECT IF ( found .AND. .NOT. resorted ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & REAL( fill_value, KIND = wp ), & BTEST( wall_flags_0(k,j,i), flag_nr ) ) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE doq_output_2d !------------------------------------------------------------------------------! ! ! Description: ! ------------ !> Subroutine defining 3D output variables !------------------------------------------------------------------------------! SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, & nzt_do ) IMPLICIT NONE CHARACTER (LEN=*) :: variable !< INTEGER(iwp) :: av !< index indicating averaged or instantaneous output INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w) INTEGER(iwp) :: i !< index variable along x-direction INTEGER(iwp) :: j !< index variable along y-direction INTEGER(iwp) :: k !< index variable along z-direction INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) LOGICAL :: found !< true if variable is in list LOGICAL :: resorted !< true if array is resorted REAL(wp) :: fill_value !< value for the _FillValue attribute REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output flag_nr = 0 found = .TRUE. resorted = .FALSE. SELECT CASE ( TRIM( variable ) ) CASE ( 'ti' ) IF ( av == 0 ) THEN to_be_resorted => ti ELSE IF ( .NOT. ALLOCATED( ti_av ) ) THEN ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ti_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => ti_av ENDIF flag_nr = 0 CASE ( 'uu' ) IF ( av == 0 ) THEN to_be_resorted => uu ELSE IF ( .NOT. ALLOCATED( uu_av ) ) THEN ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) uu_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => uu_av ENDIF flag_nr = 1 CASE ( 'vv' ) IF ( av == 0 ) THEN to_be_resorted => vv ELSE IF ( .NOT. ALLOCATED( vv_av ) ) THEN ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) vv_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => vv_av ENDIF flag_nr = 2 CASE ( 'ww' ) IF ( av == 0 ) THEN to_be_resorted => ww ELSE IF ( .NOT. ALLOCATED( ww_av ) ) THEN ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ww_av = REAL( fill_value, KIND = wp ) ENDIF to_be_resorted => ww_av ENDIF flag_nr = 3 CASE DEFAULT found = .FALSE. END SELECT IF ( found .AND. .NOT. resorted ) THEN DO i = nxl, nxr DO j = nys, nyn DO k = nzb_do, nzt_do local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & REAL( fill_value, KIND = wp ), & BTEST( wall_flags_0(k,j,i), flag_nr ) ) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE doq_output_3d ! Description: ! ------------ !> Resorts the user-defined output quantity with indices (k,j,i) to a !> temporary array with indices (i,j,k) for masked data output. !------------------------------------------------------------------------------! SUBROUTINE doq_output_mask( av, variable, found, local_pf ) USE control_parameters USE indices USE surface_mod, & ONLY: get_topography_top_index_ji IMPLICIT NONE CHARACTER (LEN=*) :: variable !< CHARACTER (LEN=5) :: grid !< flag to distinquish between staggered grids INTEGER(iwp) :: av !< index indicating averaged or instantaneous output INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w) INTEGER(iwp) :: i !< index variable along x-direction INTEGER(iwp) :: j !< index variable along y-direction INTEGER(iwp) :: k !< index variable along z-direction INTEGER(iwp) :: topo_top_ind !< k index of highest horizontal surface LOGICAL :: found !< true if variable is in list LOGICAL :: resorted !< true if array is resorted REAL(wp), & DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: & local_pf !< REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output flag_nr = 0 found = .TRUE. resorted = .FALSE. grid = 's' SELECT CASE ( TRIM( variable ) ) CASE ( 'ti' ) IF ( av == 0 ) THEN to_be_resorted => ti ELSE to_be_resorted => ti_av ENDIF grid = 's' flag_nr = 0 CASE ( 'uu' ) IF ( av == 0 ) THEN to_be_resorted => uu ELSE to_be_resorted => uu_av ENDIF grid = 'u' flag_nr = 1 CASE ( 'vv' ) IF ( av == 0 ) THEN to_be_resorted => vv ELSE to_be_resorted => vv_av ENDIF grid = 'v' flag_nr = 2 CASE ( 'ww' ) IF ( av == 0 ) THEN to_be_resorted => ww ELSE to_be_resorted => ww_av ENDIF grid = 'w' flag_nr = 3 CASE DEFAULT found = .FALSE. END SELECT IF ( found .AND. .NOT. resorted ) THEN IF ( .NOT. mask_surface(mid) ) THEN ! !-- Default masked output DO i = 1, mask_size_l(mid,1) DO j = 1, mask_size_l(mid,2) DO k = 1, mask_size_l(mid,3) local_pf(i,j,k) = to_be_resorted(mask_k(mid,k), & mask_j(mid,j), & mask_i(mid,i)) ENDDO ENDDO ENDDO ELSE ! !-- Terrain-following masked output DO i = 1, mask_size_l(mid,1) DO j = 1, mask_size_l(mid,2) ! !-- Get k index of highest horizontal surface topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), & mask_i(mid,i), & grid ) ! !-- Save output array DO k = 1, mask_size_l(mid,3) local_pf(i,j,k) = to_be_resorted( & MIN( topo_top_ind+mask_k(mid,k), & nzt+1 ), & mask_j(mid,j), & mask_i(mid,i) ) ENDDO ENDDO ENDDO ENDIF ENDIF END SUBROUTINE doq_output_mask !------------------------------------------------------------------------------! ! Description: ! ------------ !> Allocate required arrays !------------------------------------------------------------------------------! SUBROUTINE doq_init IMPLICIT NONE INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities ! !-- Next line is to avoid compiler warnings about unused variables IF ( timestep_number_at_prev_calc == 0 ) CONTINUE initialized_diagnostic_output_quantities = .FALSE. ivar = 1 DO WHILE ( ivar <= SIZE( do_all ) ) SELECT CASE ( TRIM( do_all(ivar) ) ) ! !-- Allocate array for 'turbulence intensity' CASE ( 'ti' ) IF ( .NOT. ALLOCATED( ti ) ) THEN ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) ) ti = 0.0_wp ENDIF ! !-- Allocate array for uu CASE ( 'uu' ) IF ( .NOT. ALLOCATED( uu ) ) THEN ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) ) uu = 0.0_wp ENDIF ! !-- Allocate array for vv CASE ( 'vv' ) IF ( .NOT. ALLOCATED( vv ) ) THEN ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) ) vv = 0.0_wp ENDIF ! !-- Allocate array for ww CASE ( 'ww' ) IF ( .NOT. ALLOCATED( ww ) ) THEN ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) ) ww = 0.0_wp ENDIF END SELECT ivar = ivar + 1 ENDDO initialized_diagnostic_output_quantities = .TRUE. END SUBROUTINE doq_init !--------------------------------------------------------------------------------------------------! ! Description: ! ------------ !> Calculate diagnostic quantities !--------------------------------------------------------------------------------------------------! SUBROUTINE doq_calculate IMPLICIT NONE INTEGER(iwp) :: i, j, k !< grid loop index in x-, y-, z-direction INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities ! CALL cpu_log( log_point(41), 'calculate_quantities', 'start' ) ! !-- Preparatory steps and initialization of output arrays IF ( .NOT. prepared_diagnostic_output_quantities ) CALL doq_prepare IF ( .NOT. initialized_diagnostic_output_quantities ) CALL doq_init ! !-- Save timestep number to check in time_integration if doq_calculate !-- has been called already, since the CALL occurs at two locations, but the calculations need to be !-- done only once per timestep. timestep_number_at_prev_calc = current_timestep_number ivar = 1 DO WHILE ( ivar <= SIZE( do_all ) ) SELECT CASE ( TRIM( do_all(ivar) ) ) ! !-- Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point CASE ( 'ti' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt ti(k,j,i) = 0.25_wp * SQRT( & ( ( w(k,j+1,i) + w(k-1,j+1,i) & - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy & - ( v(k+1,j,i) + v(k+1,j+1,i) & - v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2 & + ( ( u(k+1,j,i) + u(k+1,j,i+1) & - u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k) & - ( w(k,j,i+1) + w(k-1,j,i+1) & - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx )**2 & + ( ( v(k,j,i+1) + v(k,j+1,i+1) & - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx & - ( u(k,j+1,i) + u(k,j+1,i+1) & - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy )**2 ) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) ) ENDDO ENDDO ENDDO ! !-- uu CASE ( 'uu' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt uu(k,j,i) = u(k,j,i) * u(k,j,i) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1) ) ENDDO ENDDO ENDDO ! !-- vv CASE ( 'vv' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt vv(k,j,i) = v(k,j,i) * v(k,j,i) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2) ) ENDDO ENDDO ENDDO ! !-- ww CASE ( 'ww' ) DO i = nxl, nxr DO j = nys, nyn DO k = nzb+1, nzt-1 ww(k,j,i) = w(k,j,i) * w(k,j,i) & * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3) ) ENDDO ENDDO ENDDO END SELECT ivar = ivar + 1 ENDDO ! CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' ) END SUBROUTINE doq_calculate !------------------------------------------------------------------------------! ! Description: ! ------------ !> ... !------------------------------------------------------------------------------! SUBROUTINE doq_prepare USE control_parameters, & ONLY: do2d, do3d, domask, masks, mid IMPLICIT NONE CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) :: do2d_var = ' ' !< !< label array for 2d output quantities INTEGER(iwp) :: av !< index defining type of output, av=0 instantaneous, av=1 averaged INTEGER(iwp) :: ivar !< loop index INTEGER(iwp) :: ivar_all !< loop index INTEGER(iwp) :: l !< index for cutting string prepared_diagnostic_output_quantities = .FALSE. ivar = 1 ivar_all = 1 DO av = 0, 1 ! !-- Remove _xy, _xz, or _yz from string l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) ) do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3) ! !-- Gather 2d output quantity names. !-- Check for double occurrence of output quantity, e.g. by _xy, !-- _yz, _xz. DO WHILE ( do2d_var(av,ivar)(1:1) /= ' ' ) IF ( .NOT. ANY( do_all == do2d_var(av,ivar) ) ) THEN do_all(ivar_all) = do2d_var(av,ivar) ENDIF ivar = ivar + 1 ivar_all = ivar_all + 1 l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) ) do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3) ENDDO ivar = 1 ! !-- Gather 3d output quantity names DO WHILE ( do3d(av,ivar)(1:1) /= ' ' ) do_all(ivar_all) = do3d(av,ivar) ivar = ivar + 1 ivar_all = ivar_all + 1 ENDDO ivar = 1 ! !-- Gather masked output quantity names. Also check for double output !-- e.g. by different masks. DO mid = 1, masks DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' ) IF ( .NOT. ANY( do_all == domask(mid,av,ivar) ) ) THEN do_all(ivar_all) = do2d_var(av,ivar) ENDIF ivar = ivar + 1 ivar_all = ivar_all + 1 ENDDO ivar = 1 ENDDO ENDDO prepared_diagnostic_output_quantities = .TRUE. END SUBROUTINE doq_prepare !------------------------------------------------------------------------------! ! Description: ! ------------ !> Subroutine reads local (subdomain) restart data !> Note: With the current structure reading of non-standard array is not !> possible !------------------------------------------------------------------------------! ! SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, & ! nxr_on_file, nynf, nync, nyn_on_file, nysf, & ! nysc, nys_on_file, tmp_3d_non_standard, found ) ! ! ! USE control_parameters ! ! USE indices ! ! USE kinds ! ! USE pegrid ! ! ! IMPLICIT NONE ! ! INTEGER(iwp) :: k !< ! INTEGER(iwp) :: nxlc !< ! INTEGER(iwp) :: nxlf !< ! INTEGER(iwp) :: nxl_on_file !< ! INTEGER(iwp) :: nxrc !< ! INTEGER(iwp) :: nxrf !< ! INTEGER(iwp) :: nxr_on_file !< ! INTEGER(iwp) :: nync !< ! INTEGER(iwp) :: nynf !< ! INTEGER(iwp) :: nyn_on_file !< ! INTEGER(iwp) :: nysc !< ! INTEGER(iwp) :: nysf !< ! INTEGER(iwp) :: nys_on_file !< ! ! LOGICAL, INTENT(OUT) :: found ! ! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions ! ! ! !-- If temporary non-standard array for reading is already allocated, ! !-- deallocate it. ! IF ( ALLOCATED( tmp_3d_non_standard ) ) DEALLOCATE( tmp_3d_non_standard ) ! ! found = .TRUE. ! ! SELECT CASE ( restart_string(1:length) ) ! ! CASE ( 'ti_av' ) ! IF ( .NOT. ALLOCATED( ti_av ) ) THEN ! ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ! ENDIF ! IF ( k == 1 ) THEN ! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, & ! nxl_on_file:nxr_on_file) ) ! READ ( 13 ) tmp_3d_non_standard ! ENDIF ! ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) ! ! CASE ( 'uu_av' ) ! IF ( .NOT. ALLOCATED( uu_av ) ) THEN ! ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ! ENDIF ! IF ( k == 1 ) THEN ! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, & ! nxl_on_file:nxr_on_file) ) ! READ ( 13 ) tmp_3d_non_standard ! ENDIF ! uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) ! ! CASE ( 'vv_av' ) ! IF ( .NOT. ALLOCATED( vv_av ) ) THEN ! ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ! ENDIF ! IF ( k == 1 ) THEN ! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, & ! nxl_on_file:nxr_on_file) ) ! READ ( 13 ) tmp_3d_non_standard ! ENDIF ! vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) ! ! CASE ( 'ww_av' ) ! IF ( .NOT. ALLOCATED( ww_av ) ) THEN ! ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) ! ENDIF ! IF ( k == 1 ) THEN ! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, & ! nxl_on_file:nxr_on_file) ) ! READ ( 13 ) tmp_3d_non_standard ! ENDIF ! ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf) ! ! ! CASE DEFAULT ! ! found = .FALSE. ! ! END SELECT ! ! END SUBROUTINE doq_rrd_local !------------------------------------------------------------------------------! ! Description: ! ------------ !> Subroutine writes local (subdomain) restart data !------------------------------------------------------------------------------! SUBROUTINE doq_wrd_local IMPLICIT NONE IF ( ALLOCATED( ti_av ) ) THEN CALL wrd_write_string( 'ti_av' ) WRITE ( 14 ) ti_av ENDIF IF ( ALLOCATED( uu_av ) ) THEN CALL wrd_write_string( 'uu_av' ) WRITE ( 14 ) uu_av ENDIF IF ( ALLOCATED( vv_av ) ) THEN CALL wrd_write_string( 'vv_av' ) WRITE ( 14 ) vv_av ENDIF IF ( ALLOCATED( ww_av ) ) THEN CALL wrd_write_string( 'ww_av' ) WRITE ( 14 ) ww_av ENDIF END SUBROUTINE doq_wrd_local END MODULE diagnostic_output_quantities_mod