Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 3445)
+++ palm/trunk/SOURCE/Makefile	(revision 3448)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Adjustment of biometeorology dependencies
+# 
+# 3436 2018-10-26 18:35:15Z gronemeier
 # Add surface_mod to user_data_output_mask
 # 
@@ -540,4 +543,8 @@
 	header.f90 \
 	biometeorology_mod.f90 \
+	biometeorology_ipt_mod.f90 \
+	biometeorology_pt_mod.f90 \
+	biometeorology_pet_mod.f90 \
+	biometeorology_utci_mod.f90 \
 	inflow_turbulence.f90 \
 	init_3d_model.f90 \
@@ -755,7 +762,21 @@
 basic_constants_and_equations_mod.o: \
 	mod_kinds.o
+biometeorology_ipt_mod.o: \
+	biometeorology_pt_mod.o \
+	mod_kinds.o
 biometeorology_mod.o: \
-	mod_kinds.o \
-	radiation_model_mod.o
+	modules.o \
+	mod_kinds.o \
+	radiation_model_mod.o \
+	biometeorology_ipt_mod.o \
+	biometeorology_pt_mod.o \
+	biometeorology_pet_mod.o \
+	biometeorology_utci_mod.o
+biometeorology_pt_mod.o: \
+	mod_kinds.o
+biometeorology_pet_mod.o: \
+	mod_kinds.o
+biometeorology_utci_mod.o: \
+	mod_kinds.o
 boundary_conds.o: \
 	basic_constants_and_equations_mod.o \
@@ -919,4 +940,5 @@
 data_output_2d.o: \
 	basic_constants_and_equations_mod.o \
+	biometeorology_mod.o \
 	bulk_cloud_model_mod.o \
 	chemistry_model_mod.o\
@@ -1014,4 +1036,5 @@
 header.o: \
 	basic_constants_and_equations_mod.o \
+	biometeorology_mod.o \
 	bulk_cloud_model_mod.o \
 	chemistry_model_mod.o \
@@ -1042,4 +1065,5 @@
 	advec_ws.o \
 	basic_constants_and_equations_mod.o \
+	biometeorology_mod.o \
 	bulk_cloud_model_mod.o \
 	chem_emissions_mod.o \
@@ -1283,4 +1307,5 @@
 multi_agent_system_mod.o: \
 	basic_constants_and_equations_mod.o \
+	biometeorology_mod.o \
 	chemistry_model_mod.o \
 	cpulog_mod.o \
@@ -1315,6 +1340,6 @@
 	uv_exposure_model_mod.o
 nesting_offl_mod.o: \
-        cpulog_mod.o \
-        mod_kinds.o \
+	cpulog_mod.o \
+	mod_kinds.o \
 	modules.o \
 	netcdf_data_input_mod.o \
@@ -1357,4 +1382,5 @@
 	write_restart_data_mod.o
 parin.o: \
+	biometeorology_mod.o \
 	bulk_cloud_model_mod.o \
 	chemistry_model_mod.o \
@@ -1623,4 +1649,5 @@
 time_integration.o: \
 	advec_ws.o \
+	biometeorology_mod.o \
 	bulk_cloud_model_mod.o \
 	buoyancy.o \
Index: palm/trunk/SOURCE/average_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/average_3d_data.f90	(revision 3445)
+++ palm/trunk/SOURCE/average_3d_data.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Adjustment of biometeorology calls
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Renamed output variables
 ! 
@@ -153,4 +156,7 @@
     USE averaging
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_3d_data_averaging
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model, bcm_3d_data_averaging
@@ -160,6 +166,6 @@
 
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, average_count_3d, doav, doav_n, land_surface,    &
-               ocean_mode, urban_surface, varnamelength
+        ONLY:  air_chemistry, average_count_3d, biometeorology, doav, doav_n,  &
+               land_surface, ocean_mode, urban_surface, varnamelength
 
     USE cpulog,                                                                &
@@ -182,7 +188,4 @@
     USE radiation_model_mod,                                                   &
         ONLY:  radiation, radiation_3d_data_averaging
-
-    USE biometeorology_mod,                                                    &
-        ONLY:  biometeorology_3d_data_averaging
 
     USE turbulence_closure_mod,                                                &
@@ -567,6 +570,6 @@
              ENDIF
 
-             IF ( radiation  .AND.  trimvar(1:4) == 'bio_' )  THEN
-                CALL biometeorology_3d_data_averaging( 'average', doav(ii) )
+             IF ( biometeorology )  THEN
+                CALL biom_3d_data_averaging( 'average', doav(ii) )
              ENDIF
 
Index: palm/trunk/SOURCE/biometeorology_ipt_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_ipt_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/biometeorology_ipt_mod.f90	(revision 3448)
@@ -0,0 +1,615 @@
+!> @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
Index: palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3445)
+++ 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 <http://www.gnu.org/licenses/>.
 !
-! 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 <dominik.froehlich@dwd.de>
 ! @author Jaroslav Resler <resler@cs.cas.cz>
-! @author Dominik Froehlich <dominik.froehlich@dwd.de>
-!
-!------------------------------------------------------------------------------!
-
-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
Index: palm/trunk/SOURCE/biometeorology_pet_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_pet_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/biometeorology_pet_mod.f90	(revision 3448)
@@ -0,0 +1,548 @@
+!> @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 <http://www.gnu.org/licenses/>.
+!
+! 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 <dominik.froehlich@dwd.de> (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
Index: palm/trunk/SOURCE/biometeorology_pt_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_pt_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/biometeorology_pt_mod.f90	(revision 3448)
@@ -0,0 +1,1078 @@
+!> @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 <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 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 = DLOG ( 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 - DLOG ( 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 = DMOD ( 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
Index: palm/trunk/SOURCE/biometeorology_utci_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_utci_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/biometeorology_utci_mod.f90	(revision 3448)
@@ -0,0 +1,330 @@
+!> @file biometeorology_utci_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 2018, Deutscher Wetterdienst (DWD) / 
+! German Meteorological Service (DWD)
+!
+! UTCI regression equation after
+! Bröde P, Fiala D, Blazejczyk K, Holmér I, Jendritzky G, Kampmann B, Tinz B, 
+! Havenith G (2012) Deriving the operational procedure for the Universal Thermal
+! Climate Index (UTCI). International Journal of Biometeorology 56 (3):481-494. 
+! doi:10.1007/s00484-011-0454-1
+!
+! original regression source available at
+! www.utci.org
+!------------------------------------------------------------------------------!
+ MODULE biometeorology_utci_mod
+
+!-- Load required variables from existing modules
+    USE kinds
+
+    IMPLICIT NONE
+
+    PRIVATE
+
+!-- Add INTERFACES that must be available to other modules (alphabetical order)
+    PUBLIC calculate_utci_static
+
+
+    INTERFACE calculate_utci_static
+       MODULE PROCEDURE calculate_utci_static
+    END INTERFACE calculate_utci_static
+
+
+ CONTAINS
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> SUBROUTINE for calculating UTCI Temperature (UTCI)
+!> computed by a 6th order approximation
+!------------------------------------------------------------------------------!
+ SUBROUTINE calculate_utci_static( ta_in, ehpa, v_hag, tmrt, hag, utci )
+
+    IMPLICIT NONE
+!-- Type of input of the argument list
+    REAL(WP), INTENT ( IN )  ::  ta_in  !< Local air temperature (°C)
+    REAL(WP), INTENT ( IN )  ::  ehpa   !< Loacl vapour pressure (hPa)
+    REAL(WP), INTENT ( IN )  ::  v_hag  !< Incident wind speed (m/s)
+    REAL(WP), INTENT ( IN )  ::  tmrt   !< Local mean radiant temperature (°C)
+    REAL(WP), INTENT ( IN )  ::  hag    !< Height of wind speed input (m)
+!-- Type of output of the argument list
+    REAL(wp), INTENT ( OUT ) ::  utci   !< Universal Thermal Climate Index (°C)
+
+!-- Make sure precission is sufficient for regression equation
+    DOUBLE PRECISION ::  ta
+    DOUBLE PRECISION ::  pa
+    DOUBLE PRECISION ::  d_tmrt
+    DOUBLE PRECISION ::  va
+    DOUBLE PRECISION ::  utci_d
+    DOUBLE PRECISION ::  offset
+
+!-- Initialize
+    offset = 0.d0
+    ta = DBLE( ta_in )
+    d_tmrt = DBLE( tmrt - ta_in )
+
+!-- Use vapour pressure in kpa
+    pa = DBLE( ehpa / 10.0_wp )
+
+!-- Wind altitude correction from hag to 10m after Broede et al. (2012), eq.3
+!   z(0) is set to 0.01 according to UTCI definition
+    va = DBLE( v_hag *  log ( 10.0_wp / 0.01_wp ) / log ( hag / 0.01_wp ) )
+
+!-- Check if input values in range after Broede et al. (2012)
+    IF ( ( d_tmrt > 70.d0 ) .OR. ( d_tmrt < -30.d0 ) .OR.                      &
+       ( ehpa .GE. 50._wp ) ) THEN
+       utci = -999.0_wp
+       RETURN
+    ENDIF
+
+!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
+    IF ( ta > 50.d0 ) THEN
+       offset = ta - 50.d0
+       ta = 50.d0
+    ENDIF
+    IF ( ta < -50.d0 ) THEN
+       offset = ta + 50.d0
+       ta = -50.d0
+    ENDIF
+
+!-- For routine application. For wind speeds and relative
+!   humidity values below 0.5 m/s or 5%, respectively, the
+!   user is advised to use the lower bounds for the calculations.
+    IF ( va < 0.5d0 ) va  = 0.5d0
+    IF ( va > 17.0d0 ) va = 17.0d0
+
+!-- Calculate 6th order polynomial as approximation 
+    utci_d = ta +                                                              &
+       ( 6.07562052D-01 )   +                                                  &
+       ( -2.27712343D-02 ) * ta +                                              &
+       ( 8.06470249D-04 )  * ta * ta +                                         &
+       ( -1.54271372D-04 ) * ta * ta * ta +                                    &
+       ( -3.24651735D-06 ) * ta * ta * ta * ta +                               &
+       ( 7.32602852D-08 )  * ta * ta * ta * ta * ta +                          &
+       ( 1.35959073D-09 )  * ta * ta * ta * ta * ta * ta +                     &
+       ( -2.25836520D+00 ) * va +                                              &
+       ( 8.80326035D-02 )  * ta * va +                                         &
+       ( 2.16844454D-03 )  * ta * ta * va +                                    &
+       ( -1.53347087D-05 ) * ta * ta * ta * va +                               &
+       ( -5.72983704D-07 ) * ta * ta * ta * ta * va +                          &
+       ( -2.55090145D-09 ) * ta * ta * ta * ta * ta * va +                     &
+       ( -7.51269505D-01 ) * va * va +                                         &
+       ( -4.08350271D-03 ) * ta * va * va +                                    &
+       ( -5.21670675D-05 ) * ta * ta * va * va +                               &
+       ( 1.94544667D-06 )  * ta * ta * ta * va * va +                          &
+       ( 1.14099531D-08 )  * ta * ta * ta * ta * va * va +                     &
+       ( 1.58137256D-01 )  * va * va * va +                                    &
+       ( -6.57263143D-05 ) * ta * va * va * va +                               &
+       ( 2.22697524D-07 )  * ta * ta * va * va * va +                          &
+       ( -4.16117031D-08 ) * ta * ta * ta * va * va * va +                     &
+       ( -1.27762753D-02 ) * va * va * va * va +                               &
+       ( 9.66891875D-06 )  * ta * va * va * va * va +                          &
+       ( 2.52785852D-09 )  * ta * ta * va * va * va * va +                     &
+       ( 4.56306672D-04 )  * va * va * va * va * va +                          &
+       ( -1.74202546D-07 ) * ta * va * va * va * va * va +                     &
+       ( -5.91491269D-06 ) * va * va * va * va * va * va +                     &
+       ( 3.98374029D-01 )  * d_tmrt +                                          &
+       ( 1.83945314D-04 )  * ta * d_tmrt +                                     &
+       ( -1.73754510D-04 ) * ta * ta * d_tmrt +                                &
+       ( -7.60781159D-07 ) * ta * ta * ta * d_tmrt +                           &
+       ( 3.77830287D-08 )  * ta * ta * ta * ta * d_tmrt +                      &
+       ( 5.43079673D-10 )  * ta * ta * ta * ta * ta * d_tmrt +                 &
+       ( -2.00518269D-02 ) * va * d_tmrt +                                     &
+       ( 8.92859837D-04 )  * ta * va * d_tmrt +                                &
+       ( 3.45433048D-06 )  * ta * ta * va * d_tmrt +                           &
+       ( -3.77925774D-07 ) * ta * ta * ta * va * d_tmrt +                      &
+       ( -1.69699377D-09 ) * ta * ta * ta * ta * va * d_tmrt +                 &
+       ( 1.69992415D-04 )  * va * va * d_tmrt +                                &
+       ( -4.99204314D-05 ) * ta * va * va * d_tmrt +                           &
+       ( 2.47417178D-07 )  * ta * ta * va * va * d_tmrt +                      &
+       ( 1.07596466D-08 )  * ta * ta * ta * va * va * d_tmrt +                 &
+       ( 8.49242932D-05 )  * va * va * va * d_tmrt +                           &
+       ( 1.35191328D-06 )  * ta * va * va * va * d_tmrt +                      &
+       ( -6.21531254D-09 ) * ta * ta * va * va * va * d_tmrt +                 &
+       ( -4.99410301D-06 ) * va * va * va * va * d_tmrt +                      &
+       ( -1.89489258D-08 ) * ta * va * va * va * va * d_tmrt +                 &
+       ( 8.15300114D-08 )  * va * va * va * va * va * d_tmrt +                 &
+       ( 7.55043090D-04 )  * d_tmrt * d_tmrt +                                 &
+       ( -5.65095215D-05 ) * ta * d_tmrt * d_tmrt +                            &
+       ( -4.52166564D-07 ) * ta * ta * d_tmrt * d_tmrt +                       &
+       ( 2.46688878D-08 )  * ta * ta * ta * d_tmrt * d_tmrt +                  &
+       ( 2.42674348D-10 )  * ta * ta * ta * ta * d_tmrt * d_tmrt +             &
+       ( 1.54547250D-04 )  * va * d_tmrt * d_tmrt +                            &
+       ( 5.24110970D-06 )  * ta * va * d_tmrt * d_tmrt +                       &
+       ( -8.75874982D-08 ) * ta * ta * va * d_tmrt * d_tmrt +                  &
+       ( -1.50743064D-09 ) * ta * ta * ta * va * d_tmrt * d_tmrt +             &
+       ( -1.56236307D-05 ) * va * va * d_tmrt * d_tmrt +                       &
+       ( -1.33895614D-07 ) * ta * va * va * d_tmrt * d_tmrt +                  &
+       ( 2.49709824D-09 )  * ta * ta * va * va * d_tmrt * d_tmrt +             &
+       ( 6.51711721D-07 )  * va * va * va * d_tmrt * d_tmrt +                  &
+       ( 1.94960053D-09 )  * ta * va * va * va * d_tmrt * d_tmrt +             &
+       ( -1.00361113D-08 ) * va * va * va * va * d_tmrt * d_tmrt +             &
+       ( -1.21206673D-05 ) * d_tmrt * d_tmrt * d_tmrt +                        &
+       ( -2.18203660D-07 ) * ta * d_tmrt * d_tmrt * d_tmrt +                   &
+       ( 7.51269482D-09 )  * ta * ta * d_tmrt * d_tmrt * d_tmrt +              &
+       ( 9.79063848D-11 )  * ta * ta * ta * d_tmrt * d_tmrt * d_tmrt +         &
+       ( 1.25006734D-06 )  * va * d_tmrt * d_tmrt * d_tmrt +                   &
+       ( -1.81584736D-09 ) * ta * va * d_tmrt * d_tmrt * d_tmrt +              &
+       ( -3.52197671D-10 ) * ta * ta * va * d_tmrt * d_tmrt * d_tmrt +         &
+       ( -3.36514630D-08 ) * va * va * d_tmrt * d_tmrt * d_tmrt +              &
+       ( 1.35908359D-10 )  * ta * va * va * d_tmrt * d_tmrt * d_tmrt +         &
+       ( 4.17032620D-10 )  * va * va * va * d_tmrt * d_tmrt * d_tmrt +         &
+       ( -1.30369025D-09 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt +               &
+       ( 4.13908461D-10 )  * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt +          &
+       ( 9.22652254D-12 )  * ta * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
+       ( -5.08220384D-09 ) * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +          &
+       ( -2.24730961D-11 ) * ta * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
+       ( 1.17139133D-10 )  * va * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
+       ( 6.62154879D-10 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt +      &
+       ( 4.03863260D-13 )  * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + &
+       ( 1.95087203D-12 )  * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + &
+       ( -4.73602469D-12 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt *      &
+       d_tmrt +                                                                &
+       ( 5.12733497D+00 )  * pa +                                              &
+       ( -3.12788561D-01 ) * ta * pa +                                         &
+       ( -1.96701861D-02 ) * ta * ta * pa +                                    &
+       ( 9.99690870D-04 )  * ta * ta * ta * pa +                               &
+       ( 9.51738512D-06 )  * ta * ta * ta * ta * pa +                          &
+       ( -4.66426341D-07 ) * ta * ta * ta * ta * ta * pa +                     &
+       ( 5.48050612D-01 )  * va * pa +                                         &
+       ( -3.30552823D-03 ) * ta * va * pa +                                    &
+       ( -1.64119440D-03 ) * ta * ta * va * pa +                               &
+       ( -5.16670694D-06 ) * ta * ta * ta * va * pa +                          &
+       ( 9.52692432D-07 )  * ta * ta * ta * ta * va * pa +                     &
+       ( -4.29223622D-02 ) * va * va * pa +                                    &
+       ( 5.00845667D-03 )  * ta * va * va * pa +                               &
+       ( 1.00601257D-06 )  * ta * ta * va * va * pa +                          &
+       ( -1.81748644D-06 ) * ta * ta * ta * va * va * pa +                     &
+       ( -1.25813502D-03 ) * va * va * va * pa +                               &
+       ( -1.79330391D-04 ) * ta * va * va * va * pa +                          &
+       ( 2.34994441D-06 )  * ta * ta * va * va * va * pa +                     &
+       ( 1.29735808D-04 )  * va * va * va * va * pa +                          &
+       ( 1.29064870D-06 )  * ta * va * va * va * va * pa +                     &
+       ( -2.28558686D-06 ) * va * va * va * va * va * pa +                     &
+       ( -3.69476348D-02 ) * d_tmrt * pa +                                     &
+       ( 1.62325322D-03 )  * ta * d_tmrt * pa +                                &
+       ( -3.14279680D-05 ) * ta * ta * d_tmrt * pa +                           &
+       ( 2.59835559D-06 )  * ta * ta * ta * d_tmrt * pa +                      &
+       ( -4.77136523D-08 ) * ta * ta * ta * ta * d_tmrt * pa +                 &
+       ( 8.64203390D-03 )  * va * d_tmrt * pa +                                &
+       ( -6.87405181D-04 ) * ta * va * d_tmrt * pa +                           &
+       ( -9.13863872D-06 ) * ta * ta * va * d_tmrt * pa +                      &
+       ( 5.15916806D-07 )  * ta * ta * ta * va * d_tmrt * pa +                 &
+       ( -3.59217476D-05 ) * va * va * d_tmrt * pa +                           &
+       ( 3.28696511D-05 )  * ta * va * va * d_tmrt * pa +                      &
+       ( -7.10542454D-07 ) * ta * ta * va * va * d_tmrt * pa +                 &
+       ( -1.24382300D-05 ) * va * va * va * d_tmrt * pa +                      &
+       ( -7.38584400D-09 ) * ta * va * va * va * d_tmrt * pa +                 &
+       ( 2.20609296D-07 )  * va * va * va * va * d_tmrt * pa +                 &
+       ( -7.32469180D-04 ) * d_tmrt * d_tmrt * pa +                            &
+       ( -1.87381964D-05 ) * ta * d_tmrt * d_tmrt * pa +                       &
+       ( 4.80925239D-06 )  * ta * ta * d_tmrt * d_tmrt * pa +                  &
+       ( -8.75492040D-08 ) * ta * ta * ta * d_tmrt * d_tmrt * pa +             &
+       ( 2.77862930D-05 )  * va * d_tmrt * d_tmrt * pa +                       &
+       ( -5.06004592D-06 ) * ta * va * d_tmrt * d_tmrt * pa +                  &
+       ( 1.14325367D-07 )  * ta * ta * va * d_tmrt * d_tmrt * pa +             &
+       ( 2.53016723D-06 )  * va * va * d_tmrt * d_tmrt * pa +                  &
+       ( -1.72857035D-08 ) * ta * va * va * d_tmrt * d_tmrt * pa +             &
+       ( -3.95079398D-08 ) * va * va * va * d_tmrt * d_tmrt * pa +             &
+       ( -3.59413173D-07 ) * d_tmrt * d_tmrt * d_tmrt * pa +                   &
+       ( 7.04388046D-07 )  * ta * d_tmrt * d_tmrt * d_tmrt * pa +              &
+       ( -1.89309167D-08 ) * ta * ta * d_tmrt * d_tmrt * d_tmrt * pa +         &
+       ( -4.79768731D-07 ) * va * d_tmrt * d_tmrt * d_tmrt * pa +              &
+       ( 7.96079978D-09 )  * ta * va * d_tmrt * d_tmrt * d_tmrt * pa +         &
+       ( 1.62897058D-09 )  * va * va * d_tmrt * d_tmrt * d_tmrt * pa +         &
+       ( 3.94367674D-08 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +          &
+       ( -1.18566247D-09 ) * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +     &
+       ( 3.34678041D-10 )  * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +     &
+       ( -1.15606447D-10 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa + &
+       ( -2.80626406D+00 ) * pa * pa +                                         &
+       ( 5.48712484D-01 )  * ta * pa * pa +                                    &
+       ( -3.99428410D-03 ) * ta * ta * pa * pa +                               &
+       ( -9.54009191D-04 ) * ta * ta * ta * pa * pa +                          &
+       ( 1.93090978D-05 )  * ta * ta * ta * ta * pa * pa +                     &
+       ( -3.08806365D-01 ) * va * pa * pa +                                    &
+       ( 1.16952364D-02 )  * ta * va * pa * pa +                               &
+       ( 4.95271903D-04 )  * ta * ta * va * pa * pa +                          &
+       ( -1.90710882D-05 ) * ta * ta * ta * va * pa * pa +                     &
+       ( 2.10787756D-03 )  * va * va * pa * pa +                               &
+       ( -6.98445738D-04 ) * ta * va * va * pa * pa +                          &
+       ( 2.30109073D-05 )  * ta * ta * va * va * pa * pa +                     &
+       ( 4.17856590D-04 )  * va * va * va * pa * pa +                          &
+       ( -1.27043871D-05 ) * ta * va * va * va * pa * pa +                     &
+       ( -3.04620472D-06 ) * va * va * va * va * pa * pa +                     &
+       ( 5.14507424D-02 )  * d_tmrt * pa * pa +                                &
+       ( -4.32510997D-03 ) * ta * d_tmrt * pa * pa +                           &
+       ( 8.99281156D-05 )  * ta * ta * d_tmrt * pa * pa +                      &
+       ( -7.14663943D-07 ) * ta * ta * ta * d_tmrt * pa * pa +                 &
+       ( -2.66016305D-04 ) * va * d_tmrt * pa * pa +                           &
+       ( 2.63789586D-04 )  * ta * va * d_tmrt * pa * pa +                      &
+       ( -7.01199003D-06 ) * ta * ta * va * d_tmrt * pa * pa +                 &
+       ( -1.06823306D-04 ) * va * va * d_tmrt * pa * pa +                      &
+       ( 3.61341136D-06 )  * ta * va * va * d_tmrt * pa * pa +                 &
+       ( 2.29748967D-07 )  * va * va * va * d_tmrt * pa * pa +                 &
+       ( 3.04788893D-04 )  * d_tmrt * d_tmrt * pa * pa +                       &
+       ( -6.42070836D-05 ) * ta * d_tmrt * d_tmrt * pa * pa +                  &
+       ( 1.16257971D-06 )  * ta * ta * d_tmrt * d_tmrt * pa * pa +             &
+       ( 7.68023384D-06 )  * va * d_tmrt * d_tmrt * pa * pa +                  &
+       ( -5.47446896D-07 ) * ta * va * d_tmrt * d_tmrt * pa * pa +             &
+       ( -3.59937910D-08 ) * va * va * d_tmrt * d_tmrt * pa * pa +             &
+       ( -4.36497725D-06 ) * d_tmrt * d_tmrt * d_tmrt * pa * pa +              &
+       ( 1.68737969D-07 )  * ta * d_tmrt * d_tmrt * d_tmrt * pa * pa +         &
+       ( 2.67489271D-08 )  * va * d_tmrt * d_tmrt * d_tmrt * pa * pa +         &
+       ( 3.23926897D-09 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa * pa +     &
+       ( -3.53874123D-02 ) * pa * pa * pa +                                    &
+       ( -2.21201190D-01 ) * ta * pa * pa * pa +                               &
+       ( 1.55126038D-02 )  * ta * ta * pa * pa * pa +                          &
+       ( -2.63917279D-04 ) * ta * ta * ta * pa * pa * pa +                     &
+       ( 4.53433455D-02 )  * va * pa * pa * pa +                               &
+       ( -4.32943862D-03 ) * ta * va * pa * pa * pa +                          &
+       ( 1.45389826D-04 )  * ta * ta * va * pa * pa * pa +                     &
+       ( 2.17508610D-04 )  * va * va * pa * pa * pa +                          &
+       ( -6.66724702D-05 ) * ta * va * va * pa * pa * pa +                     &
+       ( 3.33217140D-05 )  * va * va * va * pa * pa * pa +                     &
+       ( -2.26921615D-03 ) * d_tmrt * pa * pa * pa +                           &
+       ( 3.80261982D-04 )  * ta * d_tmrt * pa * pa * pa +                      &
+       ( -5.45314314D-09 ) * ta * ta * d_tmrt * pa * pa * pa +                 &
+       ( -7.96355448D-04 ) * va * d_tmrt * pa * pa * pa +                      &
+       ( 2.53458034D-05 )  * ta * va * d_tmrt * pa * pa * pa +                 &
+       ( -6.31223658D-06 ) * va * va * d_tmrt * pa * pa * pa +                 &
+       ( 3.02122035D-04 )  * d_tmrt * d_tmrt * pa * pa * pa +                  &
+       ( -4.77403547D-06 ) * ta * d_tmrt * d_tmrt * pa * pa * pa +             &
+       ( 1.73825715D-06 )  * va * d_tmrt * d_tmrt * pa * pa * pa +             &
+       ( -4.09087898D-07 ) * d_tmrt * d_tmrt * d_tmrt * pa * pa * pa +         &
+       ( 6.14155345D-01 )  * pa * pa * pa * pa +                               &
+       ( -6.16755931D-02 ) * ta * pa * pa * pa * pa +                          &
+       ( 1.33374846D-03 )  * ta * ta * pa * pa * pa * pa +                     &
+       ( 3.55375387D-03 )  * va * pa * pa * pa * pa +                          &
+       ( -5.13027851D-04 ) * ta * va * pa * pa * pa * pa +                     &
+       ( 1.02449757D-04 )  * va * va * pa * pa * pa * pa +                     &
+       ( -1.48526421D-03 ) * d_tmrt * pa * pa * pa * pa +                      &
+       ( -4.11469183D-05 ) * ta * d_tmrt * pa * pa * pa * pa +                 &
+       ( -6.80434415D-06 ) * va * d_tmrt * pa * pa * pa * pa +                 &
+       ( -9.77675906D-06 ) * d_tmrt * d_tmrt * pa * pa * pa * pa +             &
+       ( 8.82773108D-02 )  * pa * pa * pa * pa * pa +                          &
+       ( -3.01859306D-03 ) * ta * pa * pa * pa * pa * pa +                     &
+       ( 1.04452989D-03 )  * va * pa * pa * pa * pa * pa +                     &
+       ( 2.47090539D-04 )  * d_tmrt * pa * pa * pa * pa * pa +                 &
+       ( 1.48348065D-03 )  * pa * pa * pa * pa * pa * pa 
+
+!-- Cast result to working precision
+    utci = REAL( utci_d + offset, wp )
+
+ END SUBROUTINE calculate_utci_static
+
+
+ END MODULE biometeorology_utci_mod
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 3445)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add biometeorology
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Renamed output variables
 ! Add checks for surface data output
@@ -704,4 +707,7 @@
     USE basic_constants_and_equations_mod
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_check_data_output, biom_check_parameters
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model, bcm_check_parameters, bcm_check_data_output,  &
@@ -769,7 +775,4 @@
         ONLY:  radiation, radiation_check_data_output,                         &
                radiation_check_data_output_pr, radiation_check_parameters
-
-    USE biometeorology_mod,                                                    &
-        ONLY:  biometeorology_check_data_output
 
     USE spectra_mod,                                                           &
@@ -1488,4 +1491,5 @@
 
 !-- Check the module settings
+    IF ( biometeorology )       CALL biom_check_parameters
     IF ( bulk_cloud_model )     CALL bcm_check_parameters
     IF ( gust_module_enabled )  CALL gust_check_parameters
@@ -3182,10 +3186,10 @@
              ENDIF
 
-             IF ( unit == 'illegal' )  THEN
-                CALL biometeorology_check_data_output( var, unit, i, ilen, k )
-             ENDIF
-
              IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
                 CALL gust_check_data_output ( var, unit )
+             ENDIF
+
+             IF ( unit == 'illegal'  .AND.  biometeorology )  THEN
+                CALL biom_check_data_output( var, unit )
              ENDIF
 
Index: palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_2d.f90	(revision 3445)
+++ palm/trunk/SOURCE/data_output_2d.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add biometeorology
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Renamed output variables
 ! 
@@ -256,4 +259,7 @@
         ONLY:  c_p, lv_d_cp, l_v
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_data_output_2d
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model, bcm_data_output_2d
@@ -263,6 +269,6 @@
 
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, data_output_2d_on_each_pe, data_output_xy,       &
-               data_output_xz, data_output_yz, do2d,                           &
+        ONLY:  air_chemistry, biometeorology, data_output_2d_on_each_pe,       &
+               data_output_xy, data_output_xz, data_output_yz, do2d,           &
                do2d_xy_last_time, do2d_xy_time_count,                          &
                do2d_xz_last_time, do2d_xz_time_count,                          &
@@ -1308,4 +1314,11 @@
                    CALL gust_data_output_2d( av, do2d(av,if), found, grid,     &
                                              local_pf, two_d, nzb_do, nzt_do )
+                ENDIF
+
+                IF ( .NOT. found  .AND.  biometeorology                        &
+                                  .AND.  mode == 'xy' )  THEN
+                   CALL biom_data_output_2d( av, do2d(av,if), found, grid,     &
+                                             local_pf, two_d, nzb_do, nzt_do,  &
+                                             fill_value )
                 ENDIF
 
Index: palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- palm/trunk/SOURCE/data_output_3d.f90	(revision 3445)
+++ palm/trunk/SOURCE/data_output_3d.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Adjustment of biometeorology calls
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Renamed output variables
 ! 
@@ -229,4 +232,7 @@
         ONLY:  lv_d_cp
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_data_output_3d
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model, bcm_data_output_3d
@@ -236,5 +242,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, do3d, do3d_no, do3d_time_count,                  &
+        ONLY:  air_chemistry, biometeorology, do3d, do3d_no, do3d_time_count,  &
                io_blocks, io_group, land_surface, message_string,              &
                ntdim_3d, nz_do3d,  ocean_mode, plant_canopy,                   &
@@ -285,7 +291,4 @@
     USE radiation_model_mod,                                                   &
         ONLY:  nzub, nzut, radiation, radiation_data_output_3d
-
-    USE biometeorology_mod,                                                    &
-        ONLY: biometeorology_data_output_3d
 
     USE turbulence_closure_mod,                                                &
@@ -767,16 +770,13 @@
              ENDIF
 
-!
-!--          Biometeorology quantity
-             IF ( .NOT. found  .AND.  radiation )  THEN
-                CALL biometeorology_data_output_3d( av, do3d(av,if), found,    &
-                                               local_pf, nzb_do, nzt_do )
-                resorted = .TRUE.
-             ENDIF
-
              IF ( .NOT. found )  THEN
                 CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
                                          nzb_do, nzt_do )
                 resorted = .TRUE.
+             ENDIF
+
+             IF ( .NOT. found  .AND.  biometeorology )  THEN
+                 CALL biom_data_output_3d( av, do3d(av,if), found, local_pf,   &
+                                           nzb_do, nzt_do )
              ENDIF
 
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 3445)
+++ palm/trunk/SOURCE/header.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add biometeorology
+! 
+! 3355 2018-10-16 14:03:34Z knoop
 ! Header output concerning offline nesting
 ! 
@@ -238,4 +241,7 @@
 ! 1551 2015-03-03 14:18:16Z maronga
 ! Added informal output for land surface model and radiation model. Removed typo.
+! 
+! 1496 2014-12-02 17:25:50Z maronga
+! Renamed: "radiation -> "cloud_top_radiation"
 ! 
 ! 1484 2014-10-21 10:53:05Z kanani
@@ -395,4 +401,7 @@
         ONLY:  g, kappa, l_v, r_d
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_header
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model, bcm_header
@@ -420,5 +429,5 @@
     USE gust_mod,                                                              &
         ONLY: gust_header, gust_module_enabled
-        
+
     USE indices,                                                               &
         ONLY:  mg_loc_ind, nnx, nny, nnz, nx, ny, nxl_mg, nxr_mg, nyn_mg,      &
@@ -438,5 +447,5 @@
     USE netcdf_interface,                                                      &
         ONLY:  netcdf_data_format, netcdf_data_format_string, netcdf_deflate
-        
+
     USE nesting_offl_mod,                                                      &
         ONLY:  nesting_offl_header
@@ -906,5 +915,5 @@
 !-- Large scale forcing and nudging
     IF ( large_scale_forcing )  CALL lsf_nudging_header( io )
-    
+
 !
 !-- Offline nesting 
@@ -1133,4 +1142,5 @@
 !
 !-- Header information from other modules
+    IF ( biometeorology      )  CALL biom_header ( io )
     IF ( gust_module_enabled )  CALL gust_header( io )
     IF ( land_surface        )  CALL lsm_header( io )
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 3445)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add biometeorology
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Initialize surface data output
 ! 
@@ -533,4 +536,7 @@
                ideal_gas_law_rho, ideal_gas_law_rho_pt, barometric_formula
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_init, biom_init_arrays
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model, bcm_init, bcm_init_arrays
@@ -552,5 +558,5 @@
     USE gust_mod,                                                              &
         ONLY:  gust_init, gust_init_arrays, gust_module_enabled
-    
+
     USE indices
     
@@ -578,5 +584,5 @@
         ONLY:  chem_emis, chem_emis_att, init_3d,                              &
                netcdf_data_input_init_3d, netcdf_data_input_interpolate
-        
+
     USE nesting_offl_mod,                                                      &
         ONLY:  nesting_offl_init
@@ -628,5 +634,5 @@
         ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
                 surf_usm_h, get_topography_top_index_ji, vertical_surfaces_exist
-                
+   
     USE surface_output_mod,                                                    &
         ONLY:  surface_output_init
@@ -1050,4 +1056,5 @@
 !
 !-- Allocate arrays for other modules
+    IF ( biometeorology      )  CALL biom_init_arrays
     IF ( bulk_cloud_model    )  CALL bcm_init_arrays
     IF ( gust_module_enabled )  CALL gust_init_arrays
@@ -2413,6 +2420,14 @@
 
 !
+!-- If required initialize biometeorology module
+    IF ( biometeorology )  THEN
+        CALL location_message( 'initializing biometeorology module', .FALSE. )
+        CALL biom_init
+        CALL location_message( 'finished', .TRUE. )
+    ENDIF
+
+!
 !-- Initialize the ws-scheme.    
-    IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init        
+    IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init
 
 !
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 3445)
+++ palm/trunk/SOURCE/modules.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add biometeorology
+! 
+! 3435 2018-10-26 18:25:44Z gronemeier
 ! +mask_k_over_surface, mask_surface
 ! 
@@ -425,4 +428,7 @@
 ! Increased pr_palm to 120. Increased length of dots_unit and dots_label to 13
 ! digits. Increased length of domask, do2d, and do3d to 20 digits.
+! 
+! 1496 2014-12-02 17:25:50Z maronga
+! Renamed "radiation" -> "cloud_top_radiation"
 ! 
 ! 1484 2014-10-21 10:53:05Z kanani
@@ -1305,4 +1311,5 @@
     LOGICAL ::  bc_radiation_r = .FALSE.                         !< radiation boundary condition for outflow at right domain boundary
     LOGICAL ::  bc_radiation_s = .FALSE.                         !< radiation boundary condition for outflow at south domain boundary
+    LOGICAL ::  biometeorology = .FALSE.                         !< biometeorology module switch
     LOGICAL ::  calc_soil_moisture_during_spinup = .FALSE.       !< namelist parameter
     LOGICAL ::  call_psolver_at_all_substeps = .TRUE.            !< namelist parameter
Index: palm/trunk/SOURCE/multi_agent_system_mod.f90
===================================================================
--- palm/trunk/SOURCE/multi_agent_system_mod.f90	(revision 3445)
+++ palm/trunk/SOURCE/multi_agent_system_mod.f90	(revision 3448)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Adjustment of biometeorology calls,
+! implement some agent biometeorology
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -84,5 +88,6 @@
 
     USE control_parameters,                                                    &
-        ONLY:  dt_3d, message_string, time_since_reference_point, dt_write_agent_data
+        ONLY:  biometeorology, dt_3d, dt_write_agent_data, message_string,     &
+               time_since_reference_point
 
     USE cpulog,                                                                &
@@ -206,4 +211,11 @@
         REAL(wp)     ::  clo                  !< clothing index
         REAL(wp)     ::  energy_storage       !< energy stored by agent
+        REAL(wp)     ::  clothing_temp        !< energy stored by agent
+        REAL(wp)     ::  actlev               !< metabolic + work energy of the person
+        REAL(wp)     ::  age_years            !< physical age of the person
+        REAL(wp)     ::  weight               !< total weight of the person (kg)
+        REAL(wp)     ::  height               !< height of the person (m)
+        REAL(wp)     ::  work                 !< workload of the agent (W)
+        INTEGER(iwp) ::  sex                  !< agents gender: 1 = male, 2 = female
         REAL(wp)     ::  force_x              !< force term x-direction
         REAL(wp)     ::  force_y              !< force term y-direction
@@ -218,5 +230,5 @@
         REAL(wp)     ::  speed_x              !< speed of agent in x
         REAL(wp)     ::  speed_y              !< speed of agent in y
-        REAL(wp)     ::  thermal_index        !< the dynamic thermal index
+        REAL(wp)     ::  ipt                  !< instationary thermal index iPT (Â°C)
         REAL(wp)     ::  windspeed            !< absolute value of windspeed at agent position
         REAL(wp)     ::  x                    !< x-position
@@ -329,6 +341,7 @@
  SUBROUTINE multi_agent_system
 
-!     USE htcm_mod,                                                              &
-!         ONLY:  htcm_dynamic
+    USE biometeorology_mod,                                                   &
+        ONLY:  biom_calc_ipt, biom_determine_input_at
+
 
     IMPLICIT NONE
@@ -341,4 +354,12 @@
     INTEGER(iwp)       ::  js                 !< counter
     INTEGER(iwp), SAVE ::  mas_count = 0      !< counts the mas-calls
+    INTEGER(iwp)                :: a     !< agent iterator
+    !-- local meteorological conditions
+    REAL(wp)                    :: tmrt  !< mean radiant temperature        (Â°C)
+    REAL(wp)                    :: ta    !< air temperature                 (Â°C)
+    REAL(wp)                    :: vp    !< vapour pressure                 (hPa)
+    REAL(wp)                    :: v     !< wind speed    (local level)     (m/s)
+    REAL(wp)                    :: pair  !< air pressure                    (hPa)
+
 
     LOGICAL       ::  first_loop_stride   !< flag for first loop stride of agent sub-timesteps
@@ -502,20 +523,45 @@
        deleted_agents = 0
 !
-!--    to be included here: call of human thermal comfort mod (and UV exposure)
-!        DO  i = nxl, nxr
-!           DO  j = nys, nyn
-! 
-!              number_of_agents = agt_count(j,i)
-! !
-! !--          If grid cell gets empty, cycle
-!              IF ( number_of_agents <= 0 ) CYCLE
-! 
-!              agents => grid_agents(j,i)%agents(1:number_of_agents)
-! !
-! !--          Evaluation of social forces
-!              CALL htcm_dynamic(i,j)
-! 
-!           ENDDO
-!        ENDDO
+       IF ( biometeorology )  THEN
+!
+!--       Call of human thermal comfort mod (and UV exposure)
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+
+                number_of_agents = agt_count(j,i)
+!
+!--             If grid cell gets empty, cycle
+                IF ( number_of_agents <= 0 )  CYCLE
+
+                agents => grid_agents(j,i)%agents(1:number_of_agents)
+!
+!--             Evaluation of social forces
+!                CALL biom_dynamic(i,j)
+!
+!--             Determine local meteorological conditions
+                CALL biom_determine_input_at ( .FALSE., i, j, ta, vp, v,      &
+                                               pair, tmrt )
+
+                DO  a = 1, number_of_agents
+!
+!--                Calculate instationary thermal indices based on local tmrt
+
+                   CALL biom_calc_ipt ( ta, vp, v, pair, tmrt,                &
+                                        agents(a)%dt_sum,                     &
+                                        agents(a)%energy_storage,             &
+                                        agents(a)%clothing_temp,              &
+                                        agents(a)%clo,                        &
+                                        agents(a)%actlev,                     &
+                                        agents(a)%age_years,                  &
+                                        agents(a)%weight,                     &
+                                        agents(a)%height,                     &
+                                        agents(a)%work,                       &
+                                        agents(a)%sex,                        &
+                                        agents(a)%ipt )
+                END DO
+
+             ENDDO
+          ENDDO
+       ENDIF
 
        IF ( dt_3d_reached_mas )  EXIT
@@ -799,6 +845,14 @@
                          tmp_agent%age_m         = 0.0_wp
                          tmp_agent%dt_sum        = 0.0_wp
-                         tmp_agent%clo           = 99999.0_wp
+                         tmp_agent%clo           = -999.0_wp
                          tmp_agent%energy_storage= 0.0_wp
+                         tmp_agent%ipt           = 99999.0_wp
+                         tmp_agent%clothing_temp = -999._wp      !< energy stored by agent (W)
+                         tmp_agent%actlev        = 134.6862_wp   !< metabolic + work energy of the person
+                         tmp_agent%age_years     = 35._wp        !< physical age of the person
+                         tmp_agent%weight        = 75._wp        !< total weight of the person (kg)
+                         tmp_agent%height        = 1.75_wp       !< height of the person (m)
+                         tmp_agent%work          = 134.6862_wp   !< workload of the agent (W)
+                         tmp_agent%sex           = 1             !< agents gender: 1 = male, 2 = female
                          tmp_agent%force_x       = 0.0_wp
                          tmp_agent%force_y       = 0.0_wp
@@ -812,5 +866,4 @@
                          tmp_agent%speed_x       = 0.0_wp
                          tmp_agent%speed_y       = 0.0_wp
-                         tmp_agent%thermal_index = 99999.0_wp
                          tmp_agent%x             = pos_x
                          tmp_agent%y             = pos_y
@@ -3185,5 +3238,5 @@
           zero_agent%speed_x       = 0.0_wp
           zero_agent%speed_y       = 0.0_wp
-          zero_agent%thermal_index = 0.0_wp
+          zero_agent%ipt           = 0.0_wp
           zero_agent%x             = 0.0_wp
           zero_agent%y             = 0.0_wp
Index: palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3445)
+++ palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Adjustment of biometeorology calls
+! 
+! 3435 2018-10-26 18:25:44Z gronemeier
 ! Bugfix: corrected order of calls to define_netcdf_grid for masked output
 ! Add vertical dimensions to masked output in case of terrain-following output
@@ -312,5 +315,6 @@
 
     USE control_parameters,                                                    &
-        ONLY:  fl_max, max_masks, multi_agent_system_end,                      &
+        ONLY:  biometeorology, fl_max,                                         &
+               max_masks, multi_agent_system_end,                              &
                multi_agent_system_start, var_fl_max, varnamelength
     USE kinds
@@ -326,5 +330,5 @@
           (/ 'ag_id           ', 'ag_x            ', 'ag_y            ',       &
              'ag_wind         ', 'ag_temp         ', 'ag_group        ',       &
-             'PM10            ', 'PM25            ', 'ag_therm_comf   ',       &
+             'PM10            ', 'PM25            ', 'ag_iPT          ',       &
              'ag_uv           ', 'not_used        ', 'not_used        ',       &
              'not_used        ' /)
@@ -513,5 +517,5 @@
                    id_var_zu_mask, id_var_zw_mask, &
                    id_var_zusi_mask, id_var_zwwi_mask
-                   
+
     INTEGER(iwp), DIMENSION(1:max_masks,0:1,0:2) ::  id_var_eutm_mask, &
                                                      id_var_nutm_mask
@@ -576,4 +580,7 @@
     USE arrays_3d,                                                             &
         ONLY:  zu, zw
+
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_define_netcdf_grid
 
     USE chemistry_model_mod,                                                   &
@@ -635,7 +642,4 @@
     USE radiation_model_mod,                                                   &
         ONLY: radiation, radiation_define_netcdf_grid
-
-    USE biometeorology_mod,                                                    &
-        ONLY: biometeorology_define_netcdf_grid
 
     USE spectra_mod,                                                           &
@@ -1132,11 +1136,4 @@
                    ENDIF
 !
-!--                Check for biometeorology quantities
-                   IF ( .NOT. found  .AND.  radiation )  THEN
-                      CALL biometeorology_define_netcdf_grid( domask(mid,av,i),&
-                                                         found, grid_x, grid_y,&
-                                                         grid_z )
-                   ENDIF
-!
 !--                Now check for user-defined quantities
                    IF ( .NOT. found )  THEN
@@ -1872,16 +1869,15 @@
                    ENDIF
 
-!
-!--                Check for biometeorology quantities
-                   IF ( .NOT. found  .AND.  radiation )  THEN
-                      CALL biometeorology_define_netcdf_grid( do3d(av,i), found,&
-                                                         grid_x, grid_y,       &
-                                                         grid_z )
-                   ENDIF
-
 !--                Check for gust module quantities
                    IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
                       CALL gust_define_netcdf_grid( do3d(av,i), found, grid_x, &
                                                     grid_y, grid_z )
+                   ENDIF
+
+!
+!--                Check for biometeorology quantities
+                   IF ( .NOT. found  .AND.  biometeorology )  THEN
+                      CALL biom_define_netcdf_grid( do3d(av,i), found,         &
+                                                    grid_x, grid_y, grid_z )
                    ENDIF
 
@@ -2469,4 +2465,12 @@
 
           ENDDO
+!
+!--       Define vars for biometeorology
+          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,                &
+                                  id_dim_time_agt /), agt_var_names(9),        &
+                                  nc_precision(8), id_var_agt(9),              &
+                                  TRIM( agt_var_units(9) ),                    &
+                                  TRIM( agt_var_names(9) ), 339, 340, 341 )
+
 !
 !--       Leave netCDF define mode
@@ -2827,12 +2831,4 @@
 
 !
-!--                      Check for biometeorology quantities
-                         IF ( .NOT. found  .AND.  radiation )  THEN
-                            CALL biometeorology_define_netcdf_grid( do2d(av,i),&
-                                                         found, grid_x, grid_y,&
-                                                         grid_z )
-                         ENDIF
-
-!
 !--                      Check for gust module quantities
                          IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
@@ -2841,5 +2837,11 @@
                                                           grid_z )
                          ENDIF
-
+!
+!--                      Check for human thermal comfort quantities
+                         IF ( .NOT. found  .AND.  biometeorology )  THEN
+                            CALL biom_define_netcdf_grid( do2d( av, i), found, &
+                                                          grid_x, grid_y,      &
+                                                          grid_z )
+                         ENDIF
 !
 !--                      Check for chemistry quantities
@@ -2894,4 +2896,6 @@
                    ELSEIF ( grid_z == 'zs' )  THEN
                       id_z = id_dim_zs_xy(av)
+                   ELSEIF ( grid_z == 'zu1' )  THEN
+                      id_z = id_dim_zu1_xy(av)
                    ENDIF
 
@@ -3188,5 +3192,5 @@
                                           count = (/ nx+1, ny+1 /) )
                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
-                
+
                 ENDDO
                 DEALLOCATE( netcdf_data_2d )
@@ -3724,13 +3728,4 @@
 
 !
-!--                   Check for biometeorology quantities
-                      IF ( .NOT. found  .AND.  radiation )  THEN
-                         CALL biometeorology_define_netcdf_grid( do2d(av,i),   &
-                                                            found,             &
-                                                            grid_x, grid_y,    &
-                                                            grid_z )
-                      ENDIF
-
-!
 !--                   Check for gust module quantities
                       IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
@@ -4582,13 +4577,4 @@
 
 !
-!--                   Check for biometeorology quantities
-                      IF ( .NOT. found  .AND.  radiation )  THEN
-                         CALL biometeorology_define_netcdf_grid( do2d(av,i),   &
-                                                            found,             &
-                                                            grid_x, grid_y,    &
-                                                            grid_z )
-                      ENDIF
-
-!
 !--                   Check for gust module quantities
                       IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 3445)
+++ palm/trunk/SOURCE/parin.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add biometeorology
+! 
+! 3435 2018-10-26 18:25:44Z gronemeier
 ! Add mask_k_over_surface
 ! 
@@ -304,4 +307,7 @@
 ! 1560 2015-03-06 10:48:54Z keck
 ! +recycling_yshift
+! 
+! 1496 2014-12-02 17:25:50Z maronga
+! Renamed: "radiation -> "cloud_top_radiation"
 ! 
 ! 1484 2014-10-21 10:53:05Z kanani
@@ -442,4 +448,7 @@
                ug, u_init, v_init, vg
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_parin
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bcm_parin
@@ -494,5 +503,5 @@
 
     USE pegrid
-               
+
     USE plant_canopy_model_mod,                                                &
          ONLY: pcm_parin
@@ -702,5 +711,5 @@
              dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_max, dt_restart,              &
              dt_run_control,end_time, force_print_header, mask_k_over_surface, &
-             mask_scale_x,        &
+             mask_scale_x,                                                     &
              mask_scale_y, mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop,  &
              mask_y_loop, mask_z_loop, netcdf_data_format, netcdf_deflate,     &
@@ -724,5 +733,5 @@
              dt_do2d_xz, dt_do2d_yz, dt_do3d, dt_max, dt_restart,              &
              dt_run_control,end_time, force_print_header, mask_k_over_surface, &
-             mask_scale_x,        &
+             mask_scale_x,                                                     &
              mask_scale_y, mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop,  &
              mask_y_loop, mask_z_loop, netcdf_data_format, netcdf_deflate,     &
@@ -884,11 +893,7 @@
 !
 !--       Check for module namelists and read them
+          CALL biom_parin
           CALL lsm_parin
-!
-!--       Check if bulk cloud model is used and read namelist if required
-          CALL bcm_parin 
-!
-!--       Check if surface output is enabled and read 
-!--       &surface_output_parameters if required          
+          CALL bcm_parin
           CALL surface_output_parin
           CALL usm_parin
Index: palm/trunk/SOURCE/sum_up_3d_data.f90
===================================================================
--- palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 3445)
+++ palm/trunk/SOURCE/sum_up_3d_data.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Adjustment of biometeorology calls
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Renamed output variables
 ! 
@@ -236,4 +239,7 @@
         ONLY:  c_p, lv_d_cp, l_v
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_3d_data_averaging
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY:  bulk_cloud_model, bcm_3d_data_averaging
@@ -243,7 +249,7 @@
 
     USE control_parameters,                                                    &
-        ONLY:  air_chemistry, average_count_3d, doav, doav_n, land_surface,    &
-               ocean_mode, rho_surface, urban_surface, uv_exposure,            &
-               varnamelength
+        ONLY:  air_chemistry, average_count_3d, biometeorology, doav, doav_n,  &
+               land_surface, ocean_mode, rho_surface, urban_surface,           &
+               uv_exposure, varnamelength
 
     USE cpulog,                                                                &
@@ -269,7 +275,4 @@
     USE radiation_model_mod,                                                   &
         ONLY:  radiation, radiation_3d_data_averaging
-
-    USE biometeorology_mod,                                                    &
-        ONLY:  biometeorology_3d_data_averaging
 
     USE surface_mod,                                                           &
@@ -517,10 +520,10 @@
                 ENDIF
 
-                IF ( radiation  .AND.  trimvar(1:4) == 'bio_' )  THEN
-                   CALL biometeorology_3d_data_averaging( 'allocate', doav(ii) )
-                ENDIF
-
                 IF ( gust_module_enabled )  THEN
                    CALL gust_3d_data_averaging( 'allocate', doav(ii) )
+                ENDIF
+
+                IF ( biometeorology  .AND.  trimvar(1:5) == 'biom_')  THEN
+                   CALL biom_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
@@ -1165,10 +1168,10 @@
              ENDIF
 
-             IF ( radiation   .AND.  trimvar(1:4) == 'bio_')  THEN
-                CALL biometeorology_3d_data_averaging( 'sum', doav(ii) )
-             ENDIF
-
              IF ( gust_module_enabled )  THEN
                 CALL gust_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
+             IF ( biometeorology  .AND.  trimvar(1:5) == 'biom_' )  THEN
+                   CALL biom_3d_data_averaging( 'sum', doav(ii) )
              ENDIF
 
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 3445)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add biometeorology
+! 
+! 3421 2018-10-24 18:39:32Z gronemeier
 ! Surface data output
 ! 
@@ -393,4 +396,7 @@
                v_p, w, w_p
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biom_calculate_static_grid, time_biom_results
+
     USE bulk_cloud_model_mod,                                                  &
         ONLY: bulk_cloud_model, calc_liquid_water_content,                     &
@@ -413,5 +419,5 @@
         ONLY:  advected_distance_x, advected_distance_y, air_chemistry,        &
                average_count_3d, averaging_interval, averaging_interval_pr,    &
-               bc_lr_cyc, bc_ns_cyc, bc_pt_t_val, bc_q_t_val,                  &
+               bc_lr_cyc, bc_ns_cyc, bc_pt_t_val, bc_q_t_val, biometeorology,  &
                call_psolver_at_all_substeps,  child_domain, cloud_droplets,    &
                constant_flux_layer, constant_heatflux,                         &
@@ -479,5 +485,5 @@
     USE multi_agent_system_mod,                                                &
         ONLY:  agents_active, multi_agent_system
-        
+
     USE nesting_offl_mod,                                                      &
         ONLY:  nesting_offl_bc, nesting_offl_mass_conservation
@@ -523,5 +529,5 @@
     USE surface_mod,                                                           &
         ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
-        
+
     USE surface_output_mod,                                                    &
         ONLY:  average_count_surf, averaging_interval_surf, dt_dosurf,         &
@@ -1285,4 +1291,11 @@
 
 !
+!--    Biometeorology calculation of stationary thermal indices
+       IF ( biometeorology  .AND.  time_do3d >= dt_do3d )  THEN
+          CALL biom_calculate_static_grid ( .FALSE. )
+          time_biom_results = time_since_reference_point
+       ENDIF
+
+!
 !--    Execute the gust module actions
        IF ( gust_module_enabled )  THEN
@@ -1548,5 +1561,5 @@
 
     ENDDO   ! time loop
-    
+
 !
 !-- Vertical nesting: Deallocate variables initialized for vertical nesting    
Index: palm/trunk/SOURCE/uv_exposure_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/uv_exposure_model_mod.f90	(revision 3445)
+++ palm/trunk/SOURCE/uv_exposure_model_mod.f90	(revision 3448)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Temporary rename of namelist until uv model moves to biometeorology module
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -454,5 +457,5 @@
     CHARACTER (LEN=80) ::  line  !< dummy string for current line in parameter file 
     
-    NAMELIST /biometeorology_parameters/  clothing
+    NAMELIST /biometeorology_parameters_uv/  clothing
     
     line = ' '
@@ -462,5 +465,5 @@
     REWIND ( 11 )
     line = ' '
-    DO WHILE ( INDEX( line, '&biometeorology_parameters' ) == 0 )
+    DO WHILE ( INDEX( line, '&biometeorology_parameters_uv' ) == 0 )
        READ ( 11, '(A)', END=20 )  line
     ENDDO
@@ -469,5 +472,5 @@
 !
 !-- Read user-defined namelist
-    READ ( 11, biometeorology_parameters, ERR = 10, END = 20 )
+    READ ( 11, biometeorology_parameters_uv, ERR = 10, END = 20 )
 
 !
@@ -477,5 +480,5 @@
  10 BACKSPACE( 11 )
     READ( 11 , '(A)') line
-    CALL parin_fail_message( 'biometeorology_parameters', line )
+    CALL parin_fail_message( 'biometeorology_parameters_uv', line )
 
 
