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

!> @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 <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2019 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: diagnostic_output_quantities_mod.f90 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:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
!                precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr, s, tend,  &
!                u, v, vpt, w, zu, zw, waterflux_output_conversion, hyrho, d_exner
! 
!     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 indices,                                                               &
!         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,       &
!                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.
    LOGICAL ::  prepared_diagnostic_output_quantities = .FALSE.

    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


    SAVE

    PRIVATE

!
!-- Public variables
    PUBLIC do_all,                                                                                 &
           initialized_diagnostic_output_quantities,                                               &
           prepared_diagnostic_output_quantities,                                                  &
           ti, ti_av,                                                                              &
           timestep_number_at_prev_calc
!
!-- Public routines
    PUBLIC diagnostic_output_quantities_calculate,                                                 &
           diagnostic_output_quantities_init,                                                      &
           diagnostic_output_quantities_prepare


    INTERFACE diagnostic_output_quantities_init
       MODULE PROCEDURE diagnostic_output_quantities_init
    END INTERFACE diagnostic_output_quantities_init

    INTERFACE diagnostic_output_quantities_calculate
       MODULE PROCEDURE diagnostic_output_quantities_calculate
    END INTERFACE diagnostic_output_quantities_calculate

    INTERFACE diagnostic_output_quantities_prepare
       MODULE PROCEDURE diagnostic_output_quantities_prepare
    END INTERFACE diagnostic_output_quantities_prepare


 CONTAINS

!------------------------------------------------------------------------------!
! Description:
! ------------
!> ...
!------------------------------------------------------------------------------!
 SUBROUTINE diagnostic_output_quantities_init


    USE indices,                                                                                   &
        ONLY:  nxl, nxr, nyn, nys, nzb, nzt

    IMPLICIT NONE
!
!-- Next line is to avoid compiler warnings about unused variables
    IF ( timestep_number_at_prev_calc == 0 )  CONTINUE

    initialized_diagnostic_output_quantities = .FALSE.

    IF ( .NOT. ALLOCATED( ti ) )  THEN
       ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
       ti = 0.0_wp
    ENDIF

    initialized_diagnostic_output_quantities = .TRUE.

 END SUBROUTINE diagnostic_output_quantities_init


!--------------------------------------------------------------------------------------------------!
! Description:
! ------------
!> ...
!--------------------------------------------------------------------------------------------------!
 SUBROUTINE diagnostic_output_quantities_calculate


    USE arrays_3d,                                                                                 &
        ONLY:  ddzu, u, v, w

    USE grid_variables,                                                                            &
        ONLY:  ddx, ddy

    USE indices,                                                                                   &
        ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0

    IMPLICIT NONE

    INTEGER(iwp) ::  i, j, k    !< grid loop index in x-, y-, z-direction
    INTEGER(iwp) ::  ivar_all   !< 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 diagnostic_output_quantities_prepare
    IF ( .NOT.  initialized_diagnostic_output_quantities )  CALL diagnostic_output_quantities_init
!
!-- Save timestep number to check in time_integration if diagnostic_output_quantities_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_all = 1

    DO  WHILE ( do_all(ivar_all)(1:1) /= ' ' )

       SELECT CASE ( TRIM( do_all(ivar_all) ) )
!
!--       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

          END SELECT

          ivar_all = ivar_all + 1
    ENDDO

!     CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )

 END SUBROUTINE diagnostic_output_quantities_calculate


!------------------------------------------------------------------------------!
! Description:
! ------------
!> ...
!------------------------------------------------------------------------------!
 SUBROUTINE diagnostic_output_quantities_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)(0:l-2) = do2d(av,ivar)(0:l-2)
!
!--    Gather 2d output quantity names, check for double occurrence of output quantity
       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)(0:l-2) = do2d(av,ivar)(0:l-2)
       ENDDO

       ivar = 1
!
!--    Gather 3d output quantity names, check for double occurrence of output quantity
       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, check for double occurrence of output quantity
       DO  mid = 1, masks
          DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
             do_all(ivar_all) = domask(mid,av,ivar)

             ivar = ivar + 1
             ivar_all = ivar_all + 1
          ENDDO
          ivar = 1
       ENDDO

    ENDDO

    prepared_diagnostic_output_quantities = .TRUE.

 END SUBROUTINE diagnostic_output_quantities_prepare

 END MODULE diagnostic_output_quantities_mod