!> @file biometeorology_pt_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 .
!
! Copyright 1997-2019, Deutscher Wetterdienst (DWD) /
! German Meteorological Service (DWD) and
! Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
!
!
! Former revisions:
! -----------------
! $Id$
! From branch resler@3462, pavelkrc:
! replace specific precision DLOG, DMOD intrinsics with generic precision for
! compatibility with 'wp'
!
!
! Initial revision
!
!
!
! Authors:
! --------
! @author Dominik Froehlich
!
!
! Description:
! ------------
!> Module for the calculation of the thermal index perceived temperature (PT).
!> The Perceived Temperature is the air temperature in a reference environ-
!> ment that results in the same Predicted Mean Vote (Fanger 1972), PMV, as
!> the real environment. In the reference environment the mean radiant
!> temperature is equal to the air temperature and the wind speed is reduced
!> to a value set at 0.1 m/s. The Perceived Temperature is linked to a
!> reference individual: male aged 35 years, body height 1.75 m, weight
!> 75 kg, walking with 4 km per hour on a horizontal plain (related to an
!> internal heat production of 172.5 W). The person varies the thermal
!> resistance of clothing (clo) in the range Icl=1.75 clo (winter) and
!> Icl=0.50 clo (summer) to achieve thermal comfort (PMV = 0) if possible.
!> If not possible the individual perceives cold stress (winter, 1.75 clo)
!> or heat load (summer, 0.5 clo).
!>
!> @todo
!>
!> @note nothing now
!>
!> @bug no known bugs by now
!------------------------------------------------------------------------------!
MODULE biometeorology_pt_mod
!-- Load required variables from existing modules
USE kinds !< to set precision of INTEGER and REAL arrays according to PALM
IMPLICIT NONE
PRIVATE !-- no private interfaces --!
PUBLIC calculate_pt_static, saturation_vapor_pressure, deltapmv, dpmv_adj, &
dpmv_cold, fanger, calc_sultr, pt_regression, ireq_neutral, iso_ridder
!-- PALM interfaces:
INTERFACE calculate_pt_static
MODULE PROCEDURE calculate_pt_static
END INTERFACE calculate_pt_static
INTERFACE saturation_vapor_pressure
MODULE PROCEDURE saturation_vapor_pressure
END INTERFACE saturation_vapor_pressure
INTERFACE deltapmv
MODULE PROCEDURE deltapmv
END INTERFACE deltapmv
INTERFACE dpmv_adj
MODULE PROCEDURE dpmv_adj
END INTERFACE dpmv_adj
INTERFACE dpmv_cold
MODULE PROCEDURE dpmv_cold
END INTERFACE dpmv_cold
INTERFACE fanger
MODULE PROCEDURE fanger
END INTERFACE fanger
INTERFACE calc_sultr
MODULE PROCEDURE calc_sultr
END INTERFACE calc_sultr
INTERFACE pt_regression
MODULE PROCEDURE pt_regression
END INTERFACE pt_regression
INTERFACE ireq_neutral
MODULE PROCEDURE ireq_neutral
END INTERFACE ireq_neutral
INTERFACE iso_ridder
MODULE PROCEDURE iso_ridder
END INTERFACE iso_ridder
CONTAINS
!------------------------------------------------------------------------------!
! Description:
! ------------
!> PT_BASIC.F90 Version of perceived temperature (PT, �C) for
!> - standard measured/predicted meteorological values and TMRT
!> as input;
!> - regressions for determination of PT;
!> - adjustment to Gagge's PMV* (2-node-model, 1986) as base of PT
!> under warm/humid conditions (Icl= 0.50 clo) and under cold
!> conditions (Icl= 1.75 clo)
!>
!> Function value is the Perceived Temperature, degree centigrade
!------------------------------------------------------------------------------!
SUBROUTINE calculate_pt_static( tt2m, el2m, vau1m, tmrt, pb, clo, pt_basic )
IMPLICIT NONE
!-- Type of input of the argument list
REAL(wp), INTENT ( IN ) :: tt2m !< Local air temperature (�C)
REAL(wp), INTENT ( IN ) :: el2m !< Local vapour pressure (hPa)
REAL(wp), INTENT ( IN ) :: tmrt !< Local mean radiant temperature (�C)
REAL(wp), INTENT ( IN ) :: vau1m !< Local wind velocitry (m/s)
REAL(wp), INTENT ( IN ) :: pb !< Local barometric air pressure (hPa)
!-- Type of output of the argument list
REAL(wp), INTENT ( OUT ) :: pt_basic !< Perceived temperature (�C)
REAL(wp), INTENT ( OUT ) :: clo !< Clothing index (dimensionless)
!-- Parameters for standard "Klima-Michel"
REAL(wp), PARAMETER :: eta = 0._wp !< Mechanical work efficiency for walking on flat ground (compare to Fanger (1972) pp 24f)
REAL(wp), PARAMETER :: actlev = 134.6862_wp !< Workload by activity per standardized surface (A_Du)
!-- Type of program variables
REAL(wp), PARAMETER :: eps = 0.0005 !< Accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
REAL(wp) :: uclo
REAL(wp) :: oclo
REAL(wp) :: d_pmv
REAL(wp) :: svp_tt
REAL(wp) :: sult_lim
REAL(wp) :: dgtcm
REAL(wp) :: dgtcstd
REAL(wp) :: clon
REAL(wp) :: ireq_minimal
REAL(wp) :: clo_fanger
REAL(wp) :: pmv_o
REAL(wp) :: pmv_u
REAL(wp) :: pmva
REAL(wp) :: ptc
REAL(wp) :: d_std
REAL(wp) :: pmvs
REAL(wp) :: top
INTEGER(iwp) :: nzaehl, nerr_kalt, nerr
LOGICAL :: sultrieness
!-- Initialise
pt_basic = 9999.0_wp
nerr = 0_iwp
nzaehl = 0_iwp
sultrieness = .FALSE.
!-- Tresholds: clothing insulation (account for model inaccuracies)
! summer clothing
uclo = 0.44453_wp
! winter clothing
oclo = 1.76267_wp
!-- decision: firstly calculate for winter or summer clothing
IF ( tt2m <= 10._wp ) THEN
!-- First guess: winter clothing insulation: cold stress
clo = oclo
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, clo, actlev, eta, pmva, top )
pmv_o = pmva
IF ( pmva > 0._wp ) THEN
!-- Case summer clothing insulation: heat load ?
clo = uclo
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, clo, actlev, eta, pmva, &
top )
pmv_u = pmva
IF ( pmva <= 0._wp ) THEN
!-- Case: comfort achievable by varying clothing insulation
!-- Between winter and summer set values
CALL iso_ridder ( tt2m, tmrt, el2m, vau1m, pb, actlev, eta, uclo, &
pmv_u, oclo, pmv_o, eps, pmva, top, nzaehl, clo )
IF ( nzaehl < 0_iwp ) THEN
nerr = -1_iwp
RETURN
ENDIF
ELSE IF ( pmva > 0.06_wp ) THEN
clo = 0.5_wp
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, clo, actlev, eta, &
pmva, top )
ENDIF
ELSE IF ( pmva < -0.11_wp ) THEN
clo = 1.75_wp
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, clo, actlev, eta, pmva, &
top )
ENDIF
ELSE
!-- First guess: summer clothing insulation: heat load
clo = uclo
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, clo, actlev, eta, pmva, top )
pmv_u = pmva
IF ( pmva < 0._wp ) THEN
!-- Case winter clothing insulation: cold stress ?
clo = oclo
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, clo, actlev, eta, pmva, &
top )
pmv_o = pmva
IF ( pmva >= 0._wp ) THEN
!-- Case: comfort achievable by varying clothing insulation
! between winter and summer set values
CALL iso_ridder ( tt2m, tmrt, el2m, vau1m, pb, actlev, eta, uclo, &
pmv_u, oclo, pmv_o, eps, pmva, top, nzaehl, clo )
IF ( nzaehl < 0_iwp ) THEN
nerr = -1_iwp
RETURN
ENDIF
ELSE IF ( pmva < -0.11_wp ) THEN
clo = 1.75_wp
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, clo, actlev, eta, &
pmva, top )
ENDIF
ELSE IF ( pmva > 0.06_wp ) THEN
clo = 0.5_wp
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, clo, actlev, eta, pmva, &
top )
ENDIF
ENDIF
!-- Determine perceived temperature by regression equation + adjustments
pmvs = pmva
CALL pt_regression ( pmva, clo, pt_basic )
ptc = pt_basic
IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
!-- Adjust for cold conditions according to Gagge 1986
CALL dpmv_cold ( pmva, tt2m, vau1m, tmrt, nerr_kalt, d_pmv )
IF ( nerr_kalt > 0_iwp ) nerr = -5_iwp
pmvs = pmva - d_pmv
IF ( pmvs > -0.11_wp ) THEN
d_pmv = 0._wp
pmvs = -0.11_wp
ENDIF
CALL pt_regression ( pmvs, clo, pt_basic )
ENDIF
clo_fanger = clo
clon = clo
IF ( clo > 0.5_wp .AND. pt_basic <= 8.73_wp ) THEN
!-- Required clothing insulation (ireq) is exclusively defined for
! operative temperatures (top) less 10 (C) for a
! reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
clon = ireq_neutral ( pt_basic, ireq_minimal, nerr )
clo = clon
ENDIF
CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
sultrieness = .FALSE.
d_std = -99._wp
IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
!-- Adjust for warm/humid conditions according to Gagge 1986
CALL saturation_vapor_pressure ( tt2m, svp_tt )
d_pmv = deltapmv ( pmva, tt2m, el2m, svp_tt, tmrt, vau1m, nerr )
pmvs = pmva + d_pmv
CALL pt_regression ( pmvs, clo, pt_basic )
IF ( sult_lim < 99._wp ) THEN
IF ( (pt_basic - ptc) > sult_lim ) sultrieness = .TRUE.
!-- Set factor to threshold for sultriness
IF ( dgtcstd /= 0_iwp ) THEN
d_std = ( ( pt_basic - ptc ) - dgtcm ) / dgtcstd
ENDIF
ENDIF
ENDIF
END SUBROUTINE calculate_pt_static
!------------------------------------------------------------------------------!
! Description:
! ------------
!> The SUBROUTINE calculates the saturation water vapour pressure
!> (hPa = hecto Pascal) for a given temperature tt (�C).
!> For example, tt can be the air temperature or the dew point temperature.
!------------------------------------------------------------------------------!
SUBROUTINE saturation_vapor_pressure( tt, p_st )
IMPLICIT NONE
REAL(wp), INTENT ( IN ) :: tt !< ambient air temperature (�C)
REAL(wp), INTENT ( OUT ) :: p_st !< saturation water vapour pressure (hPa)
REAL(wp) :: b
REAL(wp) :: c
IF ( tt < 0._wp ) THEN
!-- tt < 0 (�C): saturation water vapour pressure over ice
b = 17.84362_wp
c = 245.425_wp
ELSE
!-- tt >= 0 (�C): saturation water vapour pressure over water
b = 17.08085_wp
c = 234.175_wp
ENDIF
!-- Saturation water vapour pressure
p_st = 6.1078_wp * EXP ( b * tt / ( c + tt ) )
END SUBROUTINE saturation_vapor_pressure
!------------------------------------------------------------------------------!
! Description:
! ------------
! Find the clothing insulation value clo_res (clo) to make Fanger's Predicted
! Mean Vote (PMV) equal comfort (pmva=0) for actual meteorological conditions
! (tt2m,tmrt, el2m, vau1m, pb) and values of individual's activity level
! (Klima-Michel: 134.682 W/m2 = 2.3 met, mechanincal work load eta= 0 W/m2),
! i.e. find the root of Fanger's comfort equation
!------------------------------------------------------------------------------!
SUBROUTINE iso_ridder( tt2m, tmrt, el2m, vau1m, pb, actlev, eta, uclo, &
pmv_u, oclo, pmv_o, eps, pmva, top, nerr, &
clo_res )
IMPLICIT NONE
!-- Input variables of argument list:
REAL(wp), INTENT ( IN ) :: tt2m !< Ambient temperature (�C)
REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (�C)
REAL(wp), INTENT ( IN ) :: el2m !< Water vapour pressure (hPa)
REAL(wp), INTENT ( IN ) :: vau1m !< Wind speed (m/s) 1 m above ground
REAL(wp), INTENT ( IN ) :: pb !< Barometric pressure (hPa)
REAL(wp), INTENT ( IN ) :: actlev !< Individuals activity level per unit surface area (W/m2)
REAL(wp), INTENT ( IN ) :: eta !< Individuals work efficiency (dimensionless)
REAL(wp), INTENT ( IN ) :: uclo !< Lower threshold of bracketing clothing insulation (clo)
REAL(wp), INTENT ( IN ) :: oclo !< Upper threshold of bracketing clothing insulation (clo)
REAL(wp), INTENT ( IN ) :: eps !< (0.05) accuracy in clothing insulation (clo) for
! evaluation the root of Fanger's PMV (pmva=0)
REAL(wp), INTENT ( IN ) :: pmv_o !< Fanger's PMV corresponding to oclo
REAL(wp), INTENT ( IN ) :: pmv_u !< Fanger's PMV corresponding to uclo
! Output variables of argument list:
REAL(wp), INTENT ( OUT ) :: pmva !< 0 (set to zero, because clo is evaluated for comfort)
REAL(wp), INTENT ( OUT ) :: top !< Operative temperature (�C) at found root of Fanger's PMV
REAL(wp), INTENT ( OUT ) :: clo_res !< Resulting clothing insulation value (clo)
INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
! nerr >= 0, o.k., and nerr is the number of iterations for
! convergence
! nerr = -1: error = malfunction of Ridder's convergence method
! nerr = -2: error = maximum iterations (max_iteration) exceeded
! nerr = -3: error = root not bracketed between uclo and oclo
!-- Type of program variables
INTEGER(iwp), PARAMETER :: max_iteration = 15_iwp
REAL(wp), PARAMETER :: guess_0 = -1.11e30_wp
REAL(wp) :: x_ridder
REAL(wp) :: clo_u
REAL(wp) :: clo_o
REAL(wp) :: x_unten
REAL(wp) :: x_oben
REAL(wp) :: x_mittel
REAL(wp) :: x_neu
REAL(wp) :: y_unten
REAL(wp) :: y_oben
REAL(wp) :: y_mittel
REAL(wp) :: y_neu
REAL(wp) :: s
INTEGER(iwp) :: j
!-- Initialise
nerr = 0_iwp
!-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation
pmva = 0._wp
clo_u = uclo
y_unten = pmv_u
clo_o = oclo
y_oben = pmv_o
IF ( ( y_unten > 0._wp .AND. y_oben < 0._wp ) .OR. &
( y_unten < 0._wp .AND. y_oben > 0._wp ) ) THEN
x_unten = clo_u
x_oben = clo_o
x_ridder = guess_0
DO j = 1_iwp, max_iteration
x_mittel = 0.5_wp * ( x_unten + x_oben )
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, x_mittel, actlev, eta, &
y_mittel, top )
s = SQRT ( y_mittel**2 - y_unten * y_oben )
IF ( s == 0._wp ) THEN
clo_res = x_mittel
nerr = j
RETURN
ENDIF
x_neu = x_mittel + ( x_mittel - x_unten ) * &
( SIGN ( 1._wp, y_unten - y_oben ) * y_mittel / s )
IF ( ABS ( x_neu - x_ridder ) <= eps ) THEN
clo_res = x_ridder
nerr = j
RETURN
ENDIF
x_ridder = x_neu
CALL fanger ( tt2m, tmrt, el2m, vau1m, pb, x_ridder, actlev, eta, &
y_neu, top )
IF ( y_neu == 0._wp ) THEN
clo_res = x_ridder
nerr = j
RETURN
ENDIF
IF ( SIGN ( y_mittel, y_neu ) /= y_mittel ) THEN
x_unten = x_mittel
y_unten = y_mittel
x_oben = x_ridder
y_oben = y_neu
ELSE IF ( SIGN ( y_unten, y_neu ) /= y_unten ) THEN
x_oben = x_ridder
y_oben = y_neu
ELSE IF ( SIGN ( y_oben, y_neu ) /= y_oben ) THEN
x_unten = x_ridder
y_unten = y_neu
ELSE
!-- Never get here in x_ridder: singularity in y
nerr = -1_iwp
clo_res = x_ridder
RETURN
ENDIF
IF ( ABS ( x_oben - x_unten ) <= eps ) THEN
clo_res = x_ridder
nerr = j
RETURN
ENDIF
ENDDO
!-- x_ridder exceed maximum iterations
nerr = -2_iwp
clo_res = y_neu
RETURN
ELSE IF ( y_unten == 0. ) THEN
x_ridder = clo_u
ELSE IF ( y_oben == 0. ) THEN
x_ridder = clo_o
ELSE
!-- x_ridder not bracketed by u_clo and o_clo
nerr = -3_iwp
clo_res = x_ridder
RETURN
ENDIF
END SUBROUTINE iso_ridder
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Regression relations between perceived temperature (gt) and (adjusted)
!> PMV (21.11.1996, files F3PMVPGT.TXT und F3CLOGT.TXT). Three cases:
!>
!> - PMV < 0: gt = 5.805 + 12.6784*pmv
!> - PMV = 0 (gt depends on clothing insulation, FUNCTION iso_ridder):
!> gt = 21.258 - 9.558*clo
!> - PMV > 0: 16.826 + 6.163*pmv
!>
!> The regression presumes the Klima-Michel settings for reference
!> individual and reference environment.
!------------------------------------------------------------------------------!
SUBROUTINE pt_regression( pmv, clo, pt )
IMPLICIT NONE
REAL(wp), INTENT ( IN ) :: pmv !< Fangers predicted mean vote (dimensionless)
REAL(wp), INTENT ( IN ) :: clo !< clothing insulation index (clo)
REAL(wp), INTENT ( OUT ) :: pt !< pt (�C) corresponding to given PMV / clo
IF ( pmv <= -0.11_wp ) THEN
pt = 5.805_wp + 12.6784_wp * pmv
ELSE
IF ( pmv >= + 0.01_wp ) THEN
pt = 16.826_wp + 6.163_wp * pmv
ELSE
pt = 21.258_wp - 9.558_wp * clo
ENDIF
ENDIF
END SUBROUTINE pt_regression
!------------------------------------------------------------------------------!
! Description:
! ------------
!> FANGER.F90
!>
!> SI-VERSION: ACTLEV W m-2, DAMPFDRUCK hPa
!> Berechnet das aktuelle Predicted Mean Vote nach Fanger
!>
!> The case of free convection (vau < 0.1 m/s) is dealt with vau = 0.1 m/s
!------------------------------------------------------------------------------!
SUBROUTINE fanger( tt, tmrt, pa, in_vau, pb, in_clo, actlev, eta, pmva, top )
IMPLICIT NONE
!-- Input variables of argument list:
REAL(wp), INTENT ( IN ) :: tt !< Ambient air temperature (�C)
REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (�C)
REAL(wp), INTENT ( IN ) :: pa !< Water vapour pressure (hPa)
REAL(wp), INTENT ( IN ) :: pb !< Barometric pressure (hPa) at site
REAL(wp), INTENT ( IN ) :: in_vau !< Wind speed (m/s) 1 m above ground
REAL(wp), INTENT ( IN ) :: in_clo !< Clothing insulation (clo)
REAL(wp), INTENT ( IN ) :: actlev !< Individuals activity level per unit surface area (W/m2)
REAL(wp), INTENT ( IN ) :: eta !< Individuals mechanical work efficiency (dimensionless)
!-- Output variables of argument list:
REAL(wp), INTENT ( OUT ) :: pmva !< Actual Predicted Mean Vote (dimensionless) according
! to Fanger corresponding to meteorological (tt,tmrt,pa,vau,pb)
! and individual variables (clo, actlev, eta)
REAL(wp), INTENT ( OUT ) :: top !< operative temperature (�C)
!-- Internal variables
REAL(wp), PARAMETER :: cels_offs = 273.15_wp !< Kelvin Celsius Offset (K)
REAL(wp) :: f_cl
REAL(wp) :: heat_convection
REAL(wp) :: activity
REAL(wp) :: t_skin_aver
REAL(wp) :: bc
REAL(wp) :: cc
REAL(wp) :: dc
REAL(wp) :: ec
REAL(wp) :: gc
REAL(wp) :: t_clothing
REAL(wp) :: hr
REAL(wp) :: clo
REAL(wp) :: vau
REAL(wp) :: z1
REAL(wp) :: z2
REAL(wp) :: z3
REAL(wp) :: z4
REAL(wp) :: z5
REAL(wp) :: z6
INTEGER(iwp) :: i
!-- Clo must be > 0. to avoid div. by 0!
clo = in_clo
IF ( clo <= 0._wp ) clo = .001_wp
!-- f_cl = Increase in surface due to clothing
f_cl = 1._wp + .15_wp * clo
!-- Case of free convection (vau < 0.1 m/s ) not considered
vau = in_vau
IF ( vau < .1_wp ) THEN
vau = .1_wp
ENDIF
!-- Heat_convection = forced convection
heat_convection = 12.1_wp * SQRT ( vau * pb / 1013.25_wp )
!-- Activity = inner heat produktion per standardized surface
activity = actlev * ( 1._wp - eta )
!-- T_skin_aver = average skin temperature
t_skin_aver = 35.7_wp - .0275_wp * activity
!-- Calculation of constants for evaluation below
bc = .155_wp * clo * 3.96_wp * 10._wp**( -8 ) * f_cl
cc = f_cl * heat_convection
ec = .155_wp * clo
dc = ( 1._wp + ec * cc ) / bc
gc = ( t_skin_aver + bc * ( tmrt + cels_offs )**4 + ec * cc * tt ) / bc
!-- Calculation of clothing surface temperature (t_clothing) based on
! Newton-approximation with air temperature as initial guess
t_clothing = tt
DO I = 1, 3
t_clothing = t_clothing - ( ( t_clothing + cels_offs )**4 + t_clothing &
* dc - gc ) / ( 4._wp * ( t_clothing + cels_offs )**3 + dc )
ENDDO
!-- Empiric factor for the adaption of the heat ballance equation
! to the psycho-physical scale (Equ. 40 in FANGER)
z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
!-- Water vapour diffution through the skin
z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
!-- Sweat evaporation from the skin surface
z3 = .42_wp * ( activity - 58._wp )
!-- Loss of latent heat through respiration
z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev * &
( 34._wp - tt )
!-- Loss of radiational heat
z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + cels_offs )**4 - ( tmrt + &
cels_offs )**4 )
IF ( ABS ( t_clothing - tmrt ) > 0._wp ) THEN
hr = z5 / f_cl / ( t_clothing - tmrt )
ELSE
hr = 0._wp
ENDIF
!-- Heat loss through forced convection cc*(t_clothing-TT)
z6 = cc * ( t_clothing - tt )
!-- Predicted Mean Vote
pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
!-- Operative temperatur
top = ( hr * tmrt + heat_convection * tt ) / ( hr + heat_convection )
END SUBROUTINE fanger
!------------------------------------------------------------------------------!
! Description:
! ------------
!> For pmva > 0 and clo =0.5 the increment (deltapmv) is calculated
!> that converts pmva into Gagge's et al. (1986) PMV*.
!------------------------------------------------------------------------------!
REAL(wp) FUNCTION deltapmv( pmva, tt2m, el2m, svp_tt, tmrt, vau1m, nerr )
IMPLICIT NONE
!-- Input variables of argument list:
REAL(wp), INTENT ( IN ) :: pmva !< Actual Predicted Mean Vote (PMV) according to Fanger
REAL(wp), INTENT ( IN ) :: tt2m !< Ambient temperature (�C) at screen level
REAL(wp), INTENT ( IN ) :: el2m !< Water vapour pressure (hPa) at screen level
REAL(wp), INTENT ( IN ) :: svp_tt !< Saturation water vapour pressure (hPa) at tt2m
REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (�C) at screen level
REAL(wp), INTENT ( IN ) :: vau1m !< Wind speed (m/s) 1 m above ground
!-- Output variables of argument list:
INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
! 0 = o.k.
! -2 = pmva outside valid regression range
! -3 = rel. humidity set to 5 % or 95 %, respectively
! -4 = deltapmv set to avoid pmvs < 0
!-- Internal variable types:
REAL(wp) :: pmv
REAL(wp) :: pa_p50
REAL(wp) :: pa
REAL(wp) :: apa
REAL(wp) :: dapa
REAL(wp) :: sqvel
REAL(wp) :: ta
REAL(wp) :: dtmrt
REAL(wp) :: p10
REAL(wp) :: p95
REAL(wp) :: gew
REAL(wp) :: gew2
REAL(wp) :: dpmv_1
REAL(wp) :: dpmv_2
REAL(wp) :: pmvs
REAL(wp) :: bpmv(0:7)
REAL(wp) :: bpa_p50(0:7)
REAL(wp) :: bpa(0:7)
REAL(wp) :: bapa(0:7)
REAL(wp) :: bdapa(0:7)
REAL(wp) :: bsqvel(0:7)
REAL(wp) :: bta(0:7)
REAL(wp) :: bdtmrt(0:7)
REAL(wp) :: aconst(0:7)
INTEGER(iwp) :: nreg
DATA bpmv / &
-0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp, -0.3551048_wp,&
-0.4304076_wp, -0.4884961_wp, -0.4897495_wp /
DATA bpa_p50 / &
-0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp, +0.4285938_wp,&
+0.6281256_wp, +0.5067361_wp, +0.3965169_wp /
DATA bpa / &
+0.0580284_wp, +0.0836264_wp, +0.1009919_wp, +0.1020777_wp, +0.0898681_wp,&
+0.0839116_wp, +0.0853258_wp, +0.0866589_wp /
DATA bapa / &
-1.7838788_wp, -2.9306231_wp, -1.6350334_wp, +0.6211547_wp, +3.3918083_wp,&
+5.5521025_wp, +8.4897418_wp, +16.6265851_wp /
DATA bdapa / &
+1.6752720_wp, +2.7379504_wp, +1.2940526_wp, -1.0985759_wp, -3.9054732_wp,&
-6.0403012_wp, -8.9437119_wp, -17.0671201_wp /
DATA bsqvel / &
-0.0315598_wp, -0.0286272_wp, -0.0009228_wp, +0.0483344_wp, +0.0992366_wp,&
+0.1491379_wp, +0.1951452_wp, +0.2133949_wp /
DATA bta / &
+0.0953986_wp, +0.1524760_wp, +0.0564241_wp, -0.0893253_wp, -0.2398868_wp,&
-0.3515237_wp, -0.5095144_wp, -0.9469258_wp /
DATA bdtmrt / &
-0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp, -0.0069036_wp,&
-0.0075844_wp, -0.0079602_wp, -0.0089439_wp /
DATA aconst / &
+1.8686215_wp, +3.4260713_wp, +2.0116185_wp, -0.7777552_wp, -4.6715853_wp,&
-7.7314281_wp, -11.7602578_wp, -23.5934198_wp /
!-- Test for compliance with regression range
IF ( pmva < -1.0_wp .OR. pmva > 7.0_wp ) THEN
nerr = -2_iwp
ELSE
nerr = 0_iwp
ENDIF
!-- Initialise classic PMV
pmv = pmva
!-- Water vapour pressure of air
p10 = 0.05_wp * svp_tt
p95 = 1.00_wp * svp_tt
IF ( el2m >= p10 .AND. el2m <= p95 ) THEN
pa = el2m
ELSE
nerr = -3_iwp
IF ( el2m < p10 ) THEN
!-- Due to conditions of regress�on: r.H. >= 5 %
pa = p10
ELSE
!-- Due to conditions of regress�on: r.H. <= 95 %
pa = p95
ENDIF
ENDIF
IF ( pa > 0._wp ) THEN
!-- Natural logarithm of pa
apa = LOG ( pa )
ELSE
apa = -5._wp
ENDIF
!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
pa_p50 = 0.5_wp * svp_tt
IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN
dapa = apa - LOG ( pa_p50 )
pa_p50 = pa / pa_p50
ELSE
dapa = -5._wp
pa_p50 = 0._wp
ENDIF
!-- Square root of wind velocity
IF ( vau1m >= 0._wp ) THEN
sqvel = SQRT ( vau1m )
ELSE
sqvel = 0._wp
ENDIF
!-- Air temperature
ta = tt2m
!-- Difference mean radiation to air temperature
dtmrt = tmrt - ta
!-- Select the valid regression coefficients
nreg = INT ( pmv )
IF ( nreg < 0_iwp ) THEN
!-- value of the FUNCTION in the case pmv <= -1
deltapmv = 0._wp
RETURN
ENDIF
gew = MOD ( pmv, 1._wp )
IF ( gew < 0._wp ) gew = 0._wp
IF ( nreg > 5_iwp ) THEN
! nreg=6
nreg = 5_iwp
gew = pmv - 5._wp
gew2 = pmv - 6._wp
IF ( gew2 > 0_iwp ) THEN
gew = ( gew - gew2 ) / gew
ENDIF
ENDIF
!-- Regression valid for 0. <= pmv <= 6.
dpmv_1 = &
+ bpa ( nreg ) * pa &
+ bpmv ( nreg ) * pmv &
+ bapa ( nreg ) * apa &
+ bta ( nreg ) * ta &
+ bdtmrt ( nreg ) * dtmrt &
+ bdapa ( nreg ) * dapa &
+ bsqvel ( nreg ) * sqvel &
+ bpa_p50 ( nreg ) * pa_p50 &
+ aconst ( nreg )
dpmv_2 = 0._wp
IF ( nreg < 6_iwp ) THEN
dpmv_2 = &
+ bpa ( nreg + 1_iwp ) * pa &
+ bpmv ( nreg + 1_iwp ) * pmv &
+ bapa ( nreg + 1_iwp ) * apa &
+ bta ( nreg + 1_iwp ) * ta &
+ bdtmrt ( nreg + 1_iwp ) * dtmrt &
+ bdapa ( nreg + 1_iwp ) * dapa &
+ bsqvel ( nreg + 1_iwp ) * sqvel &
+ bpa_p50 ( nreg + 1_iwp ) * pa_p50 &
+ aconst ( nreg + 1_iwp )
ENDIF
!-- Calculate pmv modification
deltapmv = ( 1._wp - gew ) * dpmv_1 + gew * dpmv_2
pmvs = pmva + deltapmv
IF ( ( pmvs ) < 0._wp ) THEN
!-- Prevent negative pmv* due to problems with clothing insulation
nerr = -4_iwp
IF ( pmvs > -0.11_wp ) THEN
!-- Threshold from FUNCTION pt_regression for winter clothing insulation
deltapmv = deltapmv + 0.11_wp
ELSE
!-- Set pmvs to "0" for compliance with summer clothing insulation
deltapmv = -1._wp * pmva
ENDIF
ENDIF
END FUNCTION deltapmv
!------------------------------------------------------------------------------!
! Description:
! ------------
!> The subroutine "calc_sultr" returns a threshold value to perceived
!> temperature allowing to decide whether the actual perceived temperature
!> is linked to perecption of sultriness. The threshold values depends
!> on the Fanger's classical PMV, expressed here as perceived temperature
!> gtc. It is derived from all synoptic observations 1966-1995 of the sites
!> 10170 Rostock-Warnemuende and 10803 Freiburg showing heat load. The
!> threshold value is valid for subjects acclimatised to Central European
!> climate conditions. It includes all effective influences as air
!> temperature, relative humidity, mean radiant temperature, and wind
!> velocity:
!> - defined only under warm conditions: gtc based on 0.5 (clo) and
!> gtc > 16.826 (�C)
!> - undefined: sultr_res set at +99.
!------------------------------------------------------------------------------!
SUBROUTINE calc_sultr( gtc, dgtcm, dgtcstd, sultr_res )
IMPLICIT NONE
!-- Input of the argument list:
REAL(wp), INTENT ( IN ) :: gtc !< Classical perceived temperature: Base is Fanger's PMV
!-- Additional output variables of argument list:
REAL(wp), INTENT ( OUT ) :: dgtcm !< Mean deviation gtc (classical gt) to gt* (rational gt
! calculated based on Gagge's rational PMV*)
REAL(wp), INTENT ( OUT ) :: dgtcstd !< dgtcm plus its standard deviation times a factor
! determining the significance to perceive sultriness
REAL(wp), INTENT ( OUT ) :: sultr_res
!-- Types of coefficients mean deviation: third order polynomial
REAL(wp), PARAMETER :: dgtcka = +7.5776086_wp
REAL(wp), PARAMETER :: dgtckb = -0.740603_wp
REAL(wp), PARAMETER :: dgtckc = +0.0213324_wp
REAL(wp), PARAMETER :: dgtckd = -0.00027797237_wp
!-- Types of coefficients mean deviation plus standard deviation
! regression coefficients: third order polynomial
REAL(wp), PARAMETER :: dgtcsa = +0.0268918_wp
REAL(wp), PARAMETER :: dgtcsb = +0.0465957_wp
REAL(wp), PARAMETER :: dgtcsc = -0.00054709752_wp
REAL(wp), PARAMETER :: dgtcsd = +0.0000063714823_wp
!-- Factor to mean standard deviation defining SIGNificance for
! sultriness
REAL(wp), PARAMETER :: faktor = 1._wp
!-- Initialise
sultr_res = +99._wp
dgtcm = 0._wp
dgtcstd = 999999._wp
IF ( gtc < 16.826_wp .OR. gtc > 56._wp ) THEN
!-- Unallowed classical PMV/gtc
RETURN
ENDIF
!-- Mean deviation dependent on gtc
dgtcm = dgtcka + dgtckb * gtc + dgtckc * gtc**2._wp + dgtckd * gtc**3._wp
!-- Mean deviation plus its standard deviation
dgtcstd = dgtcsa + dgtcsb * gtc + dgtcsc * gtc**2._wp + dgtcsd * gtc**3._wp
!-- Value of the FUNCTION
sultr_res = dgtcm + faktor * dgtcstd
IF ( ABS ( sultr_res ) > 99._wp ) sultr_res = +99._wp
END SUBROUTINE calc_sultr
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Multiple linear regression to calculate an increment delta_kalt,
!> to adjust Fanger's classical PMV (pmva) by Gagge's 2 node model,
!> applying Fanger's convective heat transfer coefficient, hcf.
!> Wind velocitiy of the reference environment is 0.10 m/s
!------------------------------------------------------------------------------!
SUBROUTINE dpmv_cold( pmva, tt2m, vau1m, tmrt, nerr, dpmv_cold_res )
IMPLICIT NONE
!-- Type of input arguments
REAL(wp), INTENT ( IN ) :: pmva !< Fanger's classical predicted mean vote
REAL(wp), INTENT ( IN ) :: tt2m !< Air temperature (�C) 2 m above ground
REAL(wp), INTENT ( IN ) :: vau1m !< Relative wind velocity 1 m above ground (m/s)
REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (�C)
!-- Type of output argument
INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error indicator: 0 = o.k., +1 = denominator for intersection = 0
REAL(wp), INTENT ( OUT ) :: dpmv_cold_res !< Increment to adjust pmva according to the results of Gagge's
! 2 node model depending on the input
!-- Type of program variables
REAL(wp) :: delta_kalt(3)
REAL(wp) :: pmv_cross(2)
REAL(wp) :: reg_a(3)
REAL(wp) :: coeff(3,5)
REAL(wp) :: r_nenner
REAL(wp) :: pmvc
REAL(wp) :: dtmrt
REAL(wp) :: SQRT_v1m
INTEGER(iwp) :: i
INTEGER(iwp) :: j
INTEGER(iwp) :: i_bin
!-- Coefficient of the 3 regression lines
! 1:const 2:*pmvc 3:*tt2m 4:*SQRT_v1m 5:*dtmrt
coeff(1,1:5) = &
(/ +0.161_wp, +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
coeff(2,1:5) = &
(/ 0.795_wp, 0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp/)
coeff(3,1:5) = &
(/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
!-- Initialise
nerr = 0_iwp
dpmv_cold_res = 0._wp
pmvc = pmva
dtmrt = tmrt - tt2m
SQRT_v1m = vau1m
IF ( SQRT_v1m < 0.10_wp ) THEN
SQRT_v1m = 0.10_wp
ELSE
SQRT_v1m = SQRT ( SQRT_v1m )
ENDIF
DO i = 1, 3
delta_kalt (i) = 0._wp
IF ( i < 3 ) THEN
pmv_cross (i) = pmvc
ENDIF
ENDDO
DO i = 1, 3
!-- Regression constant for given meteorological variables
reg_a(i) = coeff(i, 1) + coeff(i, 3) * tt2m + coeff(i, 4) * SQRT_v1m + &
coeff(i,5)*dtmrt
delta_kalt(i) = reg_a(i) + coeff(i, 2) * pmvc
ENDDO
!-- Intersection points of regression lines in terms of Fanger's PMV
DO i = 1, 2
r_nenner = coeff (i, 2) - coeff (i + 1, 2)
IF ( ABS ( r_nenner ) > 0.00001_wp ) THEN
pmv_cross(i) = ( reg_a (i + 1) - reg_a (i) ) / r_nenner
ELSE
nerr = 1_iwp
RETURN
ENDIF
ENDDO
i_bin = 3
DO i = 1, 2
IF ( pmva > pmv_cross (i) ) THEN
i_bin = i
EXIT
ENDIF
ENDDO
!-- Adjust to operative temperature scaled according
! to classical PMV (Fanger)
dpmv_cold_res = delta_kalt(i_bin) - dpmv_adj(pmva)
END SUBROUTINE dpmv_cold
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculates the summand dpmv_adj adjusting to the operative temperature
!> scaled according to classical PMV (Fanger)
!> Reference environment: v_1m = 0.10 m/s, dTMRT = 0 K, r.h. = 50 %
!------------------------------------------------------------------------------!
REAL(wp) FUNCTION dpmv_adj( pmva )
IMPLICIT NONE
REAL(wp), INTENT ( IN ) :: pmva
INTEGER(iwp), PARAMETER :: n_bin = 3
INTEGER(iwp), PARAMETER :: n_ca = 0
INTEGER(iwp), PARAMETER :: n_ce = 3
REAL(wp), dimension (n_bin, n_ca:n_ce) :: coef
REAL(wp) :: pmv
INTEGER(iwp) :: i, i_bin
! range_1 range_2 range_3
DATA (coef(i, 0), i = 1, n_bin) /+0.0941540_wp, -0.1506620_wp, -0.0871439_wp/
DATA (coef(i, 1), i = 1, n_bin) /+0.0783162_wp, -1.0612651_wp, +0.1695040_wp/
DATA (coef(i, 2), i = 1, n_bin) /+0.1350144_wp, -1.0049144_wp, -0.0167627_wp/
DATA (coef(i, 3), i = 1, n_bin) /+0.1104037_wp, -0.2005277_wp, -0.0003230_wp/
IF ( pmva <= -2.1226_wp ) THEN
i_bin = 3_iwp
ELSE IF ( pmva <= -1.28_wp ) THEN
i_bin = 2_iwp
ELSE
i_bin = 1_iwp
ENDIF
dpmv_adj = coef( i_bin, n_ca )
pmv = 1._wp
DO i = n_ca + 1, n_ce
pmv = pmv * pmva
dpmv_adj = dpmv_adj + coef(i_bin, i) * pmv
ENDDO
RETURN
END FUNCTION dpmv_adj
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Based on perceived temperature (gt) as input, ireq_neutral determines
!> the required clothing insulation (clo) for thermally neutral conditions
!> (neither body cooling nor body heating). It is related to the Klima-
!> Michel activity level (134.682 W/m2). IREQ_neutral is only defined
!> for gt < 10 (�C)
!------------------------------------------------------------------------------!
REAL(wp) FUNCTION ireq_neutral( gt, ireq_minimal, nerr )
IMPLICIT NONE
!-- Type declaration of arguments
REAL(wp), INTENT ( IN ) :: gt
REAL(wp), INTENT ( OUT ) :: ireq_minimal
INTEGER(iwp), INTENT ( OUT ) :: nerr
!-- Type declaration for internal varables
REAL(wp) :: top02
!-- Initialise
nerr = 0_iwp
!-- Convert perceived temperature from basis 0.1 m/s to basis 0.2 m/s
top02 = 1.8788_wp + 0.9296_wp * gt
!-- IREQ neutral conditions (thermal comfort)
ireq_neutral = 1.62_wp - 0.0564_wp * top02
!-- Regression only defined for gt <= 10 (�C)
IF ( ireq_neutral < 0.5_wp ) THEN
IF ( ireq_neutral < 0.48_wp ) THEN
nerr = 1_iwp
ENDIF
ireq_neutral = 0.5_wp
ENDIF
!-- Minimal required clothing insulation: maximal acceptable body cooling
ireq_minimal = 1.26_wp - 0.0588_wp * top02
IF ( nerr > 0_iwp ) THEN
ireq_minimal = ireq_neutral
ENDIF
RETURN
END FUNCTION ireq_neutral
END MODULE biometeorology_pt_mod