!> @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