!> @file biometeorology_PET.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 2019, Deutscher Wetterdienst (DWD) /
! German Meteorological Service (DWD)
!------------------------------------------------------------------------------!
!
! Current revisions: 003
! -----------------
!
!
! Former revisions: 001
! -----------------
! $Id$
! Initial revision 001
!
!
!
! Authors:
! --------
! @author Peter Höppe (original author)
! @author Dominik Froehlich (PALM modification)
!
!------------------------------------------------------------------------------!
MODULE biometeorology_pet_mod
!-- Load required variables from existing modules
USE kinds !< to set precision of INTEGER and REAL according to PALM
IMPLICIT NONE
PRIVATE
!-- Internal constants:
REAL(wp), PARAMETER :: po = 1013.25_wp !< Air pressure at sea level (hPa)
REAL(wp), PARAMETER :: rob = 1.06_wp !<
REAL(wp), PARAMETER :: cb = 3640._wp !<
REAL(wp), PARAMETER :: food = 0._wp !< Heat gain by food (W)
REAL(wp), PARAMETER :: emsk = 0.99_wp !< Longwave emission coef. of skin
REAL(wp), PARAMETER :: emcl = 0.95_wp !< Longwave emission coef. of cloth
REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp !<
REAL(wp), PARAMETER :: sigm = 5.67_wp * 10._wp **(-8._wp) !< Stefan-Boltzmann-Const.
REAL(wp), PARAMETER :: cair = 1010._wp !<
REAL(wp), PARAMETER :: cels_offs = 273.15_wp !< Kelvin Celsius Offset (K)
!-- Internal variables:
REAL(wp) :: h !< Internal heat (W)
!-- MEMI configuration
REAL(wp) :: age !< Persons age (a)
REAL(wp) :: mbody !< Persons body mass (kg)
REAL(wp) :: ht !< Persons height (m)
REAL(wp) :: work !< Work load (W)
REAL(wp) :: eta !< Work efficiency (dimensionless)
REAL(wp) :: icl !< Clothing insulation index (clo)
REAL(wp) :: fcl !< Surface area modification by clothing (factor)
INTEGER(iwp) :: pos !< Posture: 1 = standing, 2 = sitting
INTEGER(iwp) :: sex !< Sex: 1 = male, 2 = female
PUBLIC calculate_pet_static
!-- PALM interfaces:
INTERFACE calculate_pet_static
MODULE PROCEDURE calculate_pet_static
END INTERFACE calculate_pet_static
CONTAINS
!------------------------------------------------------------------------------!
!
! Description:
! ------------
!> Physiologically Equivalent Temperature (PET),
! stationary (calculated based on MEMI),
! Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
!------------------------------------------------------------------------------!
SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, p, tx )
IMPLICIT NONE
!-- Input arguments:
REAL(wp), INTENT( IN ) :: ta !< Air temperature (°C)
REAL(wp), INTENT( IN ) :: tmrt !< Mean radiant temperature (°C)
REAL(wp), INTENT( IN ) :: v !< Wind speed (m/s)
REAL(wp), INTENT( IN ) :: vpa !< Vapor pressure (hPa)
REAL(wp), INTENT( IN ) :: p !< Air pressure (hPa)
!-- Output arguments:
REAL(wp), INTENT ( OUT ) :: tx !< PET (°C)
!-- Internal variables:
REAL(wp) :: acl
REAL(wp) :: adu
REAL(wp) :: aeff
REAL(wp) :: ere
REAL(wp) :: erel
REAL(wp) :: esw
REAL(wp) :: facl
REAL(wp) :: feff
REAL(wp) :: rdcl
REAL(wp) :: rdsk
REAL(wp) :: rtv
REAL(wp) :: vpts
REAL(wp) :: tsk
REAL(wp) :: tcl
REAL(wp) :: wetsk
!-- Configuration
age = 35._wp
mbody = 75._wp
ht = 1.75_wp
work = 80._wp
eta = 0._wp
icl = 0.9_wp
fcl = 1.15_wp
!-- Call subfunctions
CALL in_body ( ere, erel, p, rtv, ta, vpa )
CALL heat_exch ( acl, adu, aeff, ere, erel, esw, facl, feff, &
p, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, p, rdcl, &
rdsk, rtv, ta, tcl, tsk, tx, vpts, wetsk )
END SUBROUTINE calculate_pet_static
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate internal energy ballance
!------------------------------------------------------------------------------!
SUBROUTINE in_body ( ere, erel, p, rtv, ta, vpa )
!-- Input arguments:
REAL(wp), INTENT( IN ) :: p !< air pressure (hPa)
REAL(wp), INTENT( IN ) :: ta !< air temperature (°C)
REAL(wp), INTENT( IN ) :: vpa !< vapor pressure (hPa)
!-- Output arguments:
REAL(wp), INTENT( OUT ) :: ere !< energy ballance (W)
REAL(wp), INTENT( OUT ) :: erel !< latent energy ballance (W)
REAL(wp), INTENT( OUT ) :: rtv !<
!-- Internal variables:
REAL(wp) :: eres
REAL(wp) :: met
REAL(wp) :: tex
REAL(wp) :: vpex
met = 3.45_wp * mbody ** ( 3._wp / 4._wp ) * (1._wp + 0.004_wp * &
( 30._wp - age) + 0.010_wp * ( ( ht * 100._wp / &
( mbody ** ( 1._wp / 3._wp ) ) ) - 43.4_wp ) )
met = work + met
h = met * (1._wp - eta)
!-- Sensible respiration energy
tex = 0.47_wp * ta + 21.0_wp
rtv = 1.44_wp * 10._wp ** (-6._wp) * met
eres = cair * (ta - tex) * rtv
!-- Latent respiration energy
vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
erel = 0.623_wp * evap / p * ( vpa - vpex ) * rtv
!-- Sum of the results
ere = eres + erel
END SUBROUTINE in_body
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate heat gain or loss
!------------------------------------------------------------------------------!
SUBROUTINE heat_exch ( acl, adu, aeff, ere, erel, esw, facl, feff, &
p, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
!-- Input arguments:
REAL(wp), INTENT( IN ) :: ere !< Energy ballance (W)
REAL(wp), INTENT( IN ) :: erel !< Latent energy ballance (W)
REAL(wp), INTENT( IN ) :: p !< Air pressure (hPa)
REAL(wp), INTENT( IN ) :: ta !< Air temperature (°C)
REAL(wp), INTENT( IN ) :: tmrt !< Mean radiant temperature (°C)
REAL(wp), INTENT( IN ) :: v !< Wind speed (m/s)
REAL(wp), INTENT( IN ) :: vpa !< Vapor pressure (hPa)
!-- Output arguments:
REAL(wp), INTENT( OUT ) :: acl !< Clothing surface area (m²)
REAL(wp), INTENT( OUT ) :: adu !< Du-Bois area (m²)
REAL(wp), INTENT( OUT ) :: aeff !< Effective surface area (m²)
REAL(wp), INTENT( OUT ) :: esw !< Energy-loss through sweat evap. (W)
REAL(wp), INTENT( OUT ) :: facl !< Surface area extension through clothing (factor)
REAL(wp), INTENT( OUT ) :: feff !< Surface modification by posture (factor)
REAL(wp), INTENT( OUT ) :: rdcl !< Diffusion resistence of clothing (factor)
REAL(wp), INTENT( OUT ) :: rdsk !< Diffusion resistence of skin (factor)
REAL(wp), INTENT( OUT ) :: tcl !< Clothing temperature (°C)
REAL(wp), INTENT( OUT ) :: tsk !< Skin temperature (°C)
REAL(wp), INTENT( OUT ) :: vpts !< Sat. vapor pressure over skin (hPa)
REAL(wp), INTENT( OUT ) :: wetsk !< Fraction of wet skin (dimensionless)
!-- Internal variables
REAL(wp) :: c(0:10)
REAL(wp) :: cbare
REAL(wp) :: cclo
REAL(wp) :: csum
REAL(wp) :: di
REAL(wp) :: ed
REAL(wp) :: enbal
REAL(wp) :: enbal2
REAL(wp) :: eswdif
REAL(wp) :: eswphy
REAL(wp) :: eswpot
REAL(wp) :: fec
REAL(wp) :: hc
REAL(wp) :: he
REAL(wp) :: htcl
REAL(wp) :: r1
REAL(wp) :: r2
REAL(wp) :: rbare
REAL(wp) :: rcl
REAL(wp) :: rclo
REAL(wp) :: rclo2
REAL(wp) :: rsum
REAL(wp) :: sw
REAL(wp) :: swf
REAL(wp) :: swm
REAL(wp) :: tbody
REAL(wp) :: tcore(1:7)
REAL(wp) :: vb
REAL(wp) :: vb1
REAL(wp) :: vb2
REAL(wp) :: wd
REAL(wp) :: wr
REAL(wp) :: ws
REAL(wp) :: wsum
REAL(wp) :: xx
REAL(wp) :: y
INTEGER(iwp) :: count1
INTEGER(iwp) :: count3
INTEGER(iwp) :: j
INTEGER(iwp) :: i
LOGICAL :: skipIncreaseCount
wetsk = 0._wp
adu = 0.203_wp * mbody ** 0.425_wp * ht ** 0.725_wp
hc = 2.67_wp + ( 6.5_wp * v ** 0.67_wp)
hc = hc * (p /po) ** 0.55_wp
feff = 0.725_wp !< Posture: 0.725 for stading
facl = (- 2.36_wp + 173.51_wp * icl - 100.76_wp * icl * icl + 19.28_wp &
* (icl ** 3._wp)) / 100._wp
IF ( facl > 1._wp ) facl = 1._wp
rcl = ( icl / 6.45_wp ) / facl
IF ( icl >= 2._wp ) y = 1._wp
IF ( ( icl > 0.6_wp ) .AND. ( icl < 2._wp ) ) y = ( ht - 0.2_wp ) / ht
IF ( ( icl <= 0.6_wp ) .AND. ( icl > 0.3_wp ) ) y = 0.5_wp
IF ( ( icl <= 0.3_wp ) .AND. ( icl > 0._wp ) ) y = 0.1_wp
r2 = adu * (fcl - 1._wp + facl) / (2._wp * 3.14_wp * ht * y)
r1 = facl * adu / (2._wp * 3.14_wp * ht * y)
di = r2 - r1
!-- Skin temperatur
DO j = 1, 7
tsk = 34._wp
count1 = 0_iwp
tcl = ( ta + tmrt + tsk ) / 3._wp
count3 = 1_iwp
enbal2 = 0._wp
DO i = 1, 100 ! allow for 100 iterations max
acl = adu * facl + adu * ( fcl - 1._wp )
rclo2 = emcl * sigm * ( ( tcl + cels_offs )**4._wp - &
( tmrt + cels_offs )** 4._wp ) * feff
htcl = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
tsk = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
!-- Radiation saldo
aeff = adu * feff
rbare = aeff * ( 1._wp - facl ) * emsk * sigm * &
( ( tmrt + cels_offs )** 4._wp - ( tsk + cels_offs )** 4._wp )
rclo = feff * acl * emcl * sigm * &
( ( tmrt + cels_offs )** 4._wp - ( tcl + cels_offs )** 4._wp )
rsum = rbare + rclo
!-- Convection
cbare = hc * ( ta - tsk ) * adu * ( 1._wp - facl )
cclo = hc * ( ta - tcl ) * acl
csum = cbare + cclo
!-- Core temperature
c(0) = h + ere
c(1) = adu * rob * cb
c(2) = 18._wp - 0.5_wp * tsk
c(3) = 5.28_wp * adu * c(2)
c(4) = 0.0208_wp * c(1)
c(5) = 0.76075_wp * 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._wp * c(4) * c(7)
c(9) = 5.28_wp * adu - c(5) - c(4) * tsk
c(10) = c(9) * c(9) - 4._wp * c(4) * &
( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
IF ( tsk == 36._wp ) tsk = 36.01_wp
tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk
tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) / &
( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk
IF ( c(10) >= 0._wp ) THEN
tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2._wp * c(4) )
tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2._wp * c(4) )
END IF
IF ( c(8) >= 0._wp ) THEN
tcore(2) = ( - c(6) + ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) )
tcore(5) = ( - c(6) - ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) )
tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk
END IF
!-- Transpiration
tbody = 0.1_wp * tsk + 0.9_wp * tcore(j)
swm = 304.94_wp * ( tbody - 36.6_wp ) * adu / 3600000._wp
vpts = 6.11_wp * 10._wp**( 7.45_wp * tsk / ( 235._wp + tsk ) )
IF ( tbody <= 36.6_wp ) swm = 0._wp
sw = swm
eswphy = - sw * evap
he = 0.633_wp * hc / ( p * cair )
fec = 1._wp / ( 1._wp + 0.92_wp * hc * rcl )
eswpot = he * ( vpa - vpts ) * adu * evap * fec
wetsk = eswphy / eswpot
IF ( wetsk > 1._wp ) wetsk = 1._wp
eswdif = eswphy - eswpot
IF ( eswdif <= 0._wp ) esw = eswpot
IF ( eswdif > 0._wp ) esw = eswphy
IF ( esw > 0._wp ) esw = 0._wp
!-- Diffusion
rdsk = 0.79_wp * 10._wp ** 7._wp
rdcl = 0._wp
ed = evap / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa - &
vpts )
!-- Max vb
vb1 = 34._wp - tsk
vb2 = tcore(j) - 36.6_wp
IF ( vb2 < 0._wp ) vb2 = 0._wp
IF ( vb1 < 0._wp ) vb1 = 0._wp
vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 )
!-- Energy ballence
enbal = h + ed + ere + esw + csum + rsum + food
!-- Clothing temperature
xx = 0.001_wp
IF ( count1 == 0_iwp ) xx = 1._wp
IF ( count1 == 1_iwp ) xx = 0.1_wp
IF ( count1 == 2_iwp ) xx = 0.01_wp
IF ( count1 == 3_iwp ) xx = 0.001_wp
IF ( enbal > 0._wp ) tcl = tcl + xx
IF ( enbal < 0._wp ) tcl = tcl - xx
skipIncreaseCount = .FALSE.
IF ( ( (enbal <= 0._wp ) .AND. (enbal2 > 0._wp ) ) .OR. &
( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) ) THEN
skipIncreaseCount = .TRUE.
ELSE
enbal2 = enbal
count3 = count3 + 1_iwp
END IF
IF ( ( count3 > 200_iwp ) .OR. skipIncreaseCount ) THEN
IF ( count1 < 3_iwp ) THEN
count1 = count1 + 1_iwp
enbal2 = 0._wp
ELSE
EXIT
END IF
END IF
END DO
IF ( count1 == 3_iwp ) THEN
SELECT CASE ( j )
CASE ( 2, 5)
IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND. &
( tsk <= 34.050_wp ) ) ) CYCLE
CASE ( 6, 1 )
IF ( c(10) < 0._wp ) CYCLE
IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND. &
( tsk > 33.850_wp ) ) ) CYCLE
CASE ( 3 )
IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND. &
( tsk <= 34.000_wp ) ) ) CYCLE
CASE ( 7 )
IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND. &
( tsk > 34.000_wp ) ) ) CYCLE
CASE default
END SELECT
END IF
IF ( ( j /= 4_iwp ) .AND. ( vb >= 91._wp ) ) CYCLE
IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) ) CYCLE
IF ( vb > 90._wp ) vb = 90._wp
!-- Loses by water
ws = sw * 3600._wp * 1000._wp
IF ( ws > 2000._wp ) ws = 2000._wp
wd = ed / evap * 3600._wp * ( -1000._wp )
wr = erel / evap * 3600._wp * ( -1000._wp )
wsum = ws + wr + wd
RETURN
END DO
END SUBROUTINE heat_exch
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Calculate PET
!------------------------------------------------------------------------------!
SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, p, rdcl, &
rdsk, rtv, ta, tcl, tsk, tx, vpts, wetsk )
!-- Input arguments:
REAL(wp), INTENT( IN ) :: acl !< clothing surface area (m²)
REAL(wp), INTENT( IN ) :: adu !< Du-Bois area (m²)
REAL(wp), INTENT( IN ) :: esw !< energy-loss through sweat evap. (W)
REAL(wp), INTENT( IN ) :: facl !< surface area extension through clothing (factor)
REAL(wp), INTENT( IN ) :: feff !< surface modification by posture (factor)
REAL(wp), INTENT( IN ) :: p !< air pressure (hPa)
REAL(wp), INTENT( IN ) :: rdcl !< diffusion resistence of clothing (factor)
REAL(wp), INTENT( IN ) :: rdsk !< diffusion resistence of skin (factor)
REAL(wp), INTENT( IN ) :: rtv !<
REAL(wp), INTENT( IN ) :: ta !< air temperature (°C)
REAL(wp), INTENT( IN ) :: tcl !< clothing temperature (°C)
REAL(wp), INTENT( IN ) :: tsk !< skin temperature (°C)
REAL(wp), INTENT( IN ) :: vpts !< sat. vapor pressure over skin (hPa)
REAL(wp), INTENT( IN ) :: wetsk !< fraction of wet skin (dimensionless)
!-- Output arguments:
REAL(wp), INTENT( OUT ) :: aeff !< effective surface area (m²)
REAL(wp), INTENT( OUT ) :: tx !< PET (°C)
!-- Internal variables
REAL ( wp ) :: cbare
REAL ( wp ) :: cclo
REAL ( wp ) :: csum
REAL ( wp ) :: ed
REAL ( wp ) :: enbal
REAL ( wp ) :: enbal2
REAL ( wp ) :: ere
REAL ( wp ) :: erel
REAL ( wp ) :: eres
REAL ( wp ) :: hc
REAL ( wp ) :: rbare
REAL ( wp ) :: rclo
REAL ( wp ) :: rsum
REAL ( wp ) :: tex
REAL ( wp ) :: vpex
REAL ( wp ) :: xx
INTEGER ( iwp ) :: count1
INTEGER ( iwp ) :: i
tx = ta
enbal2 = 0._wp
DO count1 = 0, 3
DO i = 1, 125 ! 500 / 4
hc = 2.67_wp + 6.5_wp * 0.1_wp ** 0.67_wp
hc = hc * ( p / po ) ** 0.55_wp
!-- Radiation
aeff = adu * feff
rbare = aeff * ( 1._wp - facl ) * emsk * sigm * &
( ( tx + cels_offs ) ** 4._wp - ( tsk + cels_offs ) ** 4._wp )
rclo = feff * acl * emcl * sigm * &
( ( tx + cels_offs ) ** 4._wp - ( tcl + cels_offs ) ** 4._wp )
rsum = rbare + rclo
!-- Covection
cbare = hc * ( tx - tsk ) * adu * ( 1._wp - facl )
cclo = hc * ( tx - tcl ) * acl
csum = cbare + cclo
!-- Diffusion
ed = evap / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( 12._wp - &
vpts )
!-- Respiration
tex = 0.47_wp * tx + 21._wp
eres = cair * ( tx - tex ) * rtv
vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
erel = 0.623_wp * evap / p * ( 12._wp - vpex ) * rtv
ere = eres + erel
!-- Energy ballance
enbal = h + ed + ere + esw + csum + rsum
!-- Iteration concerning ta
IF ( count1 == 0_iwp ) xx = 1._wp
IF ( count1 == 1_iwp ) xx = 0.1_wp
IF ( count1 == 2_iwp ) xx = 0.01_wp
IF ( count1 == 3_iwp ) xx = 0.001_wp
IF ( enbal > 0._wp ) tx = tx - xx
IF ( enbal < 0._wp ) tx = tx + xx
IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT
IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT
enbal2 = enbal
END DO
END DO
END SUBROUTINE pet_iteration
END MODULE biometeorology_pet_mod