Index: palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_mod.f90 (revision 3361)
+++ palm/trunk/SOURCE/biometeorology_mod.f90 (revision 3448)
@@ -1,4 +1,4 @@
-!> @file biometeorology.f90
-!------------------------------------------------------------------------------!
+!> @file biometeorology_mod.f90
+!--------------------------------------------------------------------------------!
! This file is part of PALM-4U.
!
@@ -15,10 +15,8 @@
! PALM. If not, see .
!
-! Copyright 2018, Institute of Computer Science,Academy of Sciences, Prague
-!
-! Calculation of PET:
-! Copyright 2016, Deutscher Wetterdienst (DWD) /
-! German Meteorological Service (DWD)
-!------------------------------------------------------------------------------!
+! 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:
@@ -35,83 +33,446 @@
! Authors:
! --------
+! @author Dominik Froehlich
! @author Jaroslav Resler
-! @author Dominik Froehlich
-!
-!------------------------------------------------------------------------------!
-
-MODULE biometeorology_mod
-!
-!-- Load required variables from existing modules
+!
+!
+! 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: hyp, p, pt, q, u, v, w
+ 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: rd_d_rv
-
- USE kinds !< to set precision of INTEGER and REAL arrays according to PALM
-
- USE radiation_model_mod, &
- ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw, &
- mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation
-
-
- IMPLICIT NONE
-
- REAL(wp), DIMENSION(:), ALLOCATABLE :: bio_mrt !< biometeorology mean radiant temperature for each MRT box
- REAL(wp), DIMENSION(:), ALLOCATABLE :: bio_mrt_av !< time average mean radiant temperature for each MRT box
- REAL(wp), DIMENSION(:), ALLOCATABLE :: bio_pet !< PhysiologiCALLy Equivalent Temperature (PET) for each MRT box
- REAL(wp), DIMENSION(:), ALLOCATABLE :: bio_pet_av !< time average PhysiologiCALLy Equivalent Temperature (PET) for each MRT box
-
-
- REAL(wp), PARAMETER :: sigma_sb = 5.67037321E-8_wp, & !< Stefan-Boltzmann constant
- t_zero = -273.15_wp, & !< temperature 0K in Celsius
- human_absorb = 0.7_wp, & !< SW absorbtivity of human body
- human_emiss = 0.97_wp !< emissivity of human body (Lindberg 2008)
-
-
- INTERFACE biometeorology_3d_data_averaging
- MODULE PROCEDURE biometeorology_3d_data_averaging
- END INTERFACE biometeorology_3d_data_averaging
-
- INTERFACE biometeorology_check_data_output
- MODULE PROCEDURE biometeorology_check_data_output
- END INTERFACE biometeorology_check_data_output
-
- INTERFACE biometeorology_data_output_3d
- MODULE PROCEDURE biometeorology_data_output_3d
- END INTERFACE biometeorology_data_output_3d
-
- INTERFACE biometeorology_define_netcdf_grid
- MODULE PROCEDURE biometeorology_define_netcdf_grid
- END INTERFACE biometeorology_define_netcdf_grid
-
-
- SAVE
+ ONLY: magnus
+
+ 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, rad_sw_in, &
+ rad_sw_out, rad_lw_in, rad_lw_out
+
+ USE surface_mod, &
+ ONLY: get_topography_top_index_ji
+
+ IMPLICIT NONE
PRIVATE
-!
-!-- Public functions
- PUBLIC biometeorology_3d_data_averaging, biometeorology_check_data_output, &
- biometeorology_data_output_3d, biometeorology_define_netcdf_grid
-
-!
-!-- Public variables and constants / NEEDS SORTING
-! PUBLIC
-
-
+!-- 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 :: cels_offs = 273.15_wp !< Kelvin-Celsius offset (K)
+ REAL ( wp ), PARAMETER :: sigma_sb = 5.67037321E-8_wp !< Stefan-Boltzmann constant
+ 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 :: 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:
-! ------------
-!> Check data output for biometeorology model
-!------------------------------------------------------------------------------!
- SUBROUTINE biometeorology_check_data_output( var, unit, i, ilen, k )
+!------------------------------------------------------------------------------!
+! 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 - cels_offs
+ ELSE
+ biom_mrt_av(:) = biom_mrt_av(:) + &
+ (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp &
+ - cels_offs
+ 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
@@ -119,48 +480,73 @@
IMPLICIT NONE
- CHARACTER (LEN=*) :: unit !<
- CHARACTER (LEN=*) :: var !<
-
- INTEGER(iwp) :: i
- INTEGER(iwp) :: ilen
- INTEGER(iwp) :: k
-
- SELECT CASE ( TRIM( var ) )
-
- CASE ( 'bio_mrt', 'bio_pet' )
- IF ( .NOT. radiation ) THEN
- message_string = 'output of "' // TRIM( var ) // '" require'&
- // 's radiation = .TRUE.'
- CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
- ENDIF
- IF ( mrt_nlevels == 0 ) THEN
- message_string = 'output of "' // TRIM( var ) // '" require'&
- // 's mrt_nlevels > 0'
- CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
- ENDIF
- IF ( TRIM( var ) == 'bio_mrt' ) THEN
- unit = 'K'
- ELSE
- unit = 'W m-2'
- ENDIF
-
- CASE DEFAULT
- unit = 'illegal'
-
- END SELECT
-
- END SUBROUTINE biometeorology_check_data_output
-
-
-
-!------------------------------------------------------------------------------!
-!
-! Description:
-! ------------
-!> Subroutine for averaging 3D data
-!------------------------------------------------------------------------------!
-SUBROUTINE biometeorology_3d_data_averaging( mode, variable )
-
- USE control_parameters
+ CHARACTER (LEN=*) :: unit !< The unit for the variable var
+ CHARACTER (LEN=*) :: var !< The variable in question
+
+
+ unit = '°C'
+ IF ( .NOT. biometeorology ) THEN
+ message_string = 'output of "' // TRIM( var ) // '" req' // &
+ 'uires biometeorology = .TRUE. in initialisati' &
+ // 'on_parameters'
+ CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
+ unit = 'illegal'
+ ENDIF
+ IF ( .NOT. radiation ) THEN
+ message_string = 'output of "' // TRIM( var ) // '" require'&
+ // 's radiation = .TRUE.'
+ CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+ unit = 'illegal'
+ ENDIF
+ IF ( mrt_nlevels == 0 ) THEN
+ message_string = 'output of "' // TRIM( var ) // '" require'&
+ // 's mrt_nlevels > 0'
+ CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
+ unit = 'illegal'
+ 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
@@ -168,187 +554,40 @@
USE kinds
- IMPLICIT NONE
-
- CHARACTER (LEN=*) :: mode !<
- CHARACTER (LEN=*) :: variable !<
-
- INTEGER(iwp) :: i !<
- INTEGER(iwp) :: j !<
- INTEGER(iwp) :: k !<
- INTEGER(iwp) :: l, m !< index of current surface element
-
- REAL(wp) :: mrt, pet
-
- IF ( mode == 'allocate' ) THEN
-
- SELECT CASE ( TRIM( variable ) )
- CASE ( 'bio_mrt' )
- IF ( .NOT. ALLOCATED( bio_mrt_av ) ) THEN
- ALLOCATE( bio_mrt_av(nmrtbl) )
- ENDIF
- bio_mrt_av = 0.0_wp
-
- CASE ( 'bio_pet' )
- IF ( .NOT. ALLOCATED( bio_pet_av ) ) THEN
- ALLOCATE( bio_pet_av(nmrtbl) )
- ENDIF
- bio_pet_av = 0.0_wp
-
- CASE DEFAULT
- CONTINUE
-
- END SELECT
-
- ELSEIF ( mode == 'sum' ) THEN
-
- SELECT CASE ( TRIM( variable ) )
-
- CASE ( 'bio_mrt' )
- IF ( ALLOCATED( bio_mrt_av ) ) THEN
-
- IF ( nmrtbl > 0 ) THEN
- IF ( mrt_include_sw ) THEN
- bio_mrt_av(:) = bio_mrt_av(:) + ((human_absorb*mrtinsw(:) &
- + human_emiss*mrtinlw(:)) / (human_emiss*sigma_sb)) ** .25_wp
- ELSE
- bio_mrt_av(:) = bio_mrt_av(:) + &
- (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp
- ENDIF
- ENDIF
- ENDIF
-
- CASE ( 'bio_pet' )
- IF ( ALLOCATED( bio_pet_av ) ) THEN
- 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 = mrt + (human_emiss * mrtinlw(l) / sigma_sb) ** .25_wp
- ENDIF
-
- CALL calculate_pet_static( &
- pt(k,j,i) * (hyp(k) / 100000.0_wp )**0.286_wp + t_zero, & !< Air temperature (°C)
- q(k,j,i) * hyp(k) / ( q(k,j,i) + rd_d_rv ) / 100._wp, & !< Vapor pressure (hPa)
- SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
- ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
- ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2, &
- 0.01_wp ) ), & !< Wind speed (at scalar gridpoint) (m/s)
- mrt + t_zero, & !< Mean radiant temperature (°C)
- (hyp(k) + p(k,j,i)) * 0.01_wp, & !< Air pressure (hPa)
- pet ) !< output variable of PET
-
- bio_pet_av(l) = bio_pet_av(l) + pet
- ENDDO
- ENDIF
-
- CASE DEFAULT
- CONTINUE
-
- END SELECT
-
- ELSEIF ( mode == 'average' ) THEN
-
- SELECT CASE ( TRIM( variable ) )
-
- CASE ( 'bio_mrt' )
- IF ( ALLOCATED( bio_mrt_av ) ) THEN
- bio_mrt_av(:) = bio_mrt_av(:) / REAL( average_count_3d, KIND=wp )
- ENDIF
-
- CASE ( 'bio_pet' )
- IF ( ALLOCATED( bio_pet_av ) ) THEN
- bio_pet_av(:) = bio_pet_av(:) / REAL( average_count_3d, KIND=wp )
- ENDIF
-
- END SELECT
-
- ENDIF
-
-END SUBROUTINE biometeorology_3d_data_averaging
-
-
-!------------------------------------------------------------------------------!
-!
-! Description:
-! ------------
-!> Subroutine defining appropriate grid for netcdf variables.
-!> It is called out from subroutine netcdf.
-!------------------------------------------------------------------------------!
- SUBROUTINE biometeorology_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
-
- IMPLICIT NONE
-
- CHARACTER (LEN=*), INTENT(IN) :: var !<
- LOGICAL, INTENT(OUT) :: found !<
- CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
- CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
- CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
-
- found = .TRUE.
-
-
-!
-!-- Check for the grid
- SELECT CASE ( TRIM( var ) )
-
- CASE ( 'bio_mrt', 'bio_pet' )
- grid_x = 'x'
- grid_y = 'y'
- grid_z = 'zu'
-
- CASE DEFAULT
- found = .FALSE.
- grid_x = 'none'
- grid_y = 'none'
- grid_z = 'none'
-
- END SELECT
-
- END SUBROUTINE biometeorology_define_netcdf_grid
-
-
-!------------------------------------------------------------------------------!
-!
-! Description:
-! ------------
-!> Subroutine defining 3D output variables
-!------------------------------------------------------------------------------!
- SUBROUTINE biometeorology_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
-
-
- USE indices
-
- USE kinds
-
-
- IMPLICIT NONE
-
- CHARACTER (LEN=*) :: variable !<
-
- INTEGER(iwp) :: av !<
- INTEGER(iwp) :: i, j, k, l !<
- INTEGER(iwp) :: nzb_do !<
- INTEGER(iwp) :: nzt_do !<
-
- LOGICAL :: found !<
-
- REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
-
- REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
-
- REAL(wp) :: mrt, pet
+
+ 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.
-
-
- SELECT CASE ( TRIM( variable ) )
-
-
- CASE ( 'bio_mrt' )
+ 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 )
@@ -366,43 +605,55 @@
ENDDO
- CASE ( 'bio_pet' )
- local_pf = REAL( fill_value, KIND = wp )
+ 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 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 = mrt + (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
- ENDIF
-
- CALL calculate_pet_static( &
- pt(k,j,i) * (hyp(k) / 100000.0_wp )**0.286_wp + t_zero, & !< Air temperature (°C)
- q(k,j,i) * hyp(k) / ( q(k,j,i) + rd_d_rv ) / 100.0_wp, & !< Vapor pressure (hPa)
- SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
- ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
- ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2, &
- 0.01_wp ) ), & !< Wind speed (at scalar gridpoint) (m/s)
- mrt + t_zero, & !< Mean radiant temperature (°C)
- (hyp(k) + p(k,j,i)) * 0.01_wp, & !< Air pressure (hPa)
- pet ) !< output variable of PET
-
- local_pf(i,j,k) = pet
-
- ENDDO
- ELSE
- IF ( ALLOCATED( bio_pet_av ) ) THEN
- DO l = 1, nmrtbl
- i = mrtbl(ix,l)
- j = mrtbl(iy,l)
- k = mrtbl(iz,l)
- local_pf(i,j,k) = bio_pet_av(l)
- ENDDO
- ENDIF
- ENDIF
+ 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
@@ -411,453 +662,728 @@
END SELECT
-
- END SUBROUTINE biometeorology_data_output_3d
-
-
-
-!------------------------------------------------------------------------------!
-!
-! Description:
-! ------------
-!> PhysiologiCALLy Equivalent Temperature (PET),
-! stationary (calculated based on MEMI),
-! Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
-!
-! Input arguments:
-! ----------------
-! - ta : Air temperature (°C) REAL(wp)
-! - tmrt : Mean radiant temperature (°C) REAL(wp)
-! - v : Wind speed (m/s) REAL(wp)
-! - vpa : Vapor pressure (hPa) REAL(wp)
-! - p : Air pressure (hPa) REAL(wp)
-!
-! Output arguments:
-! ----------------
-! - tx : PET (°C) REAL(wp)
-!------------------------------------------------------------------------------!
-
- SUBROUTINE calculate_pet_static( &
-!-- Meteorological input
- ta, vpa, v, tmrt, p, &
+ 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
- tx ) !, &
-!-- Configure sample person (optional)
-! age, mbody, ht, work, eta, icl, fcl, pos, sex )
-
- IMPLICIT NONE
-
- REAL(wp), INTENT( IN ) :: ta, tmrt, v, vpa, p
-
- REAL(wp), INTENT ( OUT ) :: tx
-
-! REAL(wp), INTENT ( in ), optional :: age, mbody, ht, work, eta, icl, fcl
- REAL(wp) :: age, mbody, ht, work, eta, icl, fcl
-! INTEGER(iwp), INTENT ( in ), optional :: pos, sex
- INTEGER(iwp) :: pos, sex
-
- REAL(wp) :: acl, adu, aeff, cair, cb, emcl, emsk, ere, erel, esw, &
- evap, facl, feff, food, h, po, rdcl, rdsk, rob, rtv, &
- sigm, vpts, &
- ! former intent (out)
- ! - tsk : Skin temperature (°C) real
- ! - tcl : Clothing temperature (°C) real
- ! - ws : real
- ! - wetsk : Fraction of wet skin real
- tsk, tcl, wetsk
-
-
-!-- Person data
-! IF ( .NOT. present( age ) ) age = 35.
-! IF ( .NOT. present( mbody ) ) mbody = 75.
-! IF ( .NOT. present( ht ) ) ht = 1.75
-! IF ( .NOT. present( work ) ) work = 80.
-! IF ( .NOT. present( eta ) ) eta = 0.
-! IF ( .NOT. present( icl ) ) icl = 0.9
-! IF ( .NOT. present( fcl ) ) fcl = 1.15
-! IF ( .NOT. present( pos ) ) pos = 1
-! IF ( .NOT. present( sex ) ) sex = 1
-
- age = 35.
- mbody = 75.
- ht = 1.75
- work = 80.
- eta = 0.
- icl = 0.9
- fcl = 1.15
- pos = 1
- sex = 1
-
-!-- constants
- po = 1013.25 !< preassure at sea level
- ! p = 1013.25 !< local preassure (hPa), now defined as input variable
- rob = 1.06
- cb = 3.64 * 1000.
- food = 0.
- emsk = 0.99
- emcl = 0.95
- evap = 2.42 * 10. ** 6.
- sigm = 5.67 * 10. **(-8.)
-
-
-! write(9,*) 'Call calculate_pet_static(ta=', ta, ', vpa=', vpa, &
-! ', v=', v, ', tmrt=', tmrt, ', p=', p
-! flush(9)
-!-- call subfunctions
- CALL INKOERP ( age, cair, eta, ere, erel, evap, h, ht, mbody, &
- p, rtv, sex, ta, vpa, work )
-
-
- CALL BERECH ( acl, adu, aeff, cair, cb, emcl, emsk, &
- ere, erel, esw, evap, facl, fcl, feff, food, h, ht, icl, &
- mbody, p, po, rdcl, rdsk, rob, sex, sigm, &
- ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
-
-
- CALL PET ( acl, adu, aeff, cair, emcl, emsk, esw, evap, &
- facl, feff, h, p, po, rdcl, rdsk, &
- rtv, sigm, ta, tcl, tsk, tx, vpts, wetsk )
-
-
- END SUBROUTINE calculate_pet_static
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculate internal energy ballance
-!------------------------------------------------------------------------------!
- SUBROUTINE inkoerp( age, cair, eta, ere, erel, evap, h, ht, mbody, &
- & p, rtv, sex, ta, vpa, work )
-
-
- REAL(wp) :: age, cair, eta, ere, erel, eres, evap, h, ht, mbody, &
- & met, p, rtv, ta, tex, vpa, vpex, work
-
- INTEGER(iwp) :: sex
-
- IF ( sex .EQ. 1 ) THEN
- met = 3.45 * mbody ** ( 3. / 4. ) * (1. + 0.004 * &
- ( 30. - age) + 0.010 * ( ( ht * 100. / &
- ( mbody ** ( 1. / 3. ) ) ) - 43.4 ) )
- ELSE IF ( sex .EQ. 2 ) THEN
- met = 3.19 * mbody ** ( 3. / 4. ) * ( 1. + 0.004 * &
- ( 30. - age ) + 0.018 * ( ( ht * 100. / ( mbody ** &
- ( 1. / 3. ) ) ) - 42.1 ) )
- END IF
- met = work + met
-
- h = met * (1. - eta)
-
-!-- SENSIBLE RESPIRATION ENERGY
-
- cair = 1.01 * 1000.
- tex = 0.47 * ta + 21.0
- rtv = 1.44 * 10. ** (-6.) * met
- eres = cair * (ta - tex) * rtv
-
-!-- LATENT RESPIRATION ENERGY
-
- vpex = 6.11 * 10. ** ( 7.45 * tex / ( 235. + tex ) )
- erel = 0.623 * evap / p * ( vpa - vpex ) * rtv
-
-!-- SUM OF RESULTS
-
- ere = eres + erel
- RETURN
- END SUBROUTINE inkoerp
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculate heat gain or loss
-!------------------------------------------------------------------------------!
- SUBROUTINE BERECH( acl, adu, aeff, cair, cb, emcl, emsk, &
- ere, erel, esw, evap, facl, fcl, feff, food, h, ht, icl, &
- mbody, p, po, rdcl, rdsk, rob, sex, sigm, &
- ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
-
-
- REAL(wp) :: acl, adu, aeff, c(0:10), cair, cb, cbare, &
- cclo, csum, di, ed, emcl, emsk, enbal, &
- enbal2, ere, erel, esw, eswdif, eswphy, eswpot, &
- evap, facl, fcl, fec, feff, food, h, hc, he, ht, htcl, icl, &
- mbody, p, po, r1, r2, rbare, rcl, &
- rclo, rclo2, rdcl, rdsk, rob, rsum, sigm, sw, swf, swm, &
- ta, tbody, tcl, tcore(1:7), tmrt, tsk, v, vb, vb1, vb2, &
- vpa, vpts, wetsk, wd, wr, ws, wsum, xx, y
-
- INTEGER(iwp) :: count1, count3, j, sex
- logical :: skipIncreaseCount
-
- wetsk = 0.
- adu = 0.203 * mbody ** 0.425 * ht ** 0.725
-
- hc = 2.67 + ( 6.5 * v ** 0.67)
- hc = hc * (p /po) ** 0.55
- feff = 0.725 !< Posture: 0.725 for stading, 0.696 for sitting
- facl = (- 2.36 + 173.51 * icl - 100.76 * icl * icl + 19.28 &
- * (icl ** 3.)) / 100.
-
- IF ( facl .GT. 1. ) facl = 1.
- rcl = ( icl / 6.45) / facl
- IF ( icl .GE. 2. ) y = 1.
-
- IF ( ( icl .GT. 0.6 ) .AND. ( icl .LT. 2. ) ) y = ( ht - 0.2 ) / ht
- IF ( ( icl .LE. 0.6 ) .AND. ( icl .GT. 0.3 ) ) y = 0.5
- IF ( ( icl .LE. 0.3 ) .AND. ( icl .GT. 0. ) ) y = 0.1
-
- r2 = adu * (fcl - 1. + facl) / (2. * 3.14 * ht * y)
- r1 = facl * adu / (2. * 3.14 * ht * y)
-
- di = r2 - r1
-
-!-- SKIN TEMPERATURE
-
- DO j = 1,7
-
- tsk = 34.
- count1 = 0
- tcl = ( ta + tmrt + tsk ) / 3.
- count3 = 1
- enbal2 = 0.
-
- DO
- acl = adu * facl + adu * ( fcl - 1. )
- rclo2 = emcl * sigm * ( ( tcl + 273.2 )**4. - &
- ( tmrt + 273.2 )** 4. ) * feff
- htcl = 6.28 * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
- tsk = 1. / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
-
-!-- RADIATION SALDO
-
- aeff = adu * feff
- rbare = aeff * ( 1. - facl ) * emsk * sigm * &
- ( ( tmrt + 273.2 )** 4. - ( tsk + 273.2 )** 4. )
- rclo = feff * acl * emcl * sigm * &
- ( ( tmrt + 273.2 )** 4. - ( tcl + 273.2 )** 4. )
- rsum = rbare + rclo
-
-!-- CONVECTION
-
- cbare = hc * ( ta - tsk ) * adu * ( 1. - facl )
- cclo = hc * ( ta - tcl ) * acl
- csum = cbare + cclo
-
- !-- CORE TEMPERATUR
-
- c(0) = h + ere
- c(1) = adu * rob * cb
- c(2) = 18. - 0.5 * tsk
- c(3) = 5.28 * adu * c(2)
- c(4) = 0.0208 * c(1)
- c(5) = 0.76075 * c(1)
- c(6) = c(3) - c(5) - tsk * c(4)
- c(7) = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
- c(8) = c(6) * c(6) - 4. * c(4) * c(7)
- c(9) = 5.28 * adu - c(5) - c(4) * tsk
- c(10) = c(9) * c(9) - 4. * c(4) * &
- ( c(5) * tsk - c(0) - 5.28 * adu * tsk )
-
- IF ( tsk .EQ. 36. ) tsk = 36.01
- tcore(7) = c(0) / ( 5.28 * adu + c(1) * 6.3 / 3600. ) + tsk
- tcore(3) = c(0) / ( 5.28 * adu + ( c(1) * 6.3 / 3600. ) / &
- ( 1. + 0.5 * ( 34. - tsk ) ) ) + tsk
- IF ( c(10) .GE. 0.) THEN
- tcore(6) = ( - c(9) - c(10)**0.5 ) / ( 2. * c(4) )
- tcore(1) = ( - c(9) + c(10)**0.5 ) / ( 2. * c(4) )
- END IF
- ! 22
- IF ( c(8) .GE. 0. ) THEN
- tcore(2) = ( - c(6) + ABS( c(8) ) ** 0.5 ) / ( 2. * c(4) )
- tcore(5) = ( - c(6) - ABS( c(8) ) ** 0.5 ) / ( 2. * c(4) )
- tcore(4) = c(0) / ( 5.28 * adu + c(1) * 1. / 40. ) + tsk
- END IF
- ! 24
-
-!-- TRANSPIRATION
-
- tbody = 0.1 * tsk + 0.9 * tcore(j)
- swm = 304.94 * ( tbody - 36.6 ) * adu / 3600000.
- vpts = 6.11 * 10.**( 7.45 * tsk / ( 235. + tsk ) )
-
- IF ( tbody .LE. 36.6 ) swm = 0.
- ! swf = 0.7 * swm
-
- IF ( sex .EQ. 1 ) sw = swm
- IF ( sex .EQ. 2 ) sw = 0.7 * swm ! swf
- eswphy = - sw * evap
- he = 0.633 * hc / ( p * cair )
- fec = 1. / ( 1. + 0.92 * hc * rcl )
- eswpot = he * ( vpa - vpts ) * adu * evap * fec
- wetsk = eswphy / eswpot
-
- IF ( wetsk .GT. 1. ) wetsk = 1.
-
- eswdif = eswphy - eswpot
-
- IF ( eswdif .LE. 0. ) esw = eswpot
- IF ( eswdif .GT. 0. ) esw = eswphy
- IF ( esw .GT. 0. ) esw = 0.
-
-!-- DIFFUSION
-
- rdsk = 0.79 * 10. ** 7.
- rdcl = 0.
- ed = evap / ( rdsk + rdcl ) * adu * ( 1. - wetsk ) * ( vpa - vpts )
-
-!-- MAX VB
-
- vb1 = 34. - tsk
- vb2 = tcore(j) - 36.6
-
- IF ( vb2 .LT. 0. ) vb2 = 0.
- IF ( vb1 .LT. 0. ) vb1 = 0.
- vb = ( 6.3 + 75. * vb2 ) / ( 1. + 0.5 * vb1 )
-
-!-- ENERGY BALLANCE
-
- enbal = h + ed + ere + esw + csum + rsum + food
-
-
-!-- CLOTHING TEMPERATURE
-
- xx = 0.001
- IF ( count1 .EQ. 0 ) xx = 1.
- IF ( count1 .EQ. 1 ) xx = 0.1
- IF ( count1 .EQ. 2 ) xx = 0.01
- IF ( count1 .EQ. 3 ) xx = 0.001
-
- IF ( enbal .GT. 0. ) tcl = tcl + xx
- IF ( enbal .LT. 0. ) tcl = tcl - xx
-
- skipIncreaseCount = .FALSE.
- IF ( ( (enbal .LE. 0.) .AND. (enbal2 .GT. 0.) ) .OR. &
- ( ( enbal .GE. 0. ) .AND. ( enbal2 .LT. 0. ) ) ) THEN
- skipIncreaseCount = .TRUE.
- ELSE
- enbal2 = enbal
- count3 = count3 + 1
- END IF
-
- IF ( ( count3 .GT. 200 ) .OR. skipIncreaseCount ) THEN
- IF ( count1 .LT. 3 ) THEN
- count1 = count1 + 1
- enbal2 = 0.
- ELSE
- EXIT
- END IF
- END IF
- END DO
-
- IF ( count1 .EQ. 3 ) THEN
- SELECT CASE ( j )
- CASE ( 2, 5)
- IF ( .NOT. ( ( tcore(j) .GE. 36.6 ) .AND. &
- ( tsk .LE. 34.050 ) ) ) CYCLE
- CASE ( 6, 1 )
- IF ( c(10) .LT. 0. ) CYCLE
- IF ( .NOT. ( ( tcore(j) .GE. 36.6 ) .AND. &
- ( tsk .GT. 33.850 ) ) ) CYCLE
- CASE ( 3 )
- IF ( .NOT. ( ( tcore(j) .LT. 36.6 ) .AND. &
- ( tsk .LE. 34.000 ) ) ) CYCLE
- CASE ( 7 )
- IF ( .NOT. ( ( tcore(j) .LT. 36.6 ) .AND. &
- ( tsk .GT. 34.000 ) ) ) CYCLE
- CASE default ! = CASE ( 4 ), does actually nothing
- END SELECT
- END IF
-
- IF ( ( j .NE. 4 ) .AND. ( vb .GE. 91. ) ) CYCLE
- IF ( ( j .EQ. 4 ) .AND. ( vb .LT. 89. ) ) CYCLE
- IF ( vb .GT. 90.) vb = 90.
-
-!-- LOSSES BY WATER
-
- ws = sw * 3600. * 1000.
- IF ( ws .GT. 2000. ) ws = 2000.
- wd = ed / evap * 3600. * ( -1000. )
- wr = erel / evap * 3600. * ( -1000. )
-
- wsum = ws + wr + wd
-
- RETURN
- END DO
- RETURN
- END SUBROUTINE Berech
-
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Calculate PET
-!------------------------------------------------------------------------------!
- SUBROUTINE PET ( acl, adu, aeff, cair, emcl, emsk, esw, evap, &
- facl, feff, h, p, po, rdcl, rdsk, rtv, sigm, ta, tcl, tsk, tx, vpts, wetsk)
-
- REAL ( wp ) :: acl, adu, aeff, cair, cbare, cclo, csum, ed, &
- emcl, emsk, enbal, enbal2, ere, erel, eres, esw, evap, &
- facl, feff, h, hc, p, po, rbare, rclo, rdcl, rdsk, rsum, &
- rtv, sigm, ta, tcl, tex, tsk, tx, vpex, vpts, wetsk, xx
-
- INTEGER ( iwp ) :: count1
-
- tx = ta
- enbal2 = 0.
-
- DO count1 = 0, 3
- DO
- hc = 2.67 + 6.5 * 0.1 ** 0.67
- hc = hc * ( p / po ) ** 0.55
-
-!-- Radiation
-
- aeff = adu * feff
- rbare = aeff * ( 1. - facl ) * emsk * sigm * &
- ( ( tx + 273.2 ) ** 4. - ( tsk + 273.2 ) ** 4. )
- rclo = feff * acl * emcl * sigm * &
- ( ( tx + 273.2 ) ** 4. - ( tcl + 273.2 ) ** 4. )
- rsum = rbare + rclo
-
-!-- Covection
-
- cbare = hc * ( tx - tsk ) * adu * ( 1. - facl )
- cclo = hc * ( tx - tcl ) * acl
- csum = cbare + cclo
-
-!-- Diffusion
-
- ed = evap / ( rdsk + rdcl ) * adu * ( 1. - wetsk ) * ( 12. - vpts )
-
-!-- Respiration
-
- tex = 0.47 * tx + 21.
- eres = cair * ( tx - tex ) * rtv
- vpex = 6.11 * 10. ** ( 7.45 * tex / ( 235. + tex ) )
- erel = 0.623 * evap / p * ( 12. - vpex ) * rtv
- ere = eres + erel
-
-!-- Energy ballance
-
- enbal = h + ed + ere + esw + csum + rsum
-
-!-- Iteration concerning ta
-
- IF ( count1 .EQ. 0 ) xx = 1.
- IF ( count1 .EQ. 1 ) xx = 0.1
- IF ( count1 .EQ. 2 ) xx = 0.01
- IF ( count1 .EQ. 3 ) xx = 0.001
- IF ( enbal .GT. 0. ) tx = tx - xx
- IF ( enbal .LT. 0. ) tx = tx + xx
- IF ( ( enbal .LE. 0. ) .AND. ( enbal2 .GT. 0. ) ) EXIT
- IF ( ( enbal .GE. 0. ) .AND. ( enbal2 .LT. 0. ) ) EXIT
-
- enbal2 = enbal
- END DO
- END DO
- RETURN
- END SUBROUTINE PET
-
-END MODULE biometeorology_mod
+ 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
+
+ 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)
+
+ 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/ 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 - cels_offs
+
+ 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 + cels_offs )**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 - cels_offs
+
+ 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 ) ) - cels_offs
+
+ 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 ) ) - cels_offs
+
+ 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 + cels_offs, 333.15_wp ) )
+ vp = vp * pair / ( vp + 0.622_wp )
+ IF ( vp > vp_sat ) vp = vp_sat
+
+ tmrt = ta
+ IF ( radiation ) THEN
+ 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
+
+ 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 !< 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
+
+ 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