source: palm/trunk/SOURCE/biometeorology_mod.f90 @ 3496

!> @file biometeorology_mod.f90
!--------------------------------------------------------------------------------!
! This file is part of PALM-4U.
!
! PALM-4U 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-4U 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 2018-2018 Deutscher Wetterdienst (DWD)
! Copyright 2018-2018 Institute of Computer Science, Academy of Sciences, Prague
! Copyright 2018-2018 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
! 
! 
! Former revisions:
! -----------------
! $Id: biometeorology_mod.f90 3479 2018-11-01 16:00:30Z suehring $
! - reworked check for output quantities
! - assign new palm-error numbers
! - set unit according to data standard.
! 
! 3475 2018-10-30 21:16:31Z kanani
! Add option for using MRT from RTM instead of simplified global value
! 
! 3464 2018-10-30 18:08:55Z kanani
! From branch resler@3462, pavelkrc:
! make use of basic_constants_and_equations_mod
! 
! 3448 2018-10-29 18:14:31Z kanani
! Initial revision
!
! 
!
! Authors:
! --------
! @author Dominik Froehlich <dominik.froehlich@dwd.de>
! @author Jaroslav Resler <resler@cs.cas.cz>
!
!
! Description:
! ------------
!> Human thermal comfort module calculating thermal perception of a sample
!> human being under the current meteorological conditions.
!>
!> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
!>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
!> @todo Comments start with capital letter --> "!-- Include..." 
!> @todo Variable and routine names strictly in lowercase letters and in English
!>
!> @note nothing now
!>
!> @bug  no known bugs by now
!------------------------------------------------------------------------------!
 MODULE biometeorology_mod

    USE arrays_3d,                                                             &
       ONLY:  pt, p, u, v, w, q

    USE averaging,                                                             &
       ONLY:  pt_av, q_av, u_av, v_av, w_av

    USE basic_constants_and_equations_mod,                                     &
       ONLY:  degc_to_k, magnus, sigma_sb

    USE biometeorology_ipt_mod

    USE biometeorology_pet_mod

    USE biometeorology_pt_mod,                                                 &
       ONLY: calculate_pt_static

    USE biometeorology_utci_mod

    USE control_parameters,                                                    &
       ONLY:  average_count_3d, biometeorology, dz, dz_stretch_factor,         &
              dz_stretch_level, humidity, initializing_actions, nz_do3d,       &
              simulated_time, surface_pressure

    USE grid_variables,                                                        &
       ONLY:  ddx, dx, ddy, dy

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

    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM

!-- Import radiation model to obtain input for mean radiant temperature
    USE radiation_model_mod,                                                   &
       ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                      &
             mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation,                       &
             radiation_interactions, rad_sw_in,                                &
             rad_sw_out, rad_lw_in, rad_lw_out

    USE surface_mod,                                                           &
       ONLY: get_topography_top_index_ji

    IMPLICIT NONE

    PRIVATE

!-- Declare all global variables within the module (alphabetical order)
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (°C)
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_grid    !< PT results   (°C)
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_grid  !< UTCI results (°C)
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_grid   !< PET results  (°C)
!-- Grids for averaged thermal indices       
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_av_grid    !< PT results (aver. input)   (°C)
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av_grid  !< UTCI results (aver. input) (°C)
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av_grid   !< PET results (aver. input)  (°C)
!-- Grids for radiation_model
    REAL(wp), DIMENSION(:), ALLOCATABLE ::  biom_mrt     !< biometeorology mean radiant temperature for each MRT box
    REAL(wp), DIMENSION(:), ALLOCATABLE ::  biom_mrt_av  !< time average mean

    INTEGER( iwp ) ::  biom_cell_level     !< cell level biom calculates for
    REAL ( wp )    ::  biom_output_height  !< height output is calculated in m
    REAL ( wp )    ::  time_biom_results   !< the time, the last set of biom results have been calculated for
    REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
    REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
!--

    LOGICAL ::  aver_pt  = .FALSE.  !< switch: do pt averaging in this module? (if .FALSE. this is done globally)
    LOGICAL ::  aver_q   = .FALSE.  !< switch: do e  averaging in this module?
    LOGICAL ::  aver_u   = .FALSE.  !< switch: do u  averaging in this module?
    LOGICAL ::  aver_v   = .FALSE.  !< switch: do v  averaging in this module?
    LOGICAL ::  aver_w   = .FALSE.  !< switch: do w  averaging in this module?
    LOGICAL ::  average_trigger_pt   = .FALSE.  !< update averaged input on call to biom_pt?
    LOGICAL ::  average_trigger_utci = .FALSE.  !< update averaged input on call to biom_utci?
    LOGICAL ::  average_trigger_pet  = .FALSE.  !< update averaged input on call to biom_pet?

    LOGICAL ::  mrt_from_rtm   = .TRUE.   !< switch: use mrt calculated by RTM for calculation of thermal indices
    LOGICAL ::  biom_pt        = .TRUE.   !< Turn index PT (instant. input) on or off
    LOGICAL ::  biom_pt_av     = .TRUE.   !< Turn index PT (averaged input) on or off
    LOGICAL ::  biom_pet       = .TRUE.   !< Turn index PET (instant. input) on or off
    LOGICAL ::  biom_pet_av    = .TRUE.   !< Turn index PET (averaged input) on or off
    LOGICAL ::  biom_utci      = .TRUE.   !< Turn index UTCI (instant. input) on or off
    LOGICAL ::  biom_utci_av   = .TRUE.   !< Turn index UTCI (averaged input) on or off

!-- Add INTERFACES that must be available to other modules (alphabetical order)

    PUBLIC biom_3d_data_averaging, biom_check_data_output,                     &
    biom_calculate_static_grid, biom_calc_ipt,                                 &
    biom_check_parameters, biom_data_output_3d, biom_data_output_2d,           &
    biom_define_netcdf_grid, biom_determine_input_at, biom_header, biom_init,  &
    biom_init_arrays, biom_parin, biom_pt, biom_pt_av, biom_pet, biom_pet_av,  &
    biom_utci, biom_utci_av, time_biom_results

!
!-- PALM interfaces:
!
!-- 3D averaging for HTCM _INPUT_ variables
    INTERFACE biom_3d_data_averaging
       MODULE PROCEDURE biom_3d_data_averaging
    END INTERFACE biom_3d_data_averaging

!-- Calculate static thermal indices PT, UTCI and/or PET
    INTERFACE biom_calculate_static_grid
       MODULE PROCEDURE biom_calculate_static_grid
    END INTERFACE biom_calculate_static_grid

!-- Calculate the dynamic index iPT (to be caled by the agent model)
    INTERFACE biom_calc_ipt
       MODULE PROCEDURE biom_calc_ipt
    END INTERFACE biom_calc_ipt

!-- Data output checks for 2D/3D data to be done in check_parameters
     INTERFACE biom_check_data_output
        MODULE PROCEDURE biom_check_data_output
     END INTERFACE biom_check_data_output

!-- Input parameter checks to be done in check_parameters
    INTERFACE biom_check_parameters
       MODULE PROCEDURE biom_check_parameters
    END INTERFACE biom_check_parameters

!-- Data output of 2D quantities
    INTERFACE biom_data_output_2d
       MODULE PROCEDURE biom_data_output_2d
    END INTERFACE biom_data_output_2d

!-- no 3D data, thus, no averaging of 3D data, removed
    INTERFACE biom_data_output_3d
       MODULE PROCEDURE biom_data_output_3d
    END INTERFACE biom_data_output_3d

!-- Definition of data output quantities
    INTERFACE biom_define_netcdf_grid
       MODULE PROCEDURE biom_define_netcdf_grid
    END INTERFACE biom_define_netcdf_grid

!-- Obtains all relevant input values to estimate local thermal comfort/stress
    INTERFACE biom_determine_input_at
       MODULE PROCEDURE biom_determine_input_at
    END INTERFACE biom_determine_input_at

!-- Output of information to the header file
    INTERFACE biom_header
       MODULE PROCEDURE biom_header
    END INTERFACE biom_header

!-- Initialization actions 
    INTERFACE biom_init
       MODULE PROCEDURE biom_init
    END INTERFACE biom_init

!-- Initialization of arrays
    INTERFACE biom_init_arrays
       MODULE PROCEDURE biom_init_arrays
    END INTERFACE biom_init_arrays

!-- Reading of NAMELIST parameters
    INTERFACE biom_parin
       MODULE PROCEDURE biom_parin
    END INTERFACE biom_parin

 
 CONTAINS
 
 
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Sum up and time-average biom input quantities as well as allocate
!> the array necessary for storing the average.
!------------------------------------------------------------------------------!
 SUBROUTINE biom_3d_data_averaging( mode, variable )

    IMPLICIT NONE

    CHARACTER (LEN=*) ::  mode     !< 
    CHARACTER (LEN=*) ::  variable !< 

    INTEGER(iwp) ::  i  !< 
    INTEGER(iwp) ::  j  !< 
    INTEGER(iwp) ::  k  !< 


    IF ( mode == 'allocate' )  THEN

       SELECT CASE ( TRIM( variable ) )

          CASE ( 'biom_mrt' )
                IF ( .NOT. ALLOCATED( biom_mrt_av ) )  THEN
                   ALLOCATE( biom_mrt_av(nmrtbl) )
                ENDIF
                biom_mrt_av = 0.0_wp

          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )

!--          Indices in unknown order as depending on input file, determine 
!            first index to average und update only once
             IF ( .NOT. average_trigger_pt .AND. .NOT. average_trigger_utci    &
                .AND. .NOT. average_trigger_pet ) THEN
                IF ( TRIM( variable ) == 'biom_pt' ) THEN
                    average_trigger_pt = .TRUE.
                ENDIF
                IF ( TRIM( variable ) == 'biom_utci' ) THEN
                    average_trigger_utci = .TRUE.
                ENDIF
                IF ( TRIM( variable ) == 'biom_pet' ) THEN
                    average_trigger_pet = .TRUE.
                ENDIF
             ENDIF

!--          Only continue if updateindex
             IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')   RETURN
             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci') RETURN
             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')  RETURN

!--          Set averaging switch to .TRUE. if not done by other module before!
             IF ( .NOT. ALLOCATED( pt_av ) )  THEN
                ALLOCATE( pt_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
                aver_pt = .TRUE.
             ENDIF
             pt_av = 0.0_wp

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

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

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

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

          CASE DEFAULT
             CONTINUE

       END SELECT

    ELSEIF ( mode == 'sum' )  THEN

       SELECT CASE ( TRIM( variable ) )

          CASE ( 'biom_mrt' )
             IF ( ALLOCATED( biom_mrt_av ) )  THEN

                 IF ( nmrtbl > 0 )  THEN
                    IF ( mrt_include_sw )  THEN
                       biom_mrt_av(:) = biom_mrt_av(:) + &
                          ((human_absorb*mrtinsw(:) + human_emiss*mrtinlw(:))  &
                          / (human_emiss*sigma_sb)) ** .25_wp - degc_to_k
                    ELSE
                       biom_mrt_av(:) = biom_mrt_av(:) + &
                          (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp      &
                          - degc_to_k
                    ENDIF
                 ENDIF
             ENDIF

          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )

!--          Only continue if updateindex
             IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')    &
                RETURN
             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci')  &
                RETURN
             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')   &
                RETURN

             IF ( ALLOCATED( pt_av ) .AND. aver_pt ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

             IF ( ALLOCATED( q_av )  .AND. aver_q ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

             IF ( ALLOCATED( u_av )  .AND. aver_u ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

             IF ( ALLOCATED( v_av )  .AND. aver_v ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

             IF ( ALLOCATED( w_av )  .AND. aver_w ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

           CASE DEFAULT
             CONTINUE

       END SELECT

    ELSEIF ( mode == 'average' )  THEN

       SELECT CASE ( TRIM( variable ) )

          CASE ( 'biom_mrt' )
             IF ( ALLOCATED( biom_mrt_av ) )  THEN
                biom_mrt_av(:) = biom_mrt_av(:) / REAL( average_count_3d, KIND=wp )
             ENDIF

          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )

!--          Only continue if update index
             IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')    &
                RETURN
             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci')  &
                RETURN
             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')   &
                RETURN

             IF ( ALLOCATED( pt_av ) .AND. aver_pt ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

             IF ( ALLOCATED( q_av ) .AND. aver_q ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

             IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

             IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

             IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN
                DO  i = nxl, nxr
                   DO  j = nys, nyn
                      DO  k = nzb, nzt+1
                         w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND=wp )
                      ENDDO
                   ENDDO
                ENDDO
             ENDIF

!--          Udate thermal indices with derived averages
             CALL biom_calculate_static_grid ( .TRUE. )

        END SELECT

    ENDIF


 END SUBROUTINE biom_3d_data_averaging



!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check data output for biometeorology model
!------------------------------------------------------------------------------!
 SUBROUTINE biom_check_data_output( var, unit )

    USE control_parameters,                                                    &
        ONLY: data_output, message_string

    IMPLICIT NONE

    CHARACTER (LEN=*) ::  unit     !< The unit for the variable var
    CHARACTER (LEN=*) ::  var      !< The variable in question


    SELECT CASE ( TRIM( var ) )
     
       CASE( 'biom_tmrt', 'biom_mrt', 'biom_pet', 'biom_pt', 'biom_utci' )
          unit = 'degree_C'

       CASE DEFAULT
          unit = 'illegal'

    END SELECT

    IF ( unit /= 'illegal' )  THEN
!
!--    Futher checks if output belongs to biometeorology
       IF ( .NOT. biometeorology .OR. .NOT. radiation ) THEN
         message_string = 'output of "' // TRIM( var ) // '" req' //  &
               'uires biometeorology = .TRUE. and radiation = .TRUE. ' &
               // 'in initialisation_parameters'
         CALL message( 'biom_check_data_output', 'PA0561', 1, 2, 0, 6, 0 )
       ENDIF
       IF ( mrt_nlevels == 0 ) THEN
          message_string = 'output of "' // TRIM( var ) // '" require'&
                           // 's mrt_nlevels > 0'
          CALL message( 'biom_check_data_output', 'PA0562', 1, 2, 0, 6, 0 )
       ENDIF

    ENDIF

 END SUBROUTINE biom_check_data_output

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check parameters routine for biom module
!> check_parameters 1302
!------------------------------------------------------------------------------!
 SUBROUTINE biom_check_parameters

    USE control_parameters,                                                    &
        ONLY:  message_string


    IMPLICIT NONE


!--    Disabled as long as radiation model not available
       IF ( .NOT. radiation )  THEN
          message_string = 'The human thermal comfort module does require ' // &
                           'radiation information in terms of the mean ' //    &
                           'radiant temperature, but radiation is not enabled!'
          CALL message( 'check_parameters', 'PAHU01', 0, 0, 0, 6, 0 )
       ENDIF

       IF ( .NOT. humidity )  THEN
          message_string = 'The human thermal comfort module does require ' // &
                           'air humidity information, but humidity module ' // &
                           'is disabled!'
          CALL message( 'check_parameters', 'PAHU02', 0, 0, 0, 6, 0 )
       ENDIF

 END SUBROUTINE biom_check_parameters


!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Subroutine defining 3D output variables (dummy, always 2d!)
!> data_output_3d 709ff
!------------------------------------------------------------------------------!
 SUBROUTINE biom_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )

    USE indices

    USE kinds


    IMPLICIT NONE

!-- Input variables
    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.

!-- Output variables
    LOGICAL, INTENT(OUT) ::  found   !< Output found?
    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return

!-- Internal variables
    INTEGER(iwp) ::  l    !< Running index, radiation grid
    INTEGER(iwp) ::  i    !< Running index, x-dir
    INTEGER(iwp) ::  j    !< Running index, y-dir
    INTEGER(iwp) ::  k    !< Running index, z-dir

    CHARACTER (LEN=:), allocatable ::  variable_short  !< Trimmed version of char variable

    REAL(wp), PARAMETER ::  fill_value = -999._wp
    REAL(wp) ::  mrt  !< Buffer for mean radiant temperature

    found = .TRUE.
    variable_short = TRIM( variable )

    IF ( variable_short(1:4) /= 'biom' ) THEN
!--   Nothing to do, set found to FALSE and return immediatelly
      found = .FALSE.
      RETURN
    ENDIF

    SELECT CASE ( variable_short )

        CASE ( 'biom_mrt' )

            local_pf = REAL( fill_value, KIND = wp )
            DO  l = 1, nmrtbl
                i = mrtbl(ix,l)
                j = mrtbl(iy,l)
                k = mrtbl(iz,l)
                IF ( mrt_include_sw )  THEN
                    mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
                           / (human_emiss*sigma_sb)) ** .25_wp
                ELSE
                    mrt = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
                ENDIF
                local_pf(i,j,k) = mrt
            ENDDO

        CASE ( 'biom_tmrt' )        ! 2d-array
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb_do) = REAL( tmrt_grid(j,i), sp )
                 ENDDO
              ENDDO

        CASE ( 'biom_pt' )        ! 2d-array
            IF ( av == 0 )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb_do) = REAL( pt_grid(j,i), sp )
                 ENDDO
              ENDDO
            ELSE
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb_do) = REAL( pt_av_grid(j,i), sp )
                 ENDDO
              ENDDO
            END IF

        CASE ( 'biom_utci' )        ! 2d-array
            IF ( av == 0 )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb_do) = REAL( utci_grid(j,i), sp )
                 ENDDO
              ENDDO
            ELSE
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb_do) = REAL( utci_av_grid(j,i), sp )
                 ENDDO
              ENDDO
            END IF

        CASE ( 'biom_pet' )        ! 2d-array
            IF ( av == 0 )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb_do) = REAL( pet_grid(j,i), sp )
                 ENDDO
              ENDDO
            ELSE
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb_do) = REAL( pet_av_grid(j,i), sp )
                 ENDDO
              ENDDO
            END IF

       CASE DEFAULT
          found = .FALSE.

    END SELECT

 END SUBROUTINE biom_data_output_3d

!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Subroutine defining 2D output variables
!> data_output_2d 1188ff
!------------------------------------------------------------------------------!
 SUBROUTINE biom_data_output_2d( av, variable, found, grid, local_pf,          &
                                      two_d, nzb_do, nzt_do, fill_value )

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

    USE kinds


    IMPLICIT NONE

!-- Input variables
    CHARACTER (LEN=*), INTENT(IN) ::  variable    !< Char identifier to select var for output
    INTEGER(iwp), INTENT(IN)      ::  av          !< Use averaged data? 0 = no, 1 = yes?
    INTEGER(iwp), INTENT(IN)      ::  nzb_do      !< Unused. 2D. nz bottom to nz top
    INTEGER(iwp), INTENT(IN)      ::  nzt_do      !< Unused.
    REAL(wp), INTENT(in)          ::  fill_value  !< Fill value for unassigned locations

!-- Output variables
    CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
    LOGICAL, INTENT(OUT)           ::  found  !< Output found?
    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables, horizontal cross sections, must be .TRUE.
    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return

!-- Internal variables
    CHARACTER (LEN=:), allocatable ::  variable_short  !< Trimmed version of char variable
    INTEGER(iwp) ::  i        !< Running index, x-dir
    INTEGER(iwp) ::  j        !< Running index, y-dir
    INTEGER(iwp) ::  k        !< Running index, z-dir


    found = .TRUE.
    variable_short = TRIM( variable )
    IF ( variable_short(1:4) == 'biom' ) THEN
       two_d = .TRUE.
       grid = 'zu1'
    ELSE
       found = .FALSE.
       grid  = 'none'
       RETURN
    ENDIF

    local_pf = fill_value

    SELECT CASE ( variable_short )


        CASE ( 'biom_tmrt_xy' )        ! 2d-array
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,1) = tmrt_grid(j,i)
                 ENDDO
              ENDDO


        CASE ( 'biom_pt_xy' )        ! 2d-array
            IF ( av == 0 )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = pt_grid(j,i)
                 ENDDO
              ENDDO
            ELSE
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = pt_av_grid(j,i)
                 ENDDO
              ENDDO
            END IF


        CASE ( 'biom_utci_xy' )        ! 2d-array
            IF ( av == 0 )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = utci_grid(j,i)
                 ENDDO
              ENDDO
            ELSE
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = utci_av_grid(j,i)
                 ENDDO
              ENDDO
            END IF


        CASE ( 'biom_pet_xy' )        ! 2d-array
            IF ( av == 0 )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = pet_grid(j,i)
                 ENDDO
              ENDDO
            ELSE
              DO  i = nxl, nxr
                 DO  j = nys, nyn
                    local_pf(i,j,nzb+1) = pet_av_grid(j,i)
                 ENDDO
              ENDDO
            END IF


       CASE DEFAULT
          found = .FALSE.
          grid  = 'none'

    END SELECT


 END SUBROUTINE biom_data_output_2d


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Subroutine defining appropriate grid for netcdf variables.
!> It is called out from subroutine netcdf_interface_mod.
!> netcdf_interface_mod 918ff
!------------------------------------------------------------------------------!
 SUBROUTINE biom_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )

    IMPLICIT NONE

!-- Input variables
    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< Name of output variable

!-- Output variables
    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable

    LOGICAL, INTENT(OUT)           ::  found    !< Flag if output var is found

!-- Local variables
    LOGICAL      :: is2d  !< Var is 2d?

    INTEGER(iwp) :: l     !< Length of the var array


    found  = .FALSE.
    grid_x = 'none'
    grid_y = 'none'
    grid_z = 'none'

    l = MAX( 2, LEN_TRIM( var ) )
    is2d = ( var(l-1:l) == 'xy' )


    IF ( var(1:4) == 'biom' ) THEN
       found  = .TRUE.
       grid_x = 'x'
       grid_y = 'y'
       grid_z = 'zu'
       IF ( is2d ) grid_z = 'zu1'
    ENDIF

 END SUBROUTINE biom_define_netcdf_grid

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Header output for biom module
!> header 982
!------------------------------------------------------------------------------!
 SUBROUTINE biom_header( io )

    IMPLICIT NONE

!-- Input variables
    INTEGER(iwp), INTENT(IN) ::  io           !< Unit of the output file

!-- Internal variables
    CHARACTER (LEN=86) ::  output_height_chr  !< String for output height

    WRITE( output_height_chr, '(F8.1,7X)' )  biom_output_height
!
!-- Write biom header
    WRITE( io, 1 )
    WRITE( io, 2 )  TRIM( output_height_chr )
    WRITE( io, 3 )  TRIM( ACHAR( biom_cell_level ) )

1   FORMAT (//' Human thermal comfort module information:'/                    &
              ' ------------------------------'/)
2   FORMAT ('    --> All indices calculated for a height of (m): ', A )
3   FORMAT ('    --> This corresponds to cell level : ', A )

 END SUBROUTINE biom_header


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Initialization of the HTCM
!> init_3d_model 1987ff
!------------------------------------------------------------------------------!
 SUBROUTINE biom_init

    USE control_parameters,                                                 &
        ONLY: message_string

    IMPLICIT NONE

!-- Internal vriables
    REAL ( wp )  :: height  !< current height in meters

    INTEGER ( iwp )  :: i  !< iteration index

!-- Determine cell level corresponding to 1.1 m above ground level
!   (gravimetric center of sample human)

    time_biom_results = 0.0_wp
    biom_cell_level = 0_iwp
    biom_output_height = 0.5_wp * dz(1)
    height = 0.0_wp

    biom_cell_level = INT ( 1.099_wp / dz(1) )
    biom_output_height = biom_output_height + biom_cell_level * dz(1)

    IF ( .NOT. radiation_interactions .AND. mrt_from_rtm ) THEN
       message_string = 'The mrt_from_rtm switch require ' // &
                        'enabled radiation_interactions but it ' // &
                        'is disabled! Set mrt_from_rtm to .F.'
       CALL message( 'check_parameters', 'PAHU03', 0, 0, -1, 6, 0 )
       mrt_from_rtm = .FALSE.
    ENDIF

 END SUBROUTINE biom_init


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Allocate biom arrays and define pointers if required
!> init_3d_model 1047ff
!------------------------------------------------------------------------------!
 SUBROUTINE biom_init_arrays

    IMPLICIT NONE

!-- Allocate a temporary array with the desired output dimensions.
!   Initialization omitted for performance, will be overwritten anyway
    IF ( .NOT. ALLOCATED( tmrt_grid ) ) THEN
      ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
    ENDIF

    IF ( biom_pt ) THEN
      IF ( .NOT. ALLOCATED( pt_grid ) ) THEN
        ALLOCATE( pt_grid (nys:nyn,nxl:nxr) )
      ENDIF
    ENDIF

    IF ( biom_utci ) THEN
      IF ( .NOT. ALLOCATED( utci_grid ) ) THEN
        ALLOCATE( utci_grid (nys:nyn,nxl:nxr) )
      ENDIF
    ENDIF

    IF ( biom_pet ) THEN
      IF ( .NOT. ALLOCATED( pet_grid ) ) THEN
        ALLOCATE( pet_grid (nys:nyn,nxl:nxr) )
      ENDIF
    END IF

    IF ( biom_pt_av ) THEN
      IF ( .NOT. ALLOCATED( pt_av_grid ) ) THEN
        ALLOCATE( pt_av_grid (nys:nyn,nxl:nxr) )
      ENDIF
    ENDIF

    IF ( biom_utci_av ) THEN
      IF ( .NOT. ALLOCATED( utci_av_grid ) ) THEN
        ALLOCATE( utci_av_grid (nys:nyn,nxl:nxr) )
      ENDIF
    ENDIF

    IF ( biom_pet_av ) THEN
      IF ( .NOT. ALLOCATED( pet_av_grid ) ) THEN
        ALLOCATE( pet_av_grid (nys:nyn,nxl:nxr) )
      ENDIF

    END IF

 END SUBROUTINE biom_init_arrays


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Parin for &biometeorology_parameters for reading biomet parameters
!------------------------------------------------------------------------------!
 SUBROUTINE biom_parin

    IMPLICIT NONE

!
!-- Internal variables
    CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file 

    NAMELIST /biometeorology_parameters/  mrt_from_rtm,                        &
                                          biom_pet,                            &
                                          biom_pet_av,                         &
                                          biom_pt,                             &
                                          biom_pt_av,                          &
                                          biom_utci,                           &
                                          biom_utci_av


!-- Try to find biometeorology_parameters namelist
    REWIND ( 11 )
    line = ' '
    DO WHILE ( INDEX( line, '&biometeorology_parameters' ) == 0 )
       READ ( 11, '(A)', END = 20 )  line
    ENDDO
    BACKSPACE ( 11 )

!
!-- Read biometeorology_parameters namelist
    READ ( 11, biometeorology_parameters, ERR = 10, END = 20 )

!
!-- Set flag that indicates that the biomet_module is switched on
    biometeorology = .TRUE.

    GOTO 20

!
!-- In case of error
 10 BACKSPACE( 11 )
    READ( 11 , '(A)') line
    CALL parin_fail_message( 'biometeorology_parameters', line )

!
!-- Complete
 20 CONTINUE


 END SUBROUTINE biom_parin

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculates the mean radiant temperature (tmrt) based on the Six-directions
!> method according to VDI 3787 2.
!------------------------------------------------------------------------------!
 SUBROUTINE calculate_tmrt_6_directions( SW_N, SW_E, SW_S, SW_W,               &
    SW_U, SW_D, LW_N, LW_E, LW_S, LW_W, LW_U, LW_D, tmrt )

    IMPLICIT NONE

!-- Type of input of the argument list
!   Short- (SW_) and longwave (LW_) radiation fluxes from the six directions
!   North (N), East (E), South (S), West (W), up (U) and down (D)
    REAL(wp), INTENT ( IN )  ::  SW_N  !< Sw radflux density from N    (W/m²)
    REAL(wp), INTENT ( IN )  ::  SW_E  !< Sw radflux density from E    (W/m²)
    REAL(wp), INTENT ( IN )  ::  SW_S  !< Sw radflux density from S    (W/m²)
    REAL(wp), INTENT ( IN )  ::  SW_W  !< Sw radflux density from W    (W/m²)
    REAL(wp), INTENT ( IN )  ::  SW_U  !< Sw radflux density from U    (W/m²)
    REAL(wp), INTENT ( IN )  ::  SW_D  !< Sw radflux density from D    (W/m²)
    REAL(wp), INTENT ( IN )  ::  LW_N  !< Lw radflux density from N    (W/m²)
    REAL(wp), INTENT ( IN )  ::  LW_E  !< Lw radflux density from E    (W/m²)
    REAL(wp), INTENT ( IN )  ::  LW_S  !< Lw radflux density from S    (W/m²)
    REAL(wp), INTENT ( IN )  ::  LW_W  !< Lw radflux density from W    (W/m²)
    REAL(wp), INTENT ( IN )  ::  LW_U  !< Lw radflux density from U    (W/m²)
    REAL(wp), INTENT ( IN )  ::  LW_D  !< Lw radflux density from D    (W/m²)

!-- Type of output of the argument list
    REAL(wp), INTENT ( OUT ) ::  tmrt  !< Mean radiant temperature (°C)

!-- Directional weighting factors
    REAL(wp), PARAMETER      ::  weight_h  = 0.22_wp
    REAL(wp), PARAMETER      ::  weight_v  = 0.06_wp

    REAL(wp) ::  nrfd           !<  Net radiation flux density (W/m²)

!-- Initialise
    nrfd = 0._wp

!-- Compute mean radiation flux density absorbed by rotational symetric human
    nrfd = ( weight_h * ( human_absorb * SW_N + human_emiss * LW_N ) ) +       &
       ( weight_h * ( human_absorb * SW_E + human_emiss * LW_E ) ) +           &
       ( weight_h * ( human_absorb * SW_S + human_emiss * LW_S ) ) +           &
       ( weight_h * ( human_absorb * SW_W + human_emiss * LW_W ) ) +           &
       ( weight_v * ( human_absorb * SW_U + human_emiss * LW_U ) ) +           &
       ( weight_v * ( human_absorb * SW_D + human_emiss * LW_D ) )

!-- Compute mean radiant temperature
    tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - degc_to_k

 END SUBROUTINE calculate_tmrt_6_directions

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Very crude approximation of mean radiant temperature based on upwards and 
!> downwards radiation fluxes only (other directions curently not available, 
!> replace as soon as possible!)
!------------------------------------------------------------------------------!
 SUBROUTINE calculate_tmrt_2_directions( sw_u, sw_d, lw_u, lw_d, ta, tmrt )

    IMPLICIT NONE

!-- Type of input of the argument list
    REAL(wp), INTENT ( IN )  ::  sw_u  !< Shortwave radiation flux density from upper direction (W/m²)
    REAL(wp), INTENT ( IN )  ::  sw_d  !< Shortwave radiation flux density from lower direction (W/m²)
    REAL(wp), INTENT ( IN )  ::  lw_u  !< Longwave radiation flux density from upper direction (W/m²)
    REAL(wp), INTENT ( IN )  ::  lw_d  !< Longwave radiation flux density from lower direction (W/m²)
    REAL(wp), INTENT ( IN )  ::  ta    !< Air temperature (°C)

!-- Type of output of the argument list
    REAL(wp), INTENT ( OUT ) ::  tmrt  !< mean radiant temperature, (°C)

!-- Directional weighting factors and parameters
    REAL(wp), PARAMETER ::  weight_h  = 0.22_wp     !< Weight for horizontal radiational gain after Fanger (1972)
    REAL(wp), PARAMETER ::  weight_v  = 0.06_wp     !< Weight for vertical radiational gain after Fanger (1972)

!-- Other internal variables
    REAL(wp) ::  sw_in
    REAL(wp) ::  sw_out
    REAL(wp) ::  lw_in
    REAL(wp) ::  lw_out
    REAL(wp) ::  nrfd     !< Net radiation flux density (W/m²)
    REAL(wp) ::  lw_air   !< Longwave emission by surrounding air volume (W/m²)
    REAL(wp) ::  sw_side  !< Shortwave immission from the sides (W/m²)

    INTEGER(iwp) ::  no_input  !< Count missing input radiation fluxes

!-- initialise
    sw_in    = sw_u
    sw_out   = sw_d
    lw_in    = lw_u
    lw_out   = lw_d
    nrfd     = 0._wp
    no_input = 0_iwp

!-- test for missing input data
    IF ( sw_in <= -998._wp .OR. sw_out <= -998._wp .OR. lw_in <= -998._wp .OR. &
       lw_out <= -998._wp .OR. ta <= -998._wp ) THEN
       IF ( sw_in <= -998._wp ) THEN
          sw_in = 0.
          no_input = no_input + 1
       ENDIF
       IF ( sw_out <= -998._wp ) THEN
          sw_out = 0.
          no_input = no_input + 1
       ENDIF
       IF ( lw_in <= -998._wp ) THEN
          lw_in = 0.
          no_input = no_input + 1
       ENDIF
       IF ( lw_out <= -998._wp ) THEN
          lw_out = 0.
          no_input = no_input + 1
       ENDIF

!-- Accept two missing radiation flux directions, fail otherwise as error might become too large
       IF ( ta <= -998._wp .OR. no_input >= 3 ) THEN
          tmrt = -999._wp
          RETURN
       ENDIF
    ENDIF

    sw_side = sw_in * 0.125_wp  ! distribute half of upper sw_in to the 4 sides
    lw_air = ( sigma_sb * 0.95_wp * ( ta + degc_to_k )**4 )

!-- Compute mean radiation flux density absorbed by rotational symetric human
    nrfd = ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +         &
       ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
       ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
       ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
       ( weight_v * ( human_absorb * (sw_in * 0.5_wp) + human_emiss * lw_in ) ) +     &
       ( weight_v * ( human_absorb * sw_out + human_emiss * lw_out ) )

!-- Compute mean radiant temperature
    tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - degc_to_k

 END SUBROUTINE calculate_tmrt_2_directions

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate static thermal indices for given meteorological conditions
!------------------------------------------------------------------------------!
 SUBROUTINE calculate_static_thermal_indices ( ta, vp, ws, pair, tmrt,         &
    pt_static, utci_static, pet_static )

    IMPLICIT NONE

!-- Input parameters
    REAL(wp), INTENT ( IN )  ::  ta           !< Air temperature                  (°C)
    REAL(wp), INTENT ( IN )  ::  vp           !< Vapour pressure                  (hPa)
    REAL(wp), INTENT ( IN )  ::  ws           !< Wind speed    (local level)      (m/s)
    REAL(wp), INTENT ( IN )  ::  pair         !< Air pressure                     (hPa)
    REAL(wp), INTENT ( IN )  ::  tmrt         !< Mean radiant temperature         (°C)
!-- Output parameters
    REAL(wp), INTENT ( OUT ) ::  pt_static    !< Perceived temperature            (°C)
    REAL(wp), INTENT ( OUT ) ::  utci_static  !< Universal thermal climate index  (°C)
    REAL(wp), INTENT ( OUT ) ::  pet_static   !< Physiologically equivalent temp. (°C)
!-- Temporary field, not used here
    REAL(wp)                 ::  clo          !< Clothing index (no dim.)

    clo = -999._wp

    IF ( biom_pt ) THEN
!-- Estimate local perceived temperature
       CALL calculate_pt_static( ta, vp, ws, tmrt, pair, clo, pt_static )
    ENDIF

    IF ( biom_utci ) THEN
!-- Estimate local universal thermal climate index
       CALL calculate_utci_static( ta, vp, ws, tmrt, biom_output_height,       &
          utci_static )
    ENDIF

    IF ( biom_pet ) THEN
!-- Estimate local physiologically equivalent temperature
       CALL calculate_pet_static( ta, vp, ws, tmrt, pair, pet_static )
    ENDIF

 END SUBROUTINE calculate_static_thermal_indices


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate static thermal indices for 2D grid point i, j
!------------------------------------------------------------------------------!
 SUBROUTINE biom_determine_input_at( average_input, i, j, ta, vp, ws, pair,    &
    tmrt )

    IMPLICIT NONE

!-- Input variables
    LOGICAL,      INTENT ( IN ) ::  average_input  !< Determine averaged input conditions?
    INTEGER(iwp), INTENT ( IN ) ::  i     !< Running index, x-dir
    INTEGER(iwp), INTENT ( IN ) ::  j     !< Running index, y-dir

!-- Output parameters
    REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (°C)
    REAL(wp), INTENT ( OUT )    ::  ta    !< Air temperature                 (°C)
    REAL(wp), INTENT ( OUT )    ::  vp    !< Vapour pressure                 (hPa)
    REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
    REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)

!-- Internal variables
    INTEGER(iwp)                ::  k     !< Running index, z-dir
    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only

    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)


!-- Determine cell level closest to 1.1m above ground
!   by making use of truncation due to int cast
    k = get_topography_top_index_ji(j, i, 's') + biom_cell_level  !< Vertical cell center closest to 1.1m
    k_wind = k
    IF( k_wind < 1_iwp ) THEN  ! Avoid horizontal u and v of 0.0 m/s close to ground
       k_wind = 1_iwp
    ENDIF

!-- Determine local values:
    IF ( average_input ) THEN
!--    Calculate ta from Tp assuming dry adiabatic laps rate
       ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k

       vp = 0.034_wp
       IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
          vp = q_av(k, j, i)
       ENDIF

       ws = ( 0.5_wp * ABS( u_av(k_wind, j, i) + u_av(k_wind, j, i+1) )  +     &
          0.5_wp * ABS( v_av(k_wind, j, i) + v_av(k_wind, j+1, i) )  +         &
          0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) )
    ELSE
!-- Calculate ta from Tp assuming dry adiabatic laps rate
       ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k

       vp = q(k, j, i)

       ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) )  +           &
          0.5_wp * ABS( v(k_wind, j, i) + v(k_wind, j+1, i) )  +               &
          0.5_wp * ABS( w(k_wind, j, i) + w(k_wind+1, j, i) ) )

    ENDIF

!-- Local air pressure
    pair = surface_pressure
!
!-- Calculate water vapour pressure at saturation and convert to hPa
!-- The magnus formula is limited to temperatures up to 333.15 K to
!   avoid negative values of vp_sat
    vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
    vp  = vp * pair / ( vp + 0.622_wp )
    IF ( vp > vp_sat ) vp = vp_sat

    tmrt = ta
    IF ( radiation ) THEN
       IF ( mrt_from_rtm ) THEN
          tmrt = tmrt_grid(j, i)
       ELSE
          CALL calculate_tmrt_2_directions (rad_sw_in(0, j, i),                &
             rad_sw_out(0, j, i), rad_lw_in(0, j, i), rad_lw_out(0, j, i), ta, &
             tmrt )
       ENDIF
    ENDIF

 END SUBROUTINE biom_determine_input_at


!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate static thermal indices for any point within a 2D grid
!> time_integration.f90: 1065ff
!------------------------------------------------------------------------------!
 SUBROUTINE biom_calculate_static_grid ( average_input )

    IMPLICIT NONE

!-- Input attributes
    LOGICAL, INTENT ( IN ) ::  average_input  !< Calculate based on averaged input? conditions?

!-- Internal variables
    INTEGER(iwp) ::  i, j, k, l   !< Running index

    REAL(wp) ::  ta            !< Air temperature                  (°C)
    REAL(wp) ::  vp            !< Vapour pressure                  (hPa)
    REAL(wp) ::  ws            !< Wind speed    (local level)      (m/s)
    REAL(wp) ::  pair          !< Air pressure                     (hPa)
    REAL(wp) ::  tmrt_tmp      !< Mean radiant temperature
    REAL(wp) ::  pt_tmp        !< Perceived temperature
    REAL(wp) ::  utci_tmp      !< Universal thermal climate index
    REAL(wp) ::  pet_tmp       !< Physiologically equivalent temperature


    CALL biom_init_arrays

    IF ( mrt_from_rtm ) THEN
       tmrt_grid = -999._wp
       DO  l = 1, nmrtbl
           i = mrtbl(ix,l)
           j = mrtbl(iy,l)
           k = mrtbl(iz,l)
           IF ( k - get_topography_top_index_ji( j, i, 's' ) == 1  ) THEN
              IF ( mrt_include_sw )  THEN
                  tmrt_tmp = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
                             / (human_emiss*sigma_sb)) ** .25_wp
              ELSE
                  tmrt_tmp = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
              ENDIF
              tmrt_grid(j, i) = tmrt_tmp - degc_to_k
           ENDIF
       ENDDO
    ENDIF

    DO i = nxl, nxr
       DO j = nys, nyn
!--       Determine local input conditions
          CALL biom_determine_input_at ( average_input, i, j, ta, vp, ws,      &
               pair, tmrt_tmp )
          tmrt_grid(j, i) = tmrt_tmp

!--       Only proceed if tmrt is available
          IF ( tmrt_tmp <= -998._wp ) THEN
             pt_tmp   = -999._wp
             utci_tmp = -999._wp
             pet_tmp  = -999._wp
             CYCLE
          END IF

!--       Calculate static thermal indices based on local tmrt
          CALL calculate_static_thermal_indices ( ta, vp, ws,                  &
             pair, tmrt_tmp, pt_tmp, utci_tmp, pet_tmp )

          IF ( average_input ) THEN
!--          Write results for selected averaged indices
             IF ( biom_pt_av )  THEN
                pt_av_grid(j, i)   = pt_tmp
             END IF
             IF ( biom_utci_av ) THEN
                utci_av_grid(j, i) = utci_tmp
             END IF
             IF ( biom_pet_av ) THEN
                pet_av_grid(j, i)  = pet_tmp
             END IF
          ELSE
!--          Write result for selected indices
             IF ( biom_pt )  THEN
                pt_grid(j, i)   = pt_tmp
             END IF
             IF ( biom_utci ) THEN
                utci_grid(j, i) = utci_tmp
             END IF
             IF ( biom_pet ) THEN
                pet_grid(j, i)  = pet_tmp
             END IF
          END IF

       END DO
    END DO

 END SUBROUTINE biom_calculate_static_grid

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate dynamic thermal indices (currently only iPT, but expandable)
!------------------------------------------------------------------------------!
 SUBROUTINE biom_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,         &
    t_clo, clo, actlev, age, weight, height, work, sex, ipt )

    IMPLICIT NONE

!-- Input parameters
    REAL(wp), INTENT ( IN )  ::  ta   !< Air temperature                  (°C)
    REAL(wp), INTENT ( IN )  ::  vp   !< Vapour pressure                  (hPa)
    REAL(wp), INTENT ( IN )  ::  ws   !< Wind speed    (local level)      (m/s)
    REAL(wp), INTENT ( IN )  ::  pair !< Air pressure                     (hPa)
    REAL(wp), INTENT ( IN )  ::  tmrt !< Mean radiant temperature         (°C)
    REAL(wp), INTENT ( IN )  ::  dt   !< Time past since last calculation (s)
    REAL(wp), INTENT ( IN )  ::  age  !< Age of agent                     (y)
    REAL(wp), INTENT ( IN )  ::  weight  !< Weight of agent               (Kg)
    REAL(wp), INTENT ( IN )  ::  height  !< Height of agent               (m)
    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical workload of agent
                                         !  (without metabolism!)         (W)
    INTEGER(iwp), INTENT ( IN ) ::  sex  !< Sex of agent (1 = male, 2 = female)

!-- Both, input and output parameters
    Real(wp), INTENT ( INOUT )  ::  energy_storage    !< Energy storage   (W/m²)
    Real(wp), INTENT ( INOUT )  ::  t_clo   !< Clothing temperature       (°C)
    Real(wp), INTENT ( INOUT )  ::  clo     !< Current clothing in sulation
    Real(wp), INTENT ( INOUT )  ::  actlev  !< Individuals activity level
                                            !  per unit surface area      (W/m²)
!-- Output parameters
    REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (°C)

!-- If clo equals the initial value, this is the initial call
    IF ( clo <= -998._wp ) THEN
!--    Initialize instationary perceived temperature with personalized 
!      PT as an initial guess, set actlev and clo
       CALL ipt_init ( age, weight, height, sex, work, actlev, clo,            &
          ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,                   &
          ipt )
    ELSE
!--    Estimate local instatinoary perceived temperature
       CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,     &
          clo, actlev, work, ipt )
    ENDIF

 END SUBROUTINE biom_calc_ipt

 END MODULE biometeorology_mod