source: palm/trunk/SOURCE/biometeorology_ipt_mod.f90 @ 3499

!> @file biometeorology_ipt_mod.f90
!--------------------------------------------------------------------------------!
! This file is part of PALM-4U.
!
! PALM-4U is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2019, Deutscher Wetterdienst (DWD) / 
! German Meteorological Service (DWD) and
! Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions: 001
! -----------------
! 
! 
! Former revisions: 001
! -----------------
! $Id$
! Initial revision 001
!
! 
!
! Authors:
! --------
! @author Dominik Froehlich <dominik.froehlich@mailbox.org>
!
!
! Description:
! ------------
!> Module for the calculation of the dynamic thermal index iPT
!>
!> @todo 
!>
!> @note nothing now
!>
!> @bug  no known bugs by now
!------------------------------------------------------------------------------!
 MODULE biometeorology_ipt_mod
!
!-- Load required variables from existing modules
    USE kinds  !< to set precision of INTEGER and REAL arrays according to PALM

    USE biometeorology_pt_mod,                                                 &
       ONLY: saturation_vapor_pressure, deltapmv, dpmv_adj, dpmv_cold,         &
             calc_sultr, pt_regression, ireq_neutral, iso_ridder

    IMPLICIT NONE

    PRIVATE  !-- No private interfaces --!

    PUBLIC ipt_cycle, ipt_init

!-- PALM interfaces:
    INTERFACE ipt_cycle
       MODULE PROCEDURE ipt_cycle
    END INTERFACE ipt_cycle

    INTERFACE ipt_init
       MODULE PROCEDURE ipt_init
    END INTERFACE ipt_init


 CONTAINS



!------------------------------------------------------------------------------!
! Description:
! ------------
!> The SUBROUTINE surface area calculates the surface area of the individual
!> according to its height (m), weight (kg), and age (y)
!------------------------------------------------------------------------------!
 SUBROUTINE surface_area ( height_cm, weight, age, surf )

    IMPLICIT NONE

    REAL(wp)    , INTENT(in)  :: weight
    REAL(wp)    , INTENT(in)  :: height_cm
    INTEGER(iwp), INTENT(in)  :: age
    REAL(wp)    , INTENT(out) :: surf
    REAL(wp)                  :: height

    height = height_cm * 100._wp

!-- According to Gehan-George, for children
    IF ( age < 19_iwp ) THEN
       IF ( age < 5_iwp ) THEN
          surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp
          RETURN
       END IF
       surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp
       RETURN
    END IF

!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if 
!   height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
    surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
    RETURN

 END SUBROUTINE surface_area

!------------------------------------------------------------------------------!
! Description:
! ------------
!> The SUBROUTINE persdat calculates
!>  - the total internal energy production = metabolic + workload,
!>  - the total internal energy production for a standardized surface (actlev)
!>  - the DuBois - area (a_surf [m2])
!> from
!>  - the persons age (years),
!>  - weight (kg),
!>  - height (m),
!>  - sex (1 = male, 2 = female),
!>  - work load (W)
!> for a sample human.
!------------------------------------------------------------------------------!
 SUBROUTINE persdat ( age, weight, height, sex, work, a_surf, actlev )
!
    IMPLICIT NONE

    REAL(wp), INTENT(in) ::  age
    REAL(wp), INTENT(in) ::  weight
    REAL(wp), INTENT(in) ::  height
    REAL(wp), INTENT(in) ::  work
    INTEGER(iwp), INTENT(in) ::  sex
    REAL(wp), INTENT(out) ::  actlev
    REAL(wp) ::  a_surf
    REAL(wp) ::  energieumsatz
    REAL(wp) ::  s
    REAL(wp) ::  factor
    REAL(wp) ::  heightundumsatz

    CALL surface_area ( height, weight, INT( age ), a_surf )
    s = height * 100._wp / ( weight ** ( 1._wp / 3._wp ) )
    factor = 1._wp + .004_wp  * ( 30._wp - age )
    heightundumsatz = 0.
    IF ( sex == 1_iwp ) THEN
       heightundumsatz = 3.45_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
                     .01_wp  * ( s - 43.4_wp ) )
    ELSE IF ( sex == 2_iwp ) THEN
       heightundumsatz = 3.19_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
                    .018_wp * ( s - 42.1_wp ) )
    END IF

    energieumsatz = work + heightundumsatz
    actlev = energieumsatz / a_surf

 END SUBROUTINE persdat


!------------------------------------------------------------------------------!
! Description:
! ------------
!> SUBROUTINE ipt_init
!> initializes the instationary perceived temperature
!------------------------------------------------------------------------------!

 SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo,             &
     ta, vp, vau1m, tmrt, pb, dt, storage, t_clothing,                         &
     ipt )

    IMPLICIT NONE

!-- Input parameters
    REAL(wp), INTENT(in) ::  age        !< Persons age          (years)
    REAL(wp), INTENT(in) ::  weight     !< Persons weight       (kg)
    REAL(wp), INTENT(in) ::  height     !< Persons height       (m)
    REAL(wp), INTENT(in) ::  work       !< Current workload     (W)
    REAL(wp), INTENT(in) ::  ta         !< Air Temperature      (°C)
    REAL(wp), INTENT(in) ::  vp         !< Vapor pressure       (hPa)
    REAL(wp), INTENT(in) ::  vau1m      !< Wind speed in approx. 1.1m (m/s)
    REAL(wp), INTENT(in) ::  tmrt       !< Mean radiant temperature   (°C)
    REAL(wp), INTENT(in) ::  pb         !< Air pressure         (hPa)
    REAL(wp), INTENT(in) ::  dt         !< Timestep             (s)
    INTEGER(iwp), INTENT(in)  :: sex    !< Persons sex (1 = male, 2 = female)

!-- Output parameters
    REAL(wp), INTENT(out) ::  actlev
    REAL(wp), INTENT(out) ::  clo
    REAL(wp), INTENT(out) ::  storage
    REAL(wp), INTENT(out) ::  t_clothing
    REAL(wp), INTENT(out) ::  ipt

!-- Internal variables
    REAL(wp), PARAMETER :: eps = 0.0005_wp
    REAL(wp), PARAMETER :: eta = 0._wp
    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
    REAL(wp) ::  a_surf
    REAL(wp) ::  acti
    INTEGER(iwp) ::  nzaehl
    INTEGER(iwp) ::  nerr_kalt
    INTEGER(iwp) ::  nerr

    LOGICAL ::  sultrieness

    storage = 0._wp
    CALL persdat ( age, weight, height, sex, work, a_surf, actlev )

!-- Initialise
    t_clothing = -999.0_wp
    ipt        = -999.0_wp
    nerr       = 0_wp
    nzaehl     = 0_wp
    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 ( ta <= 10._wp ) THEN

!--    First guess: winter clothing insulation: cold stress
       clo = oclo
       t_clothing = -999.0_wp  ! force initial run
       CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,        &
          t_clothing, storage, dt, pmva )
       pmv_o = pmva

       IF ( pmva > 0._wp ) THEN

!--       Case summer clothing insulation: heat load ?            
          clo = uclo
          t_clothing = -999.0_wp  ! force initial run
          CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,     &
             t_clothing, storage, dt, pmva )
          pmv_u = pmva
          IF ( pmva <= 0._wp ) THEN
!--          Case: comfort achievable by varying clothing insulation 
!            between winter and summer set values
             CALL iso_ridder ( ta, tmrt, vp, 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
             t_clothing = -999.0_wp
             CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,  &
                t_clothing, storage, dt, pmva )
          ENDIF
       ELSE IF ( pmva < -0.11_wp ) THEN
          clo = 1.75_wp
          t_clothing = -999.0_wp
          CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,     &
             t_clothing, storage, dt, pmva )
       ENDIF

    ELSE

!--    First guess: summer clothing insulation: heat load
       clo = uclo
       t_clothing = -999.0_wp
       CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,        &
          t_clothing, storage, dt, pmva )
       pmv_u = pmva

       IF ( pmva < 0._wp ) THEN

!--       Case winter clothing insulation: cold stress ?
          clo = oclo
          t_clothing = -999.0_wp
          CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,     &
             t_clothing, storage, dt, pmva )
          pmv_o = pmva

          IF ( pmva >= 0._wp ) THEN

!--          Case: comfort achievable by varying clothing insulation 
!            between winter and summer set values
             CALL iso_ridder ( ta, tmrt, vp, vau1m, pb, actlev, eta, uclo,     &
                               pmv_u, oclo, pmv_o, eps, pmva, top, nzaehl, clo )
             IF ( nzaehl < 0_wp ) THEN
                nerr = -1_iwp
                RETURN
             ENDIF
          ELSE IF ( pmva < -0.11_wp ) THEN
             clo = 1.75_wp
             t_clothing = -999.0_wp
             CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,  &
                t_clothing, storage, dt, pmva )
          ENDIF
       ELSE IF ( pmva > 0.06_wp ) THEN
          clo = 0.5_wp
          t_clothing = -999.0_wp
          CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,     &
             t_clothing, storage, dt, pmva )
       ENDIF

    ENDIF

!-- Determine perceived temperature by regression equation + adjustments
    pmvs = pmva
    CALL pt_regression ( pmva, clo, ipt )
    ptc = ipt
    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
!--    Adjust for cold conditions according to Gagge 1986
       CALL dpmv_cold ( pmva, ta, 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, ipt )
    ENDIF
    clo_fanger = clo
    clon = clo
    IF ( clo > 0.5_wp .AND. ipt <= 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 ( ipt, 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 ( ta, svp_tt )
       d_pmv  = deltapmv ( pmva, ta, vp, svp_tt, tmrt, vau1m, nerr )
       pmvs   = pmva + d_pmv
       CALL pt_regression ( pmvs, clo, ipt )
       IF ( sult_lim < 99._wp ) THEN
          IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
       ENDIF
    ENDIF

 
 END SUBROUTINE ipt_init
 
!------------------------------------------------------------------------------!
! Description:
! ------------
!> SUBROUTINE ipt_cycle
!> Calculates one timestep for the instationary version of perceived
!> temperature (iPT, °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)
!>
!------------------------------------------------------------------------------!
 SUBROUTINE ipt_cycle( ta, vp, vau1m, tmrt, pb, dt, storage, t_clothing, clo,  &
     actlev, work, ipt )

    IMPLICIT NONE

!-- Type of input of the argument list
    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature             (°C)
    REAL(wp), INTENT ( IN )  ::  vp      !< Vapor pressure              (hPa)
    REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature    (°C)
    REAL(wp), INTENT ( IN )  ::  vau1m   !< Wind speed                  (m/s)
    REAL(wp), INTENT ( IN )  ::  pb      !< Air pressure                (hPa)
    REAL(wp), INTENT ( IN )  ::  dt      !< Timestep                    (s)
    REAL(wp), INTENT ( IN )  ::  clo     !< Clothing index              (no dim)
    REAL(wp), INTENT ( IN )  ::  actlev  !< Internal heat production    (W)
    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical work load        (W)

!-- In and output parameters
    REAL(wp), INTENT (INOUT) ::  storage     !< Heat storage            (W)
    REAL(wp), INTENT (INOUT) ::  t_clothing  !< Clothig temperature     (°C)

!-- Type of output of the argument list
    REAL(wp), INTENT ( OUT ) ::  ipt  !< Instationary perceived temperature (°C)

!-- Type of program variables
    REAL(wp), PARAMETER ::  eps = 0.0005_wp
    REAL(wp) ::  d_pmv
    REAL(wp) ::  pmv_o
    REAL(wp) ::  pmv_s
    REAL(wp) ::  storage_rate
    REAL(wp) ::  ener_in
    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_u
    REAL(wp) ::  half_hour_frac
    REAL(wp) ::  storage_delta
    REAL(wp) ::  storage_delta_max
    REAL(wp) ::  pmva
    REAL(wp) ::  ptc
    REAL(wp) ::  d_std
    REAL(wp) ::  pmvs
    REAL(wp) ::  top
    INTEGER(iwp) ::  nzaehl
    INTEGER(iwp) ::  nerr_kalt
    INTEGER(iwp) ::  nerr

    LOGICAL ::  sultrieness

!-- Initialise
    ipt = -999.0_wp

    nerr     = 0_iwp
    nzaehl   = 0_iwp
    sultrieness  = .FALSE.

!-- Determine pmv_adjusted for current conditions
    CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,           &
       t_clothing, storage, dt, pmva )
    pmv_s = pmva  !< PMV considering storage

!-- Determine perceived temperature by regression equation + adjustments
    CALL pt_regression ( pmva, clo, ipt )

    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
!--    Adjust for cold conditions according to Gagge 1986
       CALL dpmv_cold ( pmva, ta, 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, ipt )
    ENDIF

!-- Consider sultriness if appropriate
    ptc = ipt
    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 ( ta, svp_tt )
       d_pmv  = deltapmv ( pmva, ta, vp, svp_tt, tmrt, vau1m, nerr )
       pmvs   = pmva + d_pmv
       CALL pt_regression ( pmvs, clo, ipt )
       IF ( sult_lim < 99._wp ) THEN
          IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
       ENDIF
    ENDIF

 END SUBROUTINE ipt_cycle

!------------------------------------------------------------------------------!
! Description:
! ------------
!> SUBROUTINE fanger_s calculates the 
!> actual Predicted Mean Vote (dimensionless) according 
!> to Fanger corresponding to meteorological (tt,tmrt,pa,vau,pb) 
!> and individual variables (clo, actlev, eta) considering a storage
!> and clothing temperature for a given timestep.
!------------------------------------------------------------------------------!
 SUBROUTINE fanger_s_acti( tt, tmrt, pa, in_vau, pb, in_clo, actlev,           &
    activity, t_cloth, s, dt, pmva )

    IMPLICIT NONE

!--  Input argument types
    REAL(wp), INTENT ( IN )  ::  tt       !< Air temperature          (°C)
    REAL(wp), INTENT ( IN )  ::  tmrt     !< Mean radiant temperature (°C)
    REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
    REAL(wp), INTENT ( IN )  ::  pb       !< Air pressure             (hPa)
    REAL(wp), INTENT ( IN )  ::  in_vau   !< Wind speed               (m/s)
    REAL(wp), INTENT ( IN )  ::  actlev   !< Metabolic + work energy  (W/m²)
    REAL(wp), INTENT ( IN )  ::  dt       !< Timestep                 (s)
    REAL(wp), INTENT ( IN )  ::  activity !< Work load                (W/m²)
    REAL(wp), INTENT ( IN )  ::  in_clo   !< Clothing index (clo)     (no dim)

!-- Output argument types
    REAL(wp), INTENT ( OUT ) ::  pmva  !< actual Predicted Mean Vote (no dim)

    REAL(wp), INTENT (INOUT) ::  s  !< storage var. of energy balance (W/m2)
    REAL(wp), INTENT (INOUT) ::  t_cloth  !< clothing temperature (°C)

!-- Internal variables
    REAL(wp) ::  f_cl
    REAL(wp) ::  heat_convection
    REAL(wp) ::  t_skin_aver
    REAL(wp) ::  bc
    REAL(wp) ::  cc
    REAL(wp) ::  dc
    REAL(wp) ::  ec
    REAL(wp) ::  gc
    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
    REAL(wp) ::  en
    REAL(wp) ::  d_s
    REAL(wp) ::  t_clothing   !< Clothing index value   (dimensionless)
    REAL(wp) ::  adjustrate
    REAL(wp) ::  adjustrate_cloth
    REAL(wp), PARAMETER  ::  time_equil = 7200._wp
    INTEGER(iwp) ::  i
    INTEGER(iwp) ::  niter  !< Running index


!-- Clo must be > 0. to avoid div. by 0!
    clo = in_clo
    IF ( clo < 001._wp ) clo = .001_wp

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

!-- 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._wp ) * f_cl
    cc = f_cl * heat_convection
    ec = .155_wp * clo
    dc = ( 1._wp + ec * cc ) / bc
    gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4._wp + ec * cc * tt ) / bc

!-- Calculation of clothing surface temperature (t_clothing) based on
!   newton-approximation with air temperature as initial guess
    niter = dt
    adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
    IF ( adjustrate >= 1._wp ) adjustrate = 1._wp
    adjustrate_cloth = adjustrate * 30._wp
    t_clothing = t_cloth

    IF ( t_cloth <= -998.0_wp ) THEN  ! If initial run
       niter = 3_wp
       adjustrate = 1._wp
       adjustrate_cloth = 1._wp
       t_clothing = tt
    ENDIF

    DO i = 1, niter
       t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
          273.2_wp )**4._wp + t_clothing *                                     &
          dc - gc ) / ( 4._wp * ( t_clothing + 273.2_wp )**3._wp + 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 + 273.2_wp )**4 - ( tmrt +         &
       273.2_wp )**4 )

!-- Heat loss through forced convection
    z6 = cc * ( t_clothing - tt )

!-- Write together as energy ballance
    en = activity - z2 - z3 - z4 - z5 - z6

!-- Manage storage
    d_s = adjustrate * en + ( 1._wp - adjustrate ) * s

!-- Predicted Mean Vote
    pmva = z1 * d_s

!-- Update storage
    s = d_s
    t_cloth = t_clothing

 END SUBROUTINE fanger_s_acti


 END MODULE biometeorology_ipt_mod