Changeset 3525 for palm

Timestamp:

Nov 14, 2018 4:06:14 PM (5 years ago)

Author:

kanani

Message:

Changes related to clean-up of biometeorology (by dom_dwd_user)

Location:

palm/trunk/SOURCE

Files:

• Makefile (modified) (3 diffs)
• average_3d_data.f90 (modified) (3 diffs)
• biometeorology_ipt_mod.f90 (deleted)
• biometeorology_mod.f90 (modified) (60 diffs)
• biometeorology_pet_mod.f90 (deleted)
• biometeorology_pt_mod.f90 (deleted)
• biometeorology_utci_mod.f90 (deleted)
• check_parameters.f90 (modified) (4 diffs)
• data_output_2d.f90 (modified) (3 diffs)
• data_output_3d.f90 (modified) (3 diffs)
• header.f90 (modified) (3 diffs)
• init_3d_model.f90 (modified) (4 diffs)
• multi_agent_system_mod.f90 (modified) (5 diffs)
• netcdf_interface_mod.f90 (modified) (4 diffs)
• parin.f90 (modified) (3 diffs)
• sum_up_3d_data.f90 (modified) (4 diffs)
• time_integration.f90 (modified) (3 diffs)

palm/trunk/SOURCE/Makefile

@@ -25 +25 @@
 # -----------------
 # $Id$
+# Changes related to clean-up of biometeorology (dom_dwd_user)
+# 
+# 3522 2018-11-13 12:14:36Z suehring
 # Dependencies for virtual measurement module added
 # 
@@ -528 +531 @@
         basic_constants_and_equations_mod.f90 \
         biometeorology_mod.f90 \
-        biometeorology_ipt_mod.f90 \
-        biometeorology_pt_mod.f90 \
-        biometeorology_pet_mod.f90 \
-        biometeorology_utci_mod.f90 \
         boundary_conds.f90 \
         buoyancy.f90 \
@@ -791 +790 @@
 basic_constants_and_equations_mod.o: \
         mod_kinds.o
-biometeorology_ipt_mod.o: \
-        biometeorology_pt_mod.o \
-        mod_kinds.o
 biometeorology_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
+        radiation_model_mod.o
 boundary_conds.o: \
         basic_constants_and_equations_mod.o \

palm/trunk/SOURCE/average_3d_data.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
-! Implementation of a new aerosol module salsa. 
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3467 2018-10-30 19:05:21Z suehring
+! Implementation of a new aerosol module salsa.
 ! 
 ! 3448 2018-10-29 18:14:31Z kanani
@@ -160 +163 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_3d_data_averaging
+        ONLY:  bio_3d_data_averaging
 
     USE bulk_cloud_model_mod,                                                  &
@@ -580 +583 @@
              ENDIF             
 
-             IF ( biometeorology )  THEN
-                CALL biom_3d_data_averaging( 'average', doav(ii) )
+             IF ( biometeorology .AND. trimvar(1:4) == 'bio_' )  THEN
+                CALL bio_3d_data_averaging( 'average', doav(ii) )
              ENDIF
 

palm/trunk/SOURCE/biometeorology_mod.f90

@@ -27 +27 @@
 ! -----------------
 ! $Id$
+! Clean up, renaming from "biom" to "bio", summary of thermal index calculation
+! into one module (dom_dwd_user)
+! 
+! 3479 2018-11-01 16:00:30Z gronemeier
 ! - reworked check for output quantities
 ! - assign new palm-error numbers
@@ -72 +76 @@
 
     USE basic_constants_and_equations_mod,                                     &
-       ONLY:  degc_to_k, magnus, sigma_sb
-
-    USE biometeorology_ipt_mod
-
-    USE biometeorology_pet_mod
-
-    USE biometeorology_pt_mod,                                                 &
-       ONLY: calculate_pt_static
-
-    USE biometeorology_utci_mod
+       ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb
 
     USE control_parameters,                                                    &
@@ -92 +87 @@
 
     USE indices,                                                               &
-       ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr
+       ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,    &
+              nysg, nyng
 
     USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
@@ -112 +108 @@
 !-- 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 ::  perct      !< 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)
+!-- Grids for averaged thermal indices
+    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< 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
+
+
+    INTEGER( iwp ) ::  bio_cell_level     !< cell level biom calculates for
+    REAL ( wp )    ::  bio_output_height  !< height output is calculated in m
+    REAL ( wp )    ::  time_bio_results   !< the time, the last set of biom results have been calculated for
     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 ::  mrt_from_rtm   = .TRUE.   !< switch: use mrt calculated by RTM for calculation of thermal indices
-    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
+    LOGICAL ::  aver_perct = .FALSE.  !< switch: do perct 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 ::  aver_mrt   = .FALSE.  !< switch: do mrt averaging in this module?
+    LOGICAL ::  average_trigger_perct = .FALSE.  !< update averaged input on call to bio_perct?
+    LOGICAL ::  average_trigger_utci  = .FALSE.  !< update averaged input on call to bio_utci?
+    LOGICAL ::  average_trigger_pet   = .FALSE.  !< update averaged input on call to bio_pet?
+
+    LOGICAL ::  bio_perct     = .TRUE.   !< Turn index PT (instant. input) on or off
+    LOGICAL ::  bio_perct_av  = .TRUE.   !< Turn index PT (averaged input) on or off
+    LOGICAL ::  bio_pet       = .TRUE.   !< Turn index PET (instant. input) on or off
+    LOGICAL ::  bio_pet_av    = .TRUE.   !< Turn index PET (averaged input) on or off
+    LOGICAL ::  bio_utci      = .TRUE.   !< Turn index UTCI (instant. input) on or off
+    LOGICAL ::  bio_utci_av   = .TRUE.   !< Turn index UTCI (averaged input) on or off
+
+
+!
+!-- INTERFACES that must be available to other modules (alphabetical order)
+
+    PUBLIC bio_3d_data_averaging, bio_check_data_output,                       &
+    bio_calculate_mrt_grid, bio_calculate_thermal_index_maps, bio_calc_ipt,    &
+    bio_check_parameters, bio_data_output_3d, bio_data_output_2d,              &
+    bio_define_netcdf_grid, bio_get_thermal_index_input_ij, bio_header,        &
+    bio_init, bio_init_arrays, bio_parin, bio_perct, bio_perct_av, bio_pet,    &
+    bio_pet_av, bio_utci, bio_utci_av, time_bio_results
 
 !
@@ -160 +157 @@
 !
 !-- 3D averaging for HTCM _INPUT_ variables
-    INTERFACE biom_3d_data_averaging
-       MODULE PROCEDURE biom_3d_data_averaging
-    END INTERFACE biom_3d_data_averaging
+    INTERFACE bio_3d_data_averaging
+       MODULE PROCEDURE bio_3d_data_averaging
+    END INTERFACE bio_3d_data_averaging
+
+!-- Calculate mtr from rtm fluxes and assign into 2D grid
+    INTERFACE bio_calculate_mrt_grid
+       MODULE PROCEDURE bio_calculate_mrt_grid
+    END INTERFACE bio_calculate_mrt_grid
 
 !-- 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
+    INTERFACE bio_calculate_thermal_index_maps
+       MODULE PROCEDURE bio_calculate_thermal_index_maps
+    END INTERFACE bio_calculate_thermal_index_maps
 
 !-- 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
+    INTERFACE bio_calc_ipt
+       MODULE PROCEDURE bio_calc_ipt
+    END INTERFACE bio_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
+     INTERFACE bio_check_data_output
+        MODULE PROCEDURE bio_check_data_output
+     END INTERFACE bio_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
+    INTERFACE bio_check_parameters
+       MODULE PROCEDURE bio_check_parameters
+    END INTERFACE bio_check_parameters
 
 !-- Data output of 2D quantities
-    INTERFACE biom_data_output_2d
-       MODULE PROCEDURE biom_data_output_2d
-    END INTERFACE biom_data_output_2d
+    INTERFACE bio_data_output_2d
+       MODULE PROCEDURE bio_data_output_2d
+    END INTERFACE bio_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
+    INTERFACE bio_data_output_3d
+       MODULE PROCEDURE bio_data_output_3d
+    END INTERFACE bio_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
+    INTERFACE bio_define_netcdf_grid
+       MODULE PROCEDURE bio_define_netcdf_grid
+    END INTERFACE bio_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
+    INTERFACE bio_get_thermal_index_input_ij
+       MODULE PROCEDURE bio_get_thermal_index_input_ij
+    END INTERFACE bio_get_thermal_index_input_ij
 
 !-- Output of information to the header file
-    INTERFACE biom_header
-       MODULE PROCEDURE biom_header
-    END INTERFACE biom_header
+    INTERFACE bio_header
+       MODULE PROCEDURE bio_header
+    END INTERFACE bio_header
 
 !-- Initialization actions 
-    INTERFACE biom_init
-       MODULE PROCEDURE biom_init
-    END INTERFACE biom_init
+    INTERFACE bio_init
+       MODULE PROCEDURE bio_init
+    END INTERFACE bio_init
 
 !-- Initialization of arrays
-    INTERFACE biom_init_arrays
-       MODULE PROCEDURE biom_init_arrays
-    END INTERFACE biom_init_arrays
+    INTERFACE bio_init_arrays
+       MODULE PROCEDURE bio_init_arrays
+    END INTERFACE bio_init_arrays
 
 !-- Reading of NAMELIST parameters
-    INTERFACE biom_parin
-       MODULE PROCEDURE biom_parin
-    END INTERFACE biom_parin
-
- 
+    INTERFACE bio_parin
+       MODULE PROCEDURE bio_parin
+    END INTERFACE bio_parin
+
+
  CONTAINS
- 
- 
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -234 +236 @@
 !> the array necessary for storing the average.
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_3d_data_averaging( mode, variable )
+ SUBROUTINE bio_3d_data_averaging( mode, variable )
 
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  mode     !< 
-    CHARACTER (LEN=*) ::  variable !< 
-
-    INTEGER(iwp) ::  i  !< 
-    INTEGER(iwp) ::  j  !< 
-    INTEGER(iwp) ::  k  !< 
+    CHARACTER (LEN=*) ::  mode     !< averaging mode: allocate, sum, or average
+    CHARACTER (LEN=*) ::  variable !< The variable in question
+
+    INTEGER(iwp) ::  i        !< Running index, x-dir
+    INTEGER(iwp) ::  j        !< Running index, y-dir
+    INTEGER(iwp) ::  k        !< Running index, z-dir
 
 
@@ -250 +252 @@
        SELECT CASE ( TRIM( variable ) )
 
-          CASE ( 'biom_mrt' )
-                IF ( .NOT. ALLOCATED( biom_mrt_av ) )  THEN
-                   ALLOCATE( biom_mrt_av(nmrtbl) )
+          CASE ( 'bio_mrt' )
+                IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
+                   ALLOCATE( mrt_av_grid(nmrtbl) )
                 ENDIF
-                biom_mrt_av = 0.0_wp
-
-          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
+                mrt_av_grid = 0.0_wp
+
+          CASE ( 'bio_perct*', 'bio_utci*', 'bio_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    &
+             IF ( .NOT. average_trigger_perct .AND. .NOT. average_trigger_utci &
                 .AND. .NOT. average_trigger_pet ) THEN
-                IF ( TRIM( variable ) == 'biom_pt' ) THEN
-                    average_trigger_pt = .TRUE.
+                IF ( TRIM( variable ) == 'bio_perct*' ) THEN
+                    average_trigger_perct = .TRUE.
                 ENDIF
-                IF ( TRIM( variable ) == 'biom_utci' ) THEN
+                IF ( TRIM( variable ) == 'bio_utci*' ) THEN
                     average_trigger_utci = .TRUE.
                 ENDIF
-                IF ( TRIM( variable ) == 'biom_pet' ) THEN
+                IF ( TRIM( variable ) == 'bio_pet*' ) THEN
                     average_trigger_pet = .TRUE.
                 ENDIF
@@ -274 +276 @@
 
 !--          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 ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') RETURN
+             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')   RETURN
+             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_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.
+                aver_perct = .TRUE.
              ENDIF
              pt_av = 0.0_wp
@@ -292 +294 @@
 
              IF ( .NOT. ALLOCATED( u_av ) )  THEN
-                ALLOCATE( u_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 aver_u = .TRUE.
              ENDIF
@@ -298 +300 @@
 
              IF ( .NOT. ALLOCATED( v_av ) )  THEN
-                ALLOCATE( v_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 aver_v = .TRUE.
              ENDIF
@@ -304 +306 @@
 
              IF ( .NOT. ALLOCATED( w_av ) )  THEN
-                ALLOCATE( w_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
+                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 aver_w = .TRUE.
              ENDIF
@@ -318 +320 @@
        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 - degc_to_k
-                    ELSE
-                       biom_mrt_av(:) = biom_mrt_av(:) + &
-                          (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp      &
-                          - degc_to_k
-                    ENDIF
-                 ENDIF
+          CASE ( 'bio_mrt' )
+             IF ( ALLOCATED( mrt_av_grid ) )  THEN
+
+                IF ( mrt_include_sw )  THEN
+                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
+                      (( human_absorb * mrtinsw(:) + human_emiss * mrtinlw(:)) &
+                      / (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
+                ELSE
+                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
+                      (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp          &
+                      - degc_to_k
+                ENDIF
              ENDIF
 
-          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
+          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
 
 !--          Only continue if updateindex
-             IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')    &
+             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
                 RETURN
-             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci')  &
+             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  &
                 RETURN
-             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')   &
+             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')   &
                 RETURN
 
-             IF ( ALLOCATED( pt_av ) .AND. aver_pt ) THEN
+             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -365 +365 @@
 
              IF ( ALLOCATED( u_av )  .AND. aver_u ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
+                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
+                   DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
                          u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
@@ -375 +375 @@
 
              IF ( ALLOCATED( v_av )  .AND. aver_v ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
+                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
+                   DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
                          v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
@@ -385 +385 @@
 
              IF ( ALLOCATED( w_av )  .AND. aver_w ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
+                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
+                   DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
                          w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
@@ -403 +403 @@
        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 )
+          CASE ( 'bio_mrt' )
+             IF ( ALLOCATED( mrt_av_grid ) )  THEN
+                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
              ENDIF
 
-          CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
+          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
 
 !--          Only continue if update index
-             IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')    &
+             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
                 RETURN
-             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci')  &
+             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  &
                 RETURN
-             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')   &
+             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')   &
                 RETURN
 
-             IF ( ALLOCATED( pt_av ) .AND. aver_pt ) THEN
+             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -439 +439 @@
 
              IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
+                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
+                   DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
                          u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND=wp )
@@ -449 +449 @@
 
              IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
                          v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND=wp )
@@ -459 +459 @@
 
              IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
                       DO  k = nzb, nzt+1
                          w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND=wp )
@@ -469 +469 @@
 
 !--          Udate thermal indices with derived averages
-             CALL biom_calculate_static_grid ( .TRUE. )
+             CALL bio_calculate_thermal_index_maps ( .TRUE. )
 
         END SELECT
@@ -476 +476 @@
 
 
- END SUBROUTINE biom_3d_data_averaging
+ END SUBROUTINE bio_3d_data_averaging
 
 
@@ -485 +485 @@
 !> Check data output for biometeorology model
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_check_data_output( var, unit )
+ SUBROUTINE bio_check_data_output( var, unit )
 
     USE control_parameters,                                                    &
@@ -498 +498 @@
     SELECT CASE ( TRIM( var ) )
      
-       CASE( 'biom_tmrt', 'biom_mrt', 'biom_pet', 'biom_pt', 'biom_utci' )
+       CASE( 'bio_mrt', 'bio_pet*', 'bio_perct*', 'bio_utci*' )
           unit = 'degree_C'
 
@@ -509 +509 @@
 !
 !--    Futher checks if output belongs to biometeorology
-       IF ( .NOT. biometeorology .OR. .NOT. radiation ) THEN
-         message_string = 'output of "' // TRIM( var ) // '" req' //  &
-               'uires biometeorology = .TRUE. and radiation = .TRUE. ' &
-               // 'in initialisation_parameters'
-         CALL message( 'biom_check_data_output', 'PA0561', 1, 2, 0, 6, 0 )
+       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'&
+          message_string = 'output of "' // TRIM( var ) // '" require'         &
                            // 's mrt_nlevels > 0'
-          CALL message( 'biom_check_data_output', 'PA0562', 1, 2, 0, 6, 0 )
+          CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
+          unit = 'illegal'
        ENDIF
 
-    ENDIF
-
- END SUBROUTINE biom_check_data_output
+
+    ENDIF
+
+ END SUBROUTINE bio_check_data_output
 
 !------------------------------------------------------------------------------!
@@ -531 +533 @@
 !> check_parameters 1302
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_check_parameters
+ SUBROUTINE bio_check_parameters
 
     USE control_parameters,                                                    &
@@ -541 +543 @@
 
 !--    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 ' // &
+          message_string = 'The estimation of thermal comfort requires '    // &
                            'air humidity information, but humidity module ' // &
                            'is disabled!'
-          CALL message( 'check_parameters', 'PAHU02', 0, 0, 0, 6, 0 )
+          CALL message( 'check_parameters', 'PA0561', 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
-
-    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.
-
-!-- 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.
-    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 )
-            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 = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
-                ENDIF
-                local_pf(i,j,k) = mrt
-            ENDDO
-
-        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  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
-          found = .FALSE.
-
-    END SELECT
-
- END SUBROUTINE biom_data_output_3d
+ END SUBROUTINE bio_check_parameters
+
 
 !------------------------------------------------------------------------------!
@@ -688 +561 @@
 !> data_output_2d 1188ff
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_data_output_2d( av, variable, found, grid, local_pf,          &
+ SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf,           &
                                       two_d, nzb_do, nzt_do, fill_value )
 
@@ -717 +590 @@
     INTEGER(iwp) ::  j        !< Running index, y-dir
     INTEGER(iwp) ::  k        !< Running index, z-dir
-
-
-    found = .TRUE.
+    INTEGER(iwp) ::  l        !< Running index, radiation grid
+
+
     variable_short = TRIM( variable )
-    IF ( variable_short(1:4) == 'biom' ) THEN
-       two_d = .TRUE.
-       grid = 'zu1'
-    ELSE
+    IF ( variable_short(1:4) /= 'bio_' ) THEN
        found = .FALSE.
        grid  = 'none'
-       RETURN
-    ENDIF
-
+    ENDIF
+
+    found = .TRUE.
     local_pf = fill_value
 
@@ -735 +605 @@
 
 
-        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
+        CASE ( 'bio_mrt_xy' )
+            grid = 'zu1'
+            two_d = .FALSE.  !< can be calculated for several levels
+            local_pf = REAL( fill_value, KIND = wp )
+            DO  l = 1, nmrtbl
+               i = mrtbl(ix,l)
+               j = mrtbl(iy,l)
+               k = mrtbl(iz,l)
+               IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.  &
+                  i < nxl .OR. i > nxr ) CYCLE
+               IF ( av == 0 )  THEN
+                  IF ( mrt_include_sw )  THEN
+                     local_pf(i,j,k) = ((human_absorb * mrtinsw(l) +           &
+                     human_emiss * mrtinlw(l)) /                               &
+                     (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
+                  ELSE
+                     local_pf(i,j,k) = (human_emiss * mrtinlw(l) /             &
+                                        sigma_sb) ** .25_wp - degc_to_k
+                  ENDIF
+               ELSE
+                  local_pf(i,j,k) = mrt_av_grid(l)
+               ENDIF
+            ENDDO
+
+
+        CASE ( 'bio_perct*_xy' )        ! 2d-array
+            grid = 'zu1'
+            two_d = .TRUE.
             IF ( av == 0 )  THEN
               DO  i = nxl, nxr
                  DO  j = nys, nyn
-                    local_pf(i,j,nzb+1) = pt_grid(j,i)
+                    local_pf(i,j,nzb+1) = perct(j,i)
                  ENDDO
               ENDDO
@@ -753 +642 @@
               DO  i = nxl, nxr
                  DO  j = nys, nyn
-                    local_pf(i,j,nzb+1) = pt_av_grid(j,i)
+                    local_pf(i,j,nzb+1) = perct_av(j,i)
                  ENDDO
               ENDDO
@@ -759 +648 @@
 
 
-        CASE ( 'biom_utci_xy' )        ! 2d-array
+        CASE ( 'bio_utci*_xy' )        ! 2d-array
+            grid = 'zu1'
+            two_d = .TRUE.
             IF ( av == 0 )  THEN
               DO  i = nxl, nxr
@@ -775 +666 @@
 
 
-        CASE ( 'biom_pet_xy' )        ! 2d-array
+        CASE ( 'bio_pet*_xy' )        ! 2d-array
+            grid = 'zu1'
+            two_d = .TRUE.
             IF ( av == 0 )  THEN
               DO  i = nxl, nxr
@@ -798 +691 @@
 
 
- END SUBROUTINE biom_data_output_2d
-
+ END SUBROUTINE bio_data_output_2d
+
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables (dummy, always 2d!)
+!> data_output_3d 709ff
+!------------------------------------------------------------------------------!
+ SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
+
+    USE indices
+
+    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.
+
+!-- 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.
+    variable_short = TRIM( variable )
+
+    IF ( variable_short(1:4) /= 'bio_' ) THEN
+!--   Nothing to do, set found to FALSE and return immediatelly
+      found = .FALSE.
+      RETURN
+    ENDIF
+
+    SELECT CASE ( variable_short )
+
+        CASE ( 'bio_mrt' )
+            local_pf = REAL( fill_value, KIND = wp )
+            DO  l = 1, nmrtbl
+               i = mrtbl(ix,l)
+               j = mrtbl(iy,l)
+               k = mrtbl(iz,l)
+               IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.  &
+                  i < nxl .OR. i > nxr ) CYCLE
+               IF ( av == 0 )  THEN
+                  IF ( mrt_include_sw )  THEN
+                     local_pf(i,j,k) = ((human_absorb * mrtinsw(l) + human_emiss * mrtinlw(l)) /   &
+                                        (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
+                  ELSE
+                     local_pf(i,j,k) = (human_emiss * mrtinlw(l) /             &
+                                        sigma_sb) ** .25_wp - degc_to_k
+                  ENDIF
+               ELSE
+                  local_pf(i,j,k) = mrt_av_grid(l)
+               ENDIF
+            ENDDO
+
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+ END SUBROUTINE bio_data_output_3d
 
 !------------------------------------------------------------------------------!
@@ -808 +777 @@
 !> netcdf_interface_mod 918ff
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
+ SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
 
     IMPLICIT NONE
@@ -836 +805 @@
     is2d = ( var(l-1:l) == 'xy' )
 
-
-    IF ( var(1:4) == 'biom' ) THEN
+    IF ( var(1:4) == 'bio_' ) THEN
        found  = .TRUE.
        grid_x = 'x'
@@ -845 +813 @@
     ENDIF
 
- END SUBROUTINE biom_define_netcdf_grid
+ END SUBROUTINE bio_define_netcdf_grid
 
 !------------------------------------------------------------------------------!
@@ -853 +821 @@
 !> header 982
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_header( io )
+ SUBROUTINE bio_header( io )
 
     IMPLICIT NONE
@@ -863 +831 @@
     CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
 
-    WRITE( output_height_chr, '(F8.1,7X)' )  biom_output_height
+    WRITE( output_height_chr, '(F8.1,7X)' )  bio_output_height
 !
 !-- Write biom header
     WRITE( io, 1 )
     WRITE( io, 2 )  TRIM( output_height_chr )
-    WRITE( io, 3 )  TRIM( ACHAR( biom_cell_level ) )
+    WRITE( io, 3 )  TRIM( ACHAR( bio_cell_level ) )
 
 1   FORMAT (//' Human thermal comfort module information:'/                    &
@@ -875 +843 @@
 3   FORMAT ('    --> This corresponds to cell level : ', A )
 
- END SUBROUTINE biom_header
+ END SUBROUTINE bio_header
 
 
@@ -884 +852 @@
 !> init_3d_model 1987ff
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_init
-
-    USE control_parameters,                                                 &
+ SUBROUTINE bio_init
+
+    USE control_parameters,                                                    &
         ONLY: message_string
 
@@ -899 +867 @@
 !   (gravimetric center of sample human)
 
-    time_biom_results = 0.0_wp
-    biom_cell_level = 0_iwp
-    biom_output_height = 0.5_wp * dz(1)
+    time_bio_results = 0.0_wp
+    bio_cell_level = 0_iwp
+    bio_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)
-
-    IF ( .NOT. radiation_interactions .AND. mrt_from_rtm ) THEN
-       message_string = 'The mrt_from_rtm switch require ' // &
+    bio_cell_level = INT ( 1.099_wp / dz(1) )
+    bio_output_height = bio_output_height + bio_cell_level * dz(1)
+
+    IF ( .NOT. radiation_interactions ) THEN
+       message_string = 'The mrt calculation requires ' // &
                         'enabled radiation_interactions but it ' // &
-                        'is disabled! Set mrt_from_rtm to .F.'
+                        'is disabled!'
        CALL message( 'check_parameters', 'PAHU03', 0, 0, -1, 6, 0 )
-       mrt_from_rtm = .FALSE.
-    ENDIF
-
- END SUBROUTINE biom_init
+    ENDIF
+
+ END SUBROUTINE bio_init
 
 
@@ -924 +891 @@
 !> init_3d_model 1047ff
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_init_arrays
+ SUBROUTINE bio_init_arrays
 
     IMPLICIT NONE
@@ -934 +901 @@
     ENDIF
 
-    IF ( biom_pt ) THEN
-      IF ( .NOT. ALLOCATED( pt_grid ) ) THEN
-        ALLOCATE( pt_grid (nys:nyn,nxl:nxr) )
+    IF ( bio_perct ) THEN
+      IF ( .NOT. ALLOCATED( perct ) ) THEN
+        ALLOCATE( perct (nys:nyn,nxl:nxr) )
       ENDIF
     ENDIF
 
-    IF ( biom_utci ) THEN
+    IF ( bio_utci ) THEN
       IF ( .NOT. ALLOCATED( utci_grid ) ) THEN
         ALLOCATE( utci_grid (nys:nyn,nxl:nxr) )
@@ -946 +913 @@
     ENDIF
 
-    IF ( biom_pet ) THEN
+    IF ( bio_pet ) THEN
       IF ( .NOT. ALLOCATED( pet_grid ) ) THEN
         ALLOCATE( pet_grid (nys:nyn,nxl:nxr) )
@@ -952 +919 @@
     END IF
 
-    IF ( biom_pt_av ) THEN
-      IF ( .NOT. ALLOCATED( pt_av_grid ) ) THEN
-        ALLOCATE( pt_av_grid (nys:nyn,nxl:nxr) )
+    IF ( bio_perct_av ) THEN
+      IF ( .NOT. ALLOCATED( perct_av ) ) THEN
+        ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
       ENDIF
     ENDIF
 
-    IF ( biom_utci_av ) THEN
+    IF ( bio_utci_av ) THEN
       IF ( .NOT. ALLOCATED( utci_av_grid ) ) THEN
         ALLOCATE( utci_av_grid (nys:nyn,nxl:nxr) )
@@ -964 +931 @@
     ENDIF
 
-    IF ( biom_pet_av ) THEN
+    IF ( bio_pet_av ) THEN
       IF ( .NOT. ALLOCATED( pet_av_grid ) ) THEN
         ALLOCATE( pet_av_grid (nys:nyn,nxl:nxr) )
@@ -971 +938 @@
     END IF
 
- END SUBROUTINE biom_init_arrays
+ END SUBROUTINE bio_init_arrays
 
 
@@ -979 +946 @@
 !> Parin for &biometeorology_parameters for reading biomet parameters
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_parin
+ SUBROUTINE bio_parin
 
     IMPLICIT NONE
@@ -987 +954 @@
     CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file 
 
-    NAMELIST /biometeorology_parameters/  mrt_from_rtm,                        &
-                                          biom_pet,                            &
-                                          biom_pet_av,                         &
-                                          biom_pt,                             &
-                                          biom_pt_av,                          &
-                                          biom_utci,                           &
-                                          biom_utci_av
+    NAMELIST /biometeorology_parameters/  bio_pet,                             &
+                                          bio_pet_av,                          &
+                                          bio_perct,                           &
+                                          bio_perct_av,                        &
+                                          bio_utci,                            &
+                                          bio_utci_av
 
 
@@ -1025 +991 @@
 
 
- END SUBROUTINE biom_parin
+ END SUBROUTINE bio_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 )
+!> Calculate biometeorology MRT for all 2D grid
+!------------------------------------------------------------------------------!
+ SUBROUTINE bio_calculate_mrt_grid ( av )
 
     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 - degc_to_k
-
- 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
+    LOGICAL, INTENT(IN)         ::  av    !< use averaged input?
+!-- Internal variables
+    INTEGER(iwp)                ::  i     !< Running index, x-dir, radiation coordinates
+    INTEGER(iwp)                ::  j     !< Running index, y-dir, radiation coordinates
+    INTEGER(iwp)                ::  k     !< Running index, y-dir, radiation coordinates
+    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
+
+
+!-- Calculate biometeorology MRT from local radiation
+!-- fluxes calculated by RTM and assign into 2D grid
+    tmrt_grid = -999.
+    DO  l = 1, nmrtbl
+       i = mrtbl(ix,l)
+       j = mrtbl(iy,l)
+       k = mrtbl(iz,l)
+       IF ( k - get_topography_top_index_ji( j, i, 's' ) == bio_cell_level +   &
+             1_iwp) THEN
+          IF ( mrt_include_sw )  THEN
+              tmrt_grid(j,i) = ((human_absorb*mrtinsw(l) +                     &
+                                human_emiss*mrtinlw(l))  /                     &
+                               (human_emiss*sigma_sb)) ** .25_wp - degc_to_k
+          ELSE
+              tmrt_grid(j,i) = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp   &
+                                 - degc_to_k
+          ENDIF
        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 + degc_to_k )**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 - degc_to_k
-
- 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
+    ENDDO
+
+END SUBROUTINE bio_calculate_mrt_grid
 
 
@@ -1215 +1038 @@
 !> Calculate static thermal indices for 2D grid point i, j
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_determine_input_at( average_input, i, j, ta, vp, ws, pair,    &
-    tmrt )
+ SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws,   &
+    pair, tmrt )
 
     IMPLICIT NONE
@@ -1234 +1057 @@
 !-- Internal variables
     INTEGER(iwp)                ::  k     !< Running index, z-dir
+    INTEGER(iwp)                ::  ir    !< Running index, x-dir, radiation coordinates
+    INTEGER(iwp)                ::  jr    !< Running index, y-dir, radiation coordinates
+    INTEGER(iwp)                ::  kr    !< Running index, y-dir, radiation coordinates
     INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
+    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
 
     REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
@@ -1241 +1068 @@
 !-- 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 = get_topography_top_index_ji(j, i, 's') + bio_cell_level  !< Vertical cell center closest to 1.1m
+
+!
+!-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground
     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
+    IF ( bio_cell_level < 1_iwp ) THEN
+       k_wind = k + 1_iwp
     ENDIF
 
@@ -1252 +1082 @@
        ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
 
-       vp = 0.034_wp
+       vp = -999._wp
        IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
           vp = q_av(k, j, i)
@@ -1264 +1094 @@
        ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
 
-       vp = q(k, j, i)
+       vp = -999._wp
+       IF ( humidity ) THEN
+          vp = q(k, j, i)
+       ENDIF
 
        ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) )  +           &
@@ -1278 +1111 @@
 !-- 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 + degc_to_k, 333.15_wp ) )
-    vp  = vp * pair / ( vp + 0.622_wp )
-    IF ( vp > vp_sat ) vp = vp_sat
-
-    tmrt = ta
+    IF ( vp > -998._wp ) THEN
+       vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
+       vp  = vp * pair / ( vp + 0.622_wp )
+       IF ( vp > vp_sat ) vp = vp_sat
+    ENDIF
+
+!-- local mtr value at [i,j]
+    tmrt = -999.  !< this can be a valid result (e.g. for inside some ostacle)
     IF ( radiation ) THEN
-       IF ( mrt_from_rtm ) THEN
-          tmrt = tmrt_grid(j, i)
-       ELSE
-          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
-    ENDIF
-
- END SUBROUTINE biom_determine_input_at
+!--    Use MRT from RTM precalculated in tmrt_grid
+       tmrt = tmrt_grid(j,i)
+    ENDIF
+
+ END SUBROUTINE bio_get_thermal_index_input_ij
 
 
@@ -1302 +1133 @@
 !> time_integration.f90: 1065ff
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_calculate_static_grid ( average_input )
+ SUBROUTINE bio_calculate_thermal_index_maps ( av )
 
     IMPLICIT NONE
 
 !-- Input attributes
-    LOGICAL, INTENT ( IN ) ::  average_input  !< Calculate based on averaged input? conditions?
+    LOGICAL, INTENT ( IN ) ::  av  !< Calculate based on averaged input conditions?
 
 !-- Internal variables
-    INTEGER(iwp) ::  i, j, k, l   !< Running index
-
+    INTEGER(iwp) ::  i, j      !< Running index
+
+    REAL(wp) ::  clo           !< Clothing index                (no dimension)
     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
-
-    IF ( mrt_from_rtm ) THEN
-       tmrt_grid = -999._wp
-       DO  l = 1, nmrtbl
-           i = mrtbl(ix,l)
-           j = mrtbl(iy,l)
-           k = mrtbl(iz,l)
-           IF ( k - get_topography_top_index_ji( j, i, 's' ) == 1  ) THEN
-              IF ( mrt_include_sw )  THEN
-                  tmrt_tmp = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
-                             / (human_emiss*sigma_sb)) ** .25_wp
-              ELSE
-                  tmrt_tmp = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
-              ENDIF
-              tmrt_grid(j, i) = tmrt_tmp - degc_to_k
-           ENDIF
-       ENDDO
-    ENDIF
+    REAL(wp) ::  perct_tmp     !< Perceived temperature            (°C)
+    REAL(wp) ::  utci_tmp      !< Universal thermal climate index  (°C)
+    REAL(wp) ::  pet_tmp       !< Physiologically equivalent temperature  (°C)
+    REAL(wp) ::  tmrt_tmp      !< Mean radiant temperature         (°C)
+
+    CALL bio_init_arrays
+
+!-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
+    CALL bio_calculate_mrt_grid ( av )
+
 
     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
+          CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,              &
+               ws, pair, tmrt_tmp )
+!           tmrt_grid(j, i) = tmrt_tmp  !< already set by bio_calculate_mrt_grid!
+
+!--       Only proceed if input is available
+          IF ( tmrt_tmp <= -998._wp .OR. vp <= -998._wp ) THEN
+             pet_tmp = -999._wp         !< set fail value, e.g. valid for within
+             perct_tmp  = -999._wp      !< some obstacle
              utci_tmp = -999._wp
-             pet_tmp  = -999._wp
-             CYCLE
+          ELSE
+!--          Calculate static thermal indices based on local tmrt
+             clo = -999._wp
+             
+             IF ( bio_perct ) THEN
+!--          Estimate local perceived temperature
+                CALL calculate_perct_static( ta, vp, ws, tmrt_tmp, pair, clo,  &
+                   perct_tmp )
+             ENDIF
+             
+             IF ( bio_utci ) THEN
+!--          Estimate local universal thermal climate index
+                CALL calculate_utci_static( ta, vp, ws, tmrt_tmp,              &
+                   bio_output_height, utci_tmp )
+             ENDIF
+             
+             IF ( bio_pet ) THEN
+!--          Estimate local physiologically equivalent temperature
+                CALL calculate_pet_static( ta, vp, ws, tmrt_tmp, pair, pet_tmp )
+             ENDIF
           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
+
+          IF ( av ) THEN
 !--          Write results for selected averaged indices
-             IF ( biom_pt_av )  THEN
-                pt_av_grid(j, i)   = pt_tmp
+             IF ( bio_perct_av )  THEN
+                perct_av(j, i) = perct_tmp
              END IF
-             IF ( biom_utci_av ) THEN
+             IF ( bio_utci_av ) THEN
                 utci_av_grid(j, i) = utci_tmp
              END IF
-             IF ( biom_pet_av ) THEN
+             IF ( bio_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
+             IF ( bio_perct )  THEN
+                perct(j, i) = perct_tmp
              END IF
-             IF ( biom_utci ) THEN
+             IF ( bio_utci ) THEN
                 utci_grid(j, i) = utci_tmp
              END IF
-             IF ( biom_pet ) THEN
+             IF ( bio_pet ) THEN
                 pet_grid(j, i)  = pet_tmp
              END IF
@@ -1388 +1221 @@
     END DO
 
- END SUBROUTINE biom_calculate_static_grid
+ END SUBROUTINE bio_calculate_thermal_index_maps
 
 !------------------------------------------------------------------------------!
@@ -1395 +1228 @@
 !> Calculate dynamic thermal indices (currently only iPT, but expandable)
 !------------------------------------------------------------------------------!
- SUBROUTINE biom_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,         &
+ SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,          &
     t_clo, clo, actlev, age, weight, height, work, sex, ipt )
 
@@ -1436 +1269 @@
     ENDIF
 
- END SUBROUTINE biom_calc_ipt
+ END SUBROUTINE bio_calc_ipt
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> SUBROUTINE for calculating UTCI Temperature (UTCI)
+!> computed by a 6th order approximation
+!
+!> 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 source available at:
+!> www.utci.org
+!------------------------------------------------------------------------------!
+ SUBROUTINE calculate_utci_static( ta_in, vp, ws_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 )  ::  vp     !< Loacl vapour pressure (hPa)
+    REAL(WP), INTENT ( IN )  ::  ws_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
+    REAL(WP) ::  ta         !< internal air temperature (°C)
+    REAL(WP) ::  pa         !< air pressure in kPa      (kPa)
+    REAL(WP) ::  d_tmrt     !< delta-tmrt               (°C)
+    REAL(WP) ::  va         !< wind speed at 10 m above ground level    (m/s)
+    REAL(WP) ::  offset     !< utci deviation by ta cond. exceeded      (°C)
+
+!-- Initialize
+    offset = 0._wp
+    ta = ta_in
+    d_tmrt = tmrt - ta_in
+
+!-- Use vapour pressure in kpa
+    pa = vp / 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 profile definition
+    va = ws_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._wp ) .OR. ( d_tmrt < -30._wp ) .OR.                    &
+       ( vp >= 50._wp ) ) THEN
+       utci = -999._wp
+       RETURN
+    ENDIF
+
+!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
+    IF ( ta > 50._wp ) THEN
+       offset = ta - 50._wp
+       ta = 50._wp
+    ENDIF
+    IF ( ta < -50._wp ) THEN
+       offset = ta + 50._wp
+       ta = -50._wp
+    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.5_wp ) va = 0.5_wp
+    IF ( va > 17._wp ) va = 17._wp
+
+!-- Calculate 6th order polynomial as approximation 
+    utci = ta +                                                                &
+       ( 6.07562052e-01 )   +                                                  &
+       ( -2.27712343e-02 ) * ta +                                              &
+       ( 8.06470249e-04 )  * ta * ta +                                         &
+       ( -1.54271372e-04 ) * ta * ta * ta +                                    &
+       ( -3.24651735e-06 ) * ta * ta * ta * ta +                               &
+       ( 7.32602852e-08 )  * ta * ta * ta * ta * ta +                          &
+       ( 1.35959073e-09 )  * ta * ta * ta * ta * ta * ta +                     &
+       ( -2.25836520e+00 ) * va +                                              &
+       ( 8.80326035e-02 )  * ta * va +                                         &
+       ( 2.16844454e-03 )  * ta * ta * va +                                    &
+       ( -1.53347087e-05 ) * ta * ta * ta * va +                               &
+       ( -5.72983704e-07 ) * ta * ta * ta * ta * va +                          &
+       ( -2.55090145e-09 ) * ta * ta * ta * ta * ta * va +                     &
+       ( -7.51269505e-01 ) * va * va +                                         &
+       ( -4.08350271e-03 ) * ta * va * va +                                    &
+       ( -5.21670675e-05 ) * ta * ta * va * va +                               &
+       ( 1.94544667e-06 )  * ta * ta * ta * va * va +                          &
+       ( 1.14099531e-08 )  * ta * ta * ta * ta * va * va +                     &
+       ( 1.58137256e-01 )  * va * va * va +                                    &
+       ( -6.57263143e-05 ) * ta * va * va * va +                               &
+       ( 2.22697524e-07 )  * ta * ta * va * va * va +                          &
+       ( -4.16117031e-08 ) * ta * ta * ta * va * va * va +                     &
+       ( -1.27762753e-02 ) * va * va * va * va +                               &
+       ( 9.66891875e-06 )  * ta * va * va * va * va +                          &
+       ( 2.52785852e-09 )  * ta * ta * va * va * va * va +                     &
+       ( 4.56306672e-04 )  * va * va * va * va * va +                          &
+       ( -1.74202546e-07 ) * ta * va * va * va * va * va +                     &
+       ( -5.91491269e-06 ) * va * va * va * va * va * va +                     &
+       ( 3.98374029e-01 )  * d_tmrt +                                          &
+       ( 1.83945314e-04 )  * ta * d_tmrt +                                     &
+       ( -1.73754510e-04 ) * ta * ta * d_tmrt +                                &
+       ( -7.60781159e-07 ) * ta * ta * ta * d_tmrt +                           &
+       ( 3.77830287e-08 )  * ta * ta * ta * ta * d_tmrt +                      &
+       ( 5.43079673e-10 )  * ta * ta * ta * ta * ta * d_tmrt +                 &
+       ( -2.00518269e-02 ) * va * d_tmrt +                                     &
+       ( 8.92859837e-04 )  * ta * va * d_tmrt +                                &
+       ( 3.45433048e-06 )  * ta * ta * va * d_tmrt +                           &
+       ( -3.77925774e-07 ) * ta * ta * ta * va * d_tmrt +                      &
+       ( -1.69699377e-09 ) * ta * ta * ta * ta * va * d_tmrt +                 &
+       ( 1.69992415e-04 )  * va * va * d_tmrt +                                &
+       ( -4.99204314e-05 ) * ta * va * va * d_tmrt +                           &
+       ( 2.47417178e-07 )  * ta * ta * va * va * d_tmrt +                      &
+       ( 1.07596466e-08 )  * ta * ta * ta * va * va * d_tmrt +                 &
+       ( 8.49242932e-05 )  * va * va * va * d_tmrt +                           &
+       ( 1.35191328e-06 )  * ta * va * va * va * d_tmrt +                      &
+       ( -6.21531254e-09 ) * ta * ta * va * va * va * d_tmrt +                 &
+       ( -4.99410301e-06 ) * va * va * va * va * d_tmrt +                      &
+       ( -1.89489258e-08 ) * ta * va * va * va * va * d_tmrt +                 &
+       ( 8.15300114e-08 )  * va * va * va * va * va * d_tmrt +                 &
+       ( 7.55043090e-04 )  * d_tmrt * d_tmrt +                                 &
+       ( -5.65095215e-05 ) * ta * d_tmrt * d_tmrt +                            &
+       ( -4.52166564e-07 ) * ta * ta * d_tmrt * d_tmrt +                       &
+       ( 2.46688878e-08 )  * ta * ta * ta * d_tmrt * d_tmrt +                  &
+       ( 2.42674348e-10 )  * ta * ta * ta * ta * d_tmrt * d_tmrt +             &
+       ( 1.54547250e-04 )  * va * d_tmrt * d_tmrt +                            &
+       ( 5.24110970e-06 )  * ta * va * d_tmrt * d_tmrt +                       &
+       ( -8.75874982e-08 ) * ta * ta * va * d_tmrt * d_tmrt +                  &
+       ( -1.50743064e-09 ) * ta * ta * ta * va * d_tmrt * d_tmrt +             &
+       ( -1.56236307e-05 ) * va * va * d_tmrt * d_tmrt +                       &
+       ( -1.33895614e-07 ) * ta * va * va * d_tmrt * d_tmrt +                  &
+       ( 2.49709824e-09 )  * ta * ta * va * va * d_tmrt * d_tmrt +             &
+       ( 6.51711721e-07 )  * va * va * va * d_tmrt * d_tmrt +                  &
+       ( 1.94960053e-09 )  * ta * va * va * va * d_tmrt * d_tmrt +             &
+       ( -1.00361113e-08 ) * va * va * va * va * d_tmrt * d_tmrt +             &
+       ( -1.21206673e-05 ) * d_tmrt * d_tmrt * d_tmrt +                        &
+       ( -2.18203660e-07 ) * ta * d_tmrt * d_tmrt * d_tmrt +                   &
+       ( 7.51269482e-09 )  * ta * ta * d_tmrt * d_tmrt * d_tmrt +              &
+       ( 9.79063848e-11 )  * ta * ta * ta * d_tmrt * d_tmrt * d_tmrt +         &
+       ( 1.25006734e-06 )  * va * d_tmrt * d_tmrt * d_tmrt +                   &
+       ( -1.81584736e-09 ) * ta * va * d_tmrt * d_tmrt * d_tmrt +              &
+       ( -3.52197671e-10 ) * ta * ta * va * d_tmrt * d_tmrt * d_tmrt +         &
+       ( -3.36514630e-08 ) * va * va * d_tmrt * d_tmrt * d_tmrt +              &
+       ( 1.35908359e-10 )  * ta * va * va * d_tmrt * d_tmrt * d_tmrt +         &
+       ( 4.17032620e-10 )  * va * va * va * d_tmrt * d_tmrt * d_tmrt +         &
+       ( -1.30369025e-09 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt +               &
+       ( 4.13908461e-10 )  * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt +          &
+       ( 9.22652254e-12 )  * ta * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
+       ( -5.08220384e-09 ) * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +          &
+       ( -2.24730961e-11 ) * ta * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
+       ( 1.17139133e-10 )  * va * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
+       ( 6.62154879e-10 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt +      &
+       ( 4.03863260e-13 )  * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + &
+       ( 1.95087203e-12 )  * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + &
+       ( -4.73602469e-12 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt *      &
+       d_tmrt +                                                                &
+       ( 5.12733497e+00 )  * pa +                                              &
+       ( -3.12788561e-01 ) * ta * pa +                                         &
+       ( -1.96701861e-02 ) * ta * ta * pa +                                    &
+       ( 9.99690870e-04 )  * ta * ta * ta * pa +                               &
+       ( 9.51738512e-06 )  * ta * ta * ta * ta * pa +                          &
+       ( -4.66426341e-07 ) * ta * ta * ta * ta * ta * pa +                     &
+       ( 5.48050612e-01 )  * va * pa +                                         &
+       ( -3.30552823e-03 ) * ta * va * pa +                                    &
+       ( -1.64119440e-03 ) * ta * ta * va * pa +                               &
+       ( -5.16670694e-06 ) * ta * ta * ta * va * pa +                          &
+       ( 9.52692432e-07 )  * ta * ta * ta * ta * va * pa +                     &
+       ( -4.29223622e-02 ) * va * va * pa +                                    &
+       ( 5.00845667e-03 )  * ta * va * va * pa +                               &
+       ( 1.00601257e-06 )  * ta * ta * va * va * pa +                          &
+       ( -1.81748644e-06 ) * ta * ta * ta * va * va * pa +                     &
+       ( -1.25813502e-03 ) * va * va * va * pa +                               &
+       ( -1.79330391e-04 ) * ta * va * va * va * pa +                          &
+       ( 2.34994441e-06 )  * ta * ta * va * va * va * pa +                     &
+       ( 1.29735808e-04 )  * va * va * va * va * pa +                          &
+       ( 1.29064870e-06 )  * ta * va * va * va * va * pa +                     &
+       ( -2.28558686e-06 ) * va * va * va * va * va * pa +                     &
+       ( -3.69476348e-02 ) * d_tmrt * pa +                                     &
+       ( 1.62325322e-03 )  * ta * d_tmrt * pa +                                &
+       ( -3.14279680e-05 ) * ta * ta * d_tmrt * pa +                           &
+       ( 2.59835559e-06 )  * ta * ta * ta * d_tmrt * pa +                      &
+       ( -4.77136523e-08 ) * ta * ta * ta * ta * d_tmrt * pa +                 &
+       ( 8.64203390e-03 )  * va * d_tmrt * pa +                                &
+       ( -6.87405181e-04 ) * ta * va * d_tmrt * pa +                           &
+       ( -9.13863872e-06 ) * ta * ta * va * d_tmrt * pa +                      &
+       ( 5.15916806e-07 )  * ta * ta * ta * va * d_tmrt * pa +                 &
+       ( -3.59217476e-05 ) * va * va * d_tmrt * pa +                           &
+       ( 3.28696511e-05 )  * ta * va * va * d_tmrt * pa +                      &
+       ( -7.10542454e-07 ) * ta * ta * va * va * d_tmrt * pa +                 &
+       ( -1.24382300e-05 ) * va * va * va * d_tmrt * pa +                      &
+       ( -7.38584400e-09 ) * ta * va * va * va * d_tmrt * pa +                 &
+       ( 2.20609296e-07 )  * va * va * va * va * d_tmrt * pa +                 &
+       ( -7.32469180e-04 ) * d_tmrt * d_tmrt * pa +                            &
+       ( -1.87381964e-05 ) * ta * d_tmrt * d_tmrt * pa +                       &
+       ( 4.80925239e-06 )  * ta * ta * d_tmrt * d_tmrt * pa +                  &
+       ( -8.75492040e-08 ) * ta * ta * ta * d_tmrt * d_tmrt * pa +             &
+       ( 2.77862930e-05 )  * va * d_tmrt * d_tmrt * pa +                       &
+       ( -5.06004592e-06 ) * ta * va * d_tmrt * d_tmrt * pa +                  &
+       ( 1.14325367e-07 )  * ta * ta * va * d_tmrt * d_tmrt * pa +             &
+       ( 2.53016723e-06 )  * va * va * d_tmrt * d_tmrt * pa +                  &
+       ( -1.72857035e-08 ) * ta * va * va * d_tmrt * d_tmrt * pa +             &
+       ( -3.95079398e-08 ) * va * va * va * d_tmrt * d_tmrt * pa +             &
+       ( -3.59413173e-07 ) * d_tmrt * d_tmrt * d_tmrt * pa +                   &
+       ( 7.04388046e-07 )  * ta * d_tmrt * d_tmrt * d_tmrt * pa +              &
+       ( -1.89309167e-08 ) * ta * ta * d_tmrt * d_tmrt * d_tmrt * pa +         &
+       ( -4.79768731e-07 ) * va * d_tmrt * d_tmrt * d_tmrt * pa +              &
+       ( 7.96079978e-09 )  * ta * va * d_tmrt * d_tmrt * d_tmrt * pa +         &
+       ( 1.62897058e-09 )  * va * va * d_tmrt * d_tmrt * d_tmrt * pa +         &
+       ( 3.94367674e-08 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +          &
+       ( -1.18566247e-09 ) * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +     &
+       ( 3.34678041e-10 )  * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +     &
+       ( -1.15606447e-10 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa + &
+       ( -2.80626406e+00 ) * pa * pa +                                         &
+       ( 5.48712484e-01 )  * ta * pa * pa +                                    &
+       ( -3.99428410e-03 ) * ta * ta * pa * pa +                               &
+       ( -9.54009191e-04 ) * ta * ta * ta * pa * pa +                          &
+       ( 1.93090978e-05 )  * ta * ta * ta * ta * pa * pa +                     &
+       ( -3.08806365e-01 ) * va * pa * pa +                                    &
+       ( 1.16952364e-02 )  * ta * va * pa * pa +                               &
+       ( 4.95271903e-04 )  * ta * ta * va * pa * pa +                          &
+       ( -1.90710882e-05 ) * ta * ta * ta * va * pa * pa +                     &
+       ( 2.10787756e-03 )  * va * va * pa * pa +                               &
+       ( -6.98445738e-04 ) * ta * va * va * pa * pa +                          &
+       ( 2.30109073e-05 )  * ta * ta * va * va * pa * pa +                     &
+       ( 4.17856590e-04 )  * va * va * va * pa * pa +                          &
+       ( -1.27043871e-05 ) * ta * va * va * va * pa * pa +                     &
+       ( -3.04620472e-06 ) * va * va * va * va * pa * pa +                     &
+       ( 5.14507424e-02 )  * d_tmrt * pa * pa +                                &
+       ( -4.32510997e-03 ) * ta * d_tmrt * pa * pa +                           &
+       ( 8.99281156e-05 )  * ta * ta * d_tmrt * pa * pa +                      &
+       ( -7.14663943e-07 ) * ta * ta * ta * d_tmrt * pa * pa +                 &
+       ( -2.66016305e-04 ) * va * d_tmrt * pa * pa +                           &
+       ( 2.63789586e-04 )  * ta * va * d_tmrt * pa * pa +                      &
+       ( -7.01199003e-06 ) * ta * ta * va * d_tmrt * pa * pa +                 &
+       ( -1.06823306e-04 ) * va * va * d_tmrt * pa * pa +                      &
+       ( 3.61341136e-06 )  * ta * va * va * d_tmrt * pa * pa +                 &
+       ( 2.29748967e-07 )  * va * va * va * d_tmrt * pa * pa +                 &
+       ( 3.04788893e-04 )  * d_tmrt * d_tmrt * pa * pa +                       &
+       ( -6.42070836e-05 ) * ta * d_tmrt * d_tmrt * pa * pa +                  &
+       ( 1.16257971e-06 )  * ta * ta * d_tmrt * d_tmrt * pa * pa +             &
+       ( 7.68023384e-06 )  * va * d_tmrt * d_tmrt * pa * pa +                  &
+       ( -5.47446896e-07 ) * ta * va * d_tmrt * d_tmrt * pa * pa +             &
+       ( -3.59937910e-08 ) * va * va * d_tmrt * d_tmrt * pa * pa +             &
+       ( -4.36497725e-06 ) * d_tmrt * d_tmrt * d_tmrt * pa * pa +              &
+       ( 1.68737969e-07 )  * ta * d_tmrt * d_tmrt * d_tmrt * pa * pa +         &
+       ( 2.67489271e-08 )  * va * d_tmrt * d_tmrt * d_tmrt * pa * pa +         &
+       ( 3.23926897e-09 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa * pa +     &
+       ( -3.53874123e-02 ) * pa * pa * pa +                                    &
+       ( -2.21201190e-01 ) * ta * pa * pa * pa +                               &
+       ( 1.55126038e-02 )  * ta * ta * pa * pa * pa +                          &
+       ( -2.63917279e-04 ) * ta * ta * ta * pa * pa * pa +                     &
+       ( 4.53433455e-02 )  * va * pa * pa * pa +                               &
+       ( -4.32943862e-03 ) * ta * va * pa * pa * pa +                          &
+       ( 1.45389826e-04 )  * ta * ta * va * pa * pa * pa +                     &
+       ( 2.17508610e-04 )  * va * va * pa * pa * pa +                          &
+       ( -6.66724702e-05 ) * ta * va * va * pa * pa * pa +                     &
+       ( 3.33217140e-05 )  * va * va * va * pa * pa * pa +                     &
+       ( -2.26921615e-03 ) * d_tmrt * pa * pa * pa +                           &
+       ( 3.80261982e-04 )  * ta * d_tmrt * pa * pa * pa +                      &
+       ( -5.45314314e-09 ) * ta * ta * d_tmrt * pa * pa * pa +                 &
+       ( -7.96355448e-04 ) * va * d_tmrt * pa * pa * pa +                      &
+       ( 2.53458034e-05 )  * ta * va * d_tmrt * pa * pa * pa +                 &
+       ( -6.31223658e-06 ) * va * va * d_tmrt * pa * pa * pa +                 &
+       ( 3.02122035e-04 )  * d_tmrt * d_tmrt * pa * pa * pa +                  &
+       ( -4.77403547e-06 ) * ta * d_tmrt * d_tmrt * pa * pa * pa +             &
+       ( 1.73825715e-06 )  * va * d_tmrt * d_tmrt * pa * pa * pa +             &
+       ( -4.09087898e-07 ) * d_tmrt * d_tmrt * d_tmrt * pa * pa * pa +         &
+       ( 6.14155345e-01 )  * pa * pa * pa * pa +                               &
+       ( -6.16755931e-02 ) * ta * pa * pa * pa * pa +                          &
+       ( 1.33374846e-03 )  * ta * ta * pa * pa * pa * pa +                     &
+       ( 3.55375387e-03 )  * va * pa * pa * pa * pa +                          &
+       ( -5.13027851e-04 ) * ta * va * pa * pa * pa * pa +                     &
+       ( 1.02449757e-04 )  * va * va * pa * pa * pa * pa +                     &
+       ( -1.48526421e-03 ) * d_tmrt * pa * pa * pa * pa +                      &
+       ( -4.11469183e-05 ) * ta * d_tmrt * pa * pa * pa * pa +                 &
+       ( -6.80434415e-06 ) * va * d_tmrt * pa * pa * pa * pa +                 &
+       ( -9.77675906e-06 ) * d_tmrt * d_tmrt * pa * pa * pa * pa +             &
+       ( 8.82773108e-02 )  * pa * pa * pa * pa * pa +                          &
+       ( -3.01859306e-03 ) * ta * pa * pa * pa * pa * pa +                     &
+       ( 1.04452989e-03 )  * va * pa * pa * pa * pa * pa +                     &
+       ( 2.47090539e-04 )  * d_tmrt * pa * pa * pa * pa * pa +                 &
+       ( 1.48348065e-03 )  * pa * pa * pa * pa * pa * pa 
+
+!-- Consider offset in result
+    utci = utci + offset
+
+ END SUBROUTINE calculate_utci_static
+
+
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> calculate_perct_static: Estimation of perceived temperature (PT, degC)
+!> Value of perct is the Perceived Temperature, degree centigrade
+!------------------------------------------------------------------------------!
+ SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct )
+
+    IMPLICIT NONE
+
+!-- Type of input of the argument list
+    REAL(wp), INTENT ( IN )  :: ta   !< Local air temperature (degC)
+    REAL(wp), INTENT ( IN )  :: vp   !< Local vapour pressure (hPa)
+    REAL(wp), INTENT ( IN )  :: tmrt !< Local mean radiant temperature (degC)
+    REAL(wp), INTENT ( IN )  :: ws   !< Local wind velocitry (m/s)
+    REAL(wp), INTENT ( IN )  :: pair !< Local barometric air pressure (hPa)
+
+!-- Type of output of the argument list
+    REAL(wp), INTENT ( OUT ) :: perct   !< Perceived temperature (degC)
+    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) ::  sclo           !< summer clothing insulation
+    REAL(wp) ::  wclo           !< winter clothing insulation
+    REAL(wp) ::  d_pmv          !< PMV deviation (dimensionless --> PMV)
+    REAL(wp) ::  svp_ta         !< saturation vapor pressure    (hPa)
+    REAL(wp) ::  sult_lim       !< threshold for sultrieness    (hPa)
+    REAL(wp) ::  dgtcm          !< Mean deviation dependent on perct
+    REAL(wp) ::  dgtcstd        !< Mean deviation plus its standard deviation
+    REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
+    REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
+!     REAL(wp) ::  clo_fanger     !< clo for fanger subroutine, unused
+    REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
+    REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
+    REAL(wp) ::  pmva           !< adjusted predicted mean vote
+    REAL(wp) ::  ptc            !< perceived temp. for cold conditions (°C)
+    REAL(wp) ::  d_std          !< factor to threshold for sultriness
+    REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
+    REAL(wp) ::  top            !< Gagge's operative temperatures (°C)
+
+    INTEGER(iwp) :: ncount      !< running index
+    INTEGER(iwp) :: nerr_cold   !< error number (cold conditions)
+    INTEGER(iwp) :: nerr        !< error number
+
+    LOGICAL :: sultrieness
+
+!-- Initialise
+    perct = 9999.0_wp
+
+    nerr     = 0_iwp
+    ncount   = 0_iwp
+    sultrieness  = .FALSE.
+!-- Tresholds: clothing insulation (account for model inaccuracies)
+!   summer clothing
+    sclo     = 0.44453_wp
+!   winter clothing
+    wclo     = 1.76267_wp
+
+!-- decision: firstly calculate for winter or summer clothing
+    IF ( ta <= 10._wp ) THEN
+
+!--    First guess: winter clothing insulation: cold stress
+       clo = wclo
+       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
+       pmv_w = pmva
+
+       IF ( pmva > 0._wp ) THEN
+
+!--       Case summer clothing insulation: heat load ?
+          clo = sclo
+          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
+             top )
+          pmv_s = 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, ws, pair, actlev, eta, sclo,      &
+                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
+             IF ( ncount < 0_iwp ) THEN
+                nerr = -1_iwp
+                RETURN
+             ENDIF
+          ELSE IF ( pmva > 0.06_wp ) THEN
+             clo = 0.5_wp
+             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
+                           pmva, top )
+          ENDIF
+       ELSE IF ( pmva < -0.11_wp ) THEN
+          clo = 1.75_wp
+          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
+             top )
+       ENDIF
+    ELSE
+
+!--    First guess: summer clothing insulation: heat load
+       clo = sclo
+       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
+       pmv_s = pmva
+
+       IF ( pmva < 0._wp ) THEN
+
+!--       Case winter clothing insulation: cold stress ?
+          clo = wclo
+          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
+             top )
+          pmv_w = 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, ws, pair, actlev, eta, sclo,      &
+                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
+             IF ( ncount < 0_iwp ) THEN
+                nerr = -1_iwp
+                RETURN
+             ENDIF
+          ELSE IF ( pmva < -0.11_wp ) THEN
+             clo = 1.75_wp
+             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
+                           pmva, top )
+          ENDIF
+       ELSE IF ( pmva > 0.06_wp ) THEN
+          clo = 0.5_wp
+          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
+             top )
+       ENDIF
+
+    ENDIF
+
+!-- Determine perceived temperature by regression equation + adjustments
+    pmvs = pmva
+    CALL perct_regression ( pmva, clo, perct )
+    ptc = perct
+    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
+!--    Adjust for cold conditions according to Gagge 1986
+       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
+       IF ( nerr_cold > 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 perct_regression ( pmvs, clo, perct )
+    ENDIF
+!     clo_fanger = clo
+    clon = clo
+    IF ( clo > 0.5_wp .AND. perct <= 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 ( perct, 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_ta )
+       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
+       pmvs   = pmva + d_pmv
+       CALL perct_regression ( pmvs, clo, perct )
+       IF ( sult_lim < 99._wp ) THEN
+          IF ( (perct - ptc) > sult_lim ) sultrieness = .TRUE.
+!--       Set factor to threshold for sultriness
+          IF ( dgtcstd /= 0_iwp ) THEN
+             d_std = ( ( perct - ptc ) - dgtcm ) / dgtcstd
+          ENDIF
+       ENDIF
+    ENDIF
+
+ END SUBROUTINE calculate_perct_static
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> The SUBROUTINE calculates the saturation water vapour pressure 
+!> (hPa = hecto Pascal) for a given temperature ta (degC).
+!> For example, ta can be the air temperature or the dew point temperature.
+!------------------------------------------------------------------------------!
+ SUBROUTINE saturation_vapor_pressure( ta, svp_ta )
+
+    IMPLICIT NONE
+
+    REAL(wp), INTENT ( IN )  ::  ta     !< ambient air temperature (degC)
+    REAL(wp), INTENT ( OUT ) ::  svp_ta !< saturation water vapour pressure (hPa)
+
+    REAL(wp)      ::  b
+    REAL(wp)      ::  c
+
+
+    IF ( ta < 0._wp ) THEN
+!--    ta  < 0 (degC): saturation water vapour pressure over ice
+       b = 17.84362_wp
+       c = 245.425_wp
+    ELSE
+!--    ta >= 0 (degC): saturation water vapour pressure over water
+       b = 17.08085_wp
+       c = 234.175_wp
+    ENDIF
+
+!-- Saturation water vapour pressure
+    svp_ta = 6.1078_wp * EXP ( b * ta / ( c + ta ) )
+
+ 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 
+!> (ta,tmrt, vp, ws, pair) and values of individual's activity level
+!------------------------------------------------------------------------------!
+ SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo,             &
+                       pmv_s, wclo, pmv_w, eps, pmva, top, nerr,               &
+                       clo_res )
+
+    IMPLICIT NONE
+
+!-- Input variables of argument list:
+    REAL(wp), INTENT ( IN )  :: ta     !< Ambient temperature (degC)
+    REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC)
+    REAL(wp), INTENT ( IN )  :: vp     !< Water vapour pressure (hPa)
+    REAL(wp), INTENT ( IN )  :: ws    !< Wind speed (m/s) 1 m above ground
+    REAL(wp), INTENT ( IN )  :: pair       !< 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 )  :: sclo     !< Lower threshold of bracketing clothing insulation (clo) 
+    REAL(wp), INTENT ( IN )  :: wclo     !< 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_w    !< Fanger's PMV corresponding to wclo
+    REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
+
+! 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 (degC) 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 sclo and wclo
+
+!-- Type of program variables
+    INTEGER(iwp), PARAMETER  ::  max_iteration = 15_iwp       !< max number of iterations
+    REAL(wp),     PARAMETER  ::  guess_0       = -1.11e30_wp  !< initial guess
+    REAL(wp) ::  x_ridder    !< current guess for clothing insulation   (clo)
+    REAL(wp) ::  clo_lower   !< lower limit of clothing insulation      (clo)
+    REAL(wp) ::  clo_upper   !< upper limit of clothing insulation      (clo)
+    REAL(wp) ::  x_lower     !< lower guess for clothing insulation     (clo)
+    REAL(wp) ::  x_upper     !< upper guess for clothing insulation     (clo)
+    REAL(wp) ::  x_average   !< average of x_lower and x_upper          (clo)
+    REAL(wp) ::  x_new       !< preliminary result for clothing insulation (clo)
+    REAL(wp) ::  y_lower     !< predicted mean vote for summer clothing
+    REAL(wp) ::  y_upper     !< predicted mean vote for winter clothing
+    REAL(wp) ::  y_average   !< average of y_lower and y_upper
+    REAL(wp) ::  y_new       !< preliminary result for pred. mean vote
+    REAL(wp) ::  sroot       !< sqrt of PMV-guess
+    INTEGER(iwp) ::  j       !< running index
+
+!-- Initialise
+    nerr    = 0_iwp
+
+!-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation
+    pmva        = 0._wp
+    clo_lower   = sclo
+    y_lower     = pmv_s
+    clo_upper   = wclo
+    y_upper     = pmv_w
+    IF ( ( y_lower > 0._wp .AND. y_upper < 0._wp ) .OR.                        &
+         ( y_lower < 0._wp .AND. y_upper > 0._wp ) ) THEN
+       x_lower  = clo_lower
+       x_upper  = clo_upper
+       x_ridder = guess_0
+
+       DO j = 1_iwp, max_iteration
+          x_average = 0.5_wp * ( x_lower + x_upper )
+          CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta,        &
+                        y_average, top )
+          sroot = SQRT ( y_average**2 - y_lower * y_upper )
+          IF ( sroot == 0._wp ) THEN
+             clo_res = x_average
+             nerr = j
+             RETURN
+          ENDIF
+          x_new = x_average + ( x_average - x_lower ) *                        &
+                      ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot )
+          IF ( ABS ( x_new - x_ridder ) <= eps ) THEN
+             clo_res = x_ridder
+             nerr       = j
+             RETURN
+          ENDIF
+          x_ridder = x_new
+          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
+                        y_new, top )
+          IF ( y_new == 0._wp ) THEN
+             clo_res = x_ridder
+             nerr       = j
+             RETURN
+          ENDIF
+          IF ( SIGN ( y_average, y_new ) /= y_average ) THEN
+             x_lower = x_average
+             y_lower = y_average
+             x_upper  = x_ridder
+             y_upper  = y_new
+          ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower ) THEN
+             x_upper  = x_ridder
+             y_upper  = y_new
+          ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper ) THEN
+             x_lower = x_ridder
+             y_lower = y_new
+          ELSE
+!--          Never get here in x_ridder: singularity in y
+             nerr       = -1_iwp
+             clo_res = x_ridder
+             RETURN
+          ENDIF
+          IF ( ABS ( x_upper - x_lower ) <= eps ) THEN
+             clo_res = x_ridder
+             nerr       = j
+             RETURN
+          ENDIF
+       ENDDO
+!--    x_ridder exceed maximum iterations
+       nerr       = -2_iwp
+       clo_res = y_new
+       RETURN
+    ELSE IF ( y_lower == 0. ) THEN
+       x_ridder = clo_lower
+    ELSE IF ( y_upper == 0. ) THEN
+       x_ridder = clo_upper
+    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 (perct) and (adjusted) 
+!> PMV. The regression presumes the Klima-Michel settings for reference 
+!> individual and reference environment.
+!------------------------------------------------------------------------------!
+ SUBROUTINE perct_regression( pmv, clo, perct )
+
+    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 ) ::  perct   !< perct (degC) corresponding to given PMV / clo
+
+    IF ( pmv <= -0.11_wp ) THEN
+       perct = 5.805_wp + 12.6784_wp * pmv
+    ELSE
+       IF ( pmv >= + 0.01_wp ) THEN
+          perct = 16.826_wp + 6.163_wp * pmv
+       ELSE
+          perct = 21.258_wp - 9.558_wp * clo
+       ENDIF
+    ENDIF
+
+ END SUBROUTINE perct_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 (ws < 0.1 m/s) is dealt with ws = 0.1 m/s
+!------------------------------------------------------------------------------!
+ SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva, top )
+
+    IMPLICIT NONE
+
+!-- Input variables of argument list:
+    REAL(wp), INTENT ( IN ) ::  ta       !< Ambient air temperature (degC)
+    REAL(wp), INTENT ( IN ) ::  tmrt     !< Mean radiant temperature (degC)
+    REAL(wp), INTENT ( IN ) ::  pa       !< Water vapour pressure (hPa)
+    REAL(wp), INTENT ( IN ) ::  pair     !< Barometric pressure (hPa) at site
+    REAL(wp), INTENT ( IN ) ::  in_ws    !< 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 (PMV, 
+!            dimensionless) according to Fanger corresponding to meteorological 
+!            (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
+    REAL(wp), INTENT ( OUT ) ::  top     !< operative temperature (degC)
+
+!-- Internal variables
+    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
+    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
+    REAL(wp) ::  activity     !< persons activity  (must stay == actlev, W)
+    REAL(wp) ::  t_skin_aver  !< average skin temperature               (°C)
+    REAL(wp) ::  bc           !< preliminary result storage
+    REAL(wp) ::  cc           !< preliminary result storage
+    REAL(wp) ::  dc           !< preliminary result storage
+    REAL(wp) ::  ec           !< preliminary result storage
+    REAL(wp) ::  gc           !< preliminary result storage
+    REAL(wp) ::  t_clothing   !< clothing temperature                   (°C)
+    REAL(wp) ::  hr           !< radiational heat resistence
+    REAL(wp) ::  clo          !< clothing insulation index              (clo)
+    REAL(wp) ::  ws           !< wind speed                             (m/s)
+    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat 
+!             ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
+    REAL(wp) ::  z2           !< Water vapour diffution through the skin
+    REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
+    REAL(wp) ::  z4           !< Loss of latent heat through respiration
+    REAL(wp) ::  z5           !< Loss of radiational heat
+    REAL(wp) ::  z6           !< Heat loss through forced convection
+    INTEGER(iwp) :: i         !< running index
+
+!-- 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 (ws < 0.1 m/s ) not considered
+    ws = in_ws
+    IF ( ws < .1_wp ) THEN
+       ws = .1_wp
+    ENDIF
+
+!-- Heat_convection = forced convection
+    heat_convection = 12.1_wp * SQRT ( ws * pair / 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 + degc_to_k )**4 + ec * cc * ta ) / bc
+
+!-- Calculation of clothing surface temperature (t_clothing) based on
+!   Newton-approximation with air temperature as initial guess
+    t_clothing = ta
+    DO i = 1, 3
+       t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing  &
+          * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**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 - ta )
+
+!-- Loss of radiational heat
+    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
+       degc_to_k )**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 - ta )
+
+!-- Predicted Mean Vote
+    pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
+
+!-- Operative temperatur
+    top = ( hr * tmrt + heat_convection * ta ) / ( 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, ta, vp, svp_ta, tmrt, ws, 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 )  :: ta       !< Ambient temperature (degC) at screen level
+    REAL(wp),     INTENT ( IN )  :: vp       !< Water vapour pressure (hPa) at screen level
+    REAL(wp),     INTENT ( IN )  :: svp_ta   !< Saturation water vapour pressure (hPa) at ta
+    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
+    REAL(wp),     INTENT ( IN )  :: ws       !< 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) ::  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_ta
+    p95  = 1.00_wp * svp_ta
+    IF ( vp >= p10 .AND. vp <= p95 ) THEN
+       pa = vp
+    ELSE
+       nerr = -3_iwp
+       IF ( vp < p10 ) THEN
+!--       Due to conditions of regression: r.H. >= 5 %
+          pa = p10
+       ELSE
+!--       Due to conditions of regression: r.H. <= 95 %
+          pa = p95
+       ENDIF
+    ENDIF
+    IF ( pa > 0._wp ) THEN
+!--    Natural logarithm of pa
+       apa = LOG ( pa )
+    ELSE
+       apa = -5._wp
+    ENDIF
+
+!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
+    pa_p50   = 0.5_wp * svp_ta
+    IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN
+       dapa   = apa - LOG ( pa_p50 )
+       pa_p50 = pa / pa_p50
+    ELSE
+       dapa   = -5._wp
+       pa_p50 = 0._wp
+    ENDIF
+!-- Square root of wind velocity
+    IF ( ws >= 0._wp ) THEN
+       sqvel = SQRT ( ws )
+    ELSE
+       sqvel = 0._wp
+    ENDIF
+!-- Air temperature
+!    ta = ta
+!-- Difference mean radiation to air temperature
+    dtmrt = tmrt - ta
+
+!-- Select the valid regression coefficients
+    nreg = INT ( pmv )
+    IF ( nreg < 0_iwp ) THEN
+!--    value of the FUNCTION in the case pmv <= -1
+       deltapmv = 0._wp
+       RETURN
+    ENDIF
+    gew = MOD ( pmv, 1._wp )
+    IF ( gew < 0._wp ) gew = 0._wp
+    IF ( nreg > 5_iwp ) THEN
+       ! nreg=6
+       nreg  = 5_iwp
+       gew   = pmv - 5._wp
+       gew2  = pmv - 6._wp
+       IF ( gew2 > 0_iwp ) THEN
+          gew = ( gew - gew2 ) / gew
+       ENDIF
+    ENDIF
+
+!-- Regression valid for 0. <= pmv <= 6.
+    dpmv_1 =                                                                   &
+       + bpa ( nreg ) * pa                                                     &
+       + bpmv ( nreg ) * pmv                                                   &
+       + bapa ( nreg ) * apa                                                   &
+       + bta ( nreg ) * ta                                                     &
+       + bdtmrt ( nreg ) * dtmrt                                               &
+       + bdapa ( nreg ) * dapa                                                 &
+       + bsqvel ( nreg ) * sqvel                                               &
+       + bpa_p50 ( nreg ) * pa_p50                                             &
+       + aconst ( nreg )
+
+    dpmv_2 = 0._wp
+    IF ( nreg < 6_iwp ) THEN
+       dpmv_2 =                                                                &
+          + bpa ( nreg + 1_iwp )     * pa                                      &
+          + bpmv ( nreg + 1_iwp )    * pmv                                     &
+          + bapa ( nreg + 1_iwp )    * apa                                     &
+          + bta ( nreg + 1_iwp )     * ta                                      &
+          + bdtmrt ( nreg + 1_iwp )  * dtmrt                                   &
+          + bdapa ( nreg + 1_iwp )   * dapa                                    &
+          + bsqvel ( nreg + 1_iwp )  * sqvel                                   &
+          + bpa_p50 ( nreg + 1_iwp ) * pa_p50                                  &
+          + aconst ( nreg + 1_iwp )
+    ENDIF
+
+!-- Calculate pmv modification
+    deltapmv = ( 1._wp - gew ) * dpmv_1 + gew * dpmv_2
+    pmvs = pmva + deltapmv
+    IF ( ( pmvs ) < 0._wp ) THEN
+!--    Prevent negative pmv* due to problems with clothing insulation
+       nerr = -4_iwp
+       IF ( pmvs > -0.11_wp ) THEN
+!--       Threshold from FUNCTION perct_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 
+!> perct. 
+!------------------------------------------------------------------------------!
+ SUBROUTINE calc_sultr( perct, dperctm, dperctstd, sultr_res )
+
+    IMPLICIT NONE
+
+!-- Input of the argument list:
+    REAL(wp), INTENT ( IN )  ::  perct      !< Classical perceived temperature: Base is Fanger's PMV
+
+!-- Additional output variables of argument list:
+    REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational gt 
+!             calculated based on Gagge's rational PMV*)
+    REAL(wp), INTENT ( OUT ) ::  dperctstd  !< dperctm 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 ::  dperctka = +7.5776086_wp
+    REAL(wp), PARAMETER ::  dperctkb = -0.740603_wp
+    REAL(wp), PARAMETER ::  dperctkc = +0.0213324_wp
+    REAL(wp), PARAMETER ::  dperctkd = -0.00027797237_wp
+
+!-- Types of coefficients mean deviation plus standard deviation
+!   regression coefficients: third order polynomial
+    REAL(wp), PARAMETER ::  dperctsa = +0.0268918_wp
+    REAL(wp), PARAMETER ::  dperctsb = +0.0465957_wp
+    REAL(wp), PARAMETER ::  dperctsc = -0.00054709752_wp
+    REAL(wp), PARAMETER ::  dperctsd = +0.0000063714823_wp
+
+!-- Factor to mean standard deviation defining SIGNificance for 
+!   sultriness
+    REAL(wp), PARAMETER :: faktor = 1._wp
+
+!-- Initialise
+    sultr_res = +99._wp
+    dperctm   = 0._wp
+    dperctstd = 999999._wp
+
+    IF ( perct < 16.826_wp .OR. perct > 56._wp ) THEN
+!--    Unallowed classical PMV/perct
+       RETURN
+    ENDIF
+
+!-- Mean deviation dependent on perct
+    dperctm = dperctka + dperctkb * perct + dperctkc * perct**2._wp + dperctkd * perct**3._wp
+
+!-- Mean deviation plus its standard deviation
+    dperctstd = dperctsa + dperctsb * perct + dperctsc * perct**2._wp + dperctsd * perct**3._wp
+
+!-- Value of the FUNCTION
+    sultr_res = dperctm + faktor * dperctstd
+    IF ( ABS ( sultr_res ) > 99._wp ) sultr_res = +99._wp
+
+ END SUBROUTINE calc_sultr
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Multiple linear regression to calculate an increment delta_cold, 
+!> 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, ta, ws, 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 ) ::  ta     !< Air temperature 2 m above ground (degC)
+    REAL(wp), INTENT ( IN ) ::  ws     !< Relative wind velocity 1 m above ground (m/s)
+    REAL(wp), INTENT ( IN ) ::  tmrt   !< Mean radiant temperature (degC)
+
+!-- 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_cold(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_ws
+    INTEGER(iwp) ::  i
+    INTEGER(iwp) ::  j
+    INTEGER(iwp) ::  i_bin
+
+!-- Coefficient of the 3 regression lines
+    !     1:const   2:*pmvc  3:*ta      4:*sqrt_ws  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 - ta
+    sqrt_ws   = ws
+    IF ( sqrt_ws < 0.10_wp ) THEN
+       sqrt_ws = 0.10_wp
+    ELSE
+       sqrt_ws = SQRT ( sqrt_ws )
+    ENDIF
+
+    DO i = 1, 3
+       delta_cold (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) * ta + coeff(i, 4) *             &
+                  sqrt_ws + coeff(i,5)*dtmrt 
+       delta_cold(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_cold(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 (perct) 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 perct < 10 (degC)
+!------------------------------------------------------------------------------!
+ REAL(wp) FUNCTION ireq_neutral( perct, ireq_minimal, nerr )
+
+    IMPLICIT NONE
+
+!-- Type declaration of arguments
+    REAL(wp),     INTENT ( IN )  ::  perct
+    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 * perct
+
+!-- IREQ neutral conditions (thermal comfort)
+    ireq_neutral = 1.62_wp - 0.0564_wp * top02
+
+!-- Regression only defined for perct <= 10 (degC)
+    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
+
+
+!------------------------------------------------------------------------------!
+! 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) ::  energy_prod
+    REAL(wp) ::  s
+    REAL(wp) ::  factor
+    REAL(wp) ::  basic_heat_prod
+
+    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 )
+    basic_heat_prod = 0.
+    IF ( sex == 1_iwp ) THEN
+       basic_heat_prod = 3.45_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
+                     .01_wp  * ( s - 43.4_wp ) )
+    ELSE IF ( sex == 2_iwp ) THEN
+       basic_heat_prod = 3.19_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
+                    .018_wp * ( s - 42.1_wp ) )
+    END IF
+
+    energy_prod = work + basic_heat_prod
+    actlev = energy_prod / 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, ws, tmrt, pair, 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) ::  ws         !< Wind speed in approx. 1.1m (m/s)
+    REAL(wp), INTENT(in) ::  tmrt       !< Mean radiant temperature   (°C)
+    REAL(wp), INTENT(in) ::  pair       !< 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) ::  sclo
+    REAL(wp) ::  wclo
+    REAL(wp) ::  d_pmv
+    REAL(wp) ::  svp_ta
+    REAL(wp) ::  sult_lim
+    REAL(wp) ::  dgtcm
+    REAL(wp) ::  dgtcstd
+    REAL(wp) ::  clon
+    REAL(wp) ::  ireq_minimal
+!     REAL(wp) ::  clo_fanger
+    REAL(wp) ::  pmv_w
+    REAL(wp) ::  pmv_s
+    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) ::  ncount
+    INTEGER(iwp) ::  nerr_cold
+    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
+    ncount     = 0_wp
+    sultrieness    = .FALSE.
+!-- Tresholds: clothing insulation (account for model inaccuracies)
+!-- Summer clothing
+    sclo     = 0.44453_wp
+!-- Winter clothing
+    wclo     = 1.76267_wp
+
+!-- Decision: firstly calculate for winter or summer clothing
+    IF ( ta <= 10._wp ) THEN
+
+!--    First guess: winter clothing insulation: cold stress
+       clo = wclo
+       t_clothing = -999.0_wp  ! force initial run
+       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
+          t_clothing, storage, dt, pmva )
+       pmv_w = pmva
+
+       IF ( pmva > 0._wp ) THEN
+
+!--       Case summer clothing insulation: heat load ?            
+          clo = sclo
+          t_clothing = -999.0_wp  ! force initial run
+          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
+             t_clothing, storage, dt, pmva )
+          pmv_s = 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, ws, pair, actlev, eta , sclo,     &
+                            pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
+             IF ( ncount < 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, ws, pair, 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, ws, pair, clo, actlev, work,      &
+             t_clothing, storage, dt, pmva )
+       ENDIF
+
+    ELSE
+
+!--    First guess: summer clothing insulation: heat load
+       clo = sclo
+       t_clothing = -999.0_wp
+       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
+          t_clothing, storage, dt, pmva )
+       pmv_s = pmva
+
+       IF ( pmva < 0._wp ) THEN
+
+!--       Case winter clothing insulation: cold stress ?
+          clo = wclo
+          t_clothing = -999.0_wp
+          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
+             t_clothing, storage, dt, pmva )
+          pmv_w = 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, ws, pair, actlev, eta, sclo,      &
+                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
+             IF ( ncount < 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, ws, pair, 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, ws, pair, clo, actlev, work,      &
+             t_clothing, storage, dt, pmva )
+       ENDIF
+
+    ENDIF
+
+!-- Determine perceived temperature by regression equation + adjustments
+    pmvs = pmva
+    CALL perct_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, ws, tmrt, nerr_cold, d_pmv )
+       IF ( nerr_cold > 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 perct_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_ta )
+       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
+       pmvs   = pmva + d_pmv
+       CALL perct_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, ws, tmrt, pair, 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 )  ::  ws      !< Wind speed                  (m/s)
+    REAL(wp), INTENT ( IN )  ::  pair    !< 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 internal variables
+    REAL(wp) ::  d_pmv
+    REAL(wp) ::  svp_ta
+    REAL(wp) ::  sult_lim
+    REAL(wp) ::  dgtcm
+    REAL(wp) ::  dgtcstd
+    REAL(wp) ::  pmva
+    REAL(wp) ::  ptc
+    REAL(wp) ::  d_std
+    REAL(wp) ::  pmvs
+    INTEGER(iwp) ::  nerr_cold
+    INTEGER(iwp) ::  nerr
+
+    LOGICAL ::  sultrieness
+
+!-- Initialise
+    ipt = -999.0_wp
+
+    nerr     = 0_iwp
+    sultrieness  = .FALSE.
+
+!-- Determine pmv_adjusted for current conditions
+    CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,            &
+       t_clothing, storage, dt, pmva )
+
+!-- Determine perceived temperature by regression equation + adjustments
+    CALL perct_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, ws, tmrt, nerr_cold, d_pmv )
+       IF ( nerr_cold > 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 perct_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_ta )
+       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
+       pmvs   = pmva + d_pmv
+       CALL perct_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 (ta,tmrt,pa,ws,pair) 
+!> and individual variables (clo, actlev, eta) considering a storage
+!> and clothing temperature for a given timestep.
+!------------------------------------------------------------------------------!
+ SUBROUTINE fanger_s_acti( ta, tmrt, pa, in_ws, pair, in_clo, actlev,          &
+    activity, t_cloth, s, dt, pmva )
+
+    IMPLICIT NONE
+
+!--  Input argument types
+    REAL(wp), INTENT ( IN )  ::  ta       !< Air temperature          (°C)
+    REAL(wp), INTENT ( IN )  ::  tmrt     !< Mean radiant temperature (°C)
+    REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
+    REAL(wp), INTENT ( IN )  ::  pair     !< Air pressure             (hPa)
+    REAL(wp), INTENT ( IN )  ::  in_ws   !< 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), PARAMETER  ::  time_equil = 7200._wp
+
+    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
+    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
+    REAL(wp) ::  t_skin_aver  !< average skin temperature               (°C)
+    REAL(wp) ::  bc           !< preliminary result storage
+    REAL(wp) ::  cc           !< preliminary result storage
+    REAL(wp) ::  dc           !< preliminary result storage
+    REAL(wp) ::  ec           !< preliminary result storage
+    REAL(wp) ::  gc           !< preliminary result storage
+    REAL(wp) ::  t_clothing   !< clothing temperature                   (°C)
+    REAL(wp) ::  hr           !< radiational heat resistence
+    REAL(wp) ::  clo          !< clothing insulation index              (clo)
+    REAL(wp) ::  ws           !< wind speed                             (m/s)
+    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat 
+!             ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
+    REAL(wp) ::  z2           !< Water vapour diffution through the skin
+    REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
+    REAL(wp) ::  z4           !< Loss of latent heat through respiration
+    REAL(wp) ::  z5           !< Loss of radiational heat
+    REAL(wp) ::  z6           !< Heat loss through forced convection
+    REAL(wp) ::  en           !< Energy ballance                        (W)
+    REAL(wp) ::  d_s          !< Storage delta                          (W)
+    REAL(wp) ::  adjustrate   !< Max storage adjustment rate
+    REAL(wp) ::  adjustrate_cloth  !< max clothing temp. adjustment rate
+
+    INTEGER(iwp) :: i         !< running index
+    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 (ws < 0.1 m/s ) not considered
+    ws = in_ws
+    IF ( ws < .1_wp ) THEN
+       ws = .1_wp
+    ENDIF
+
+!-- Heat_convection = forced convection
+    heat_convection = 12.1_wp * SQRT ( ws * pair / 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 * ta ) / 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 = ta
+    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 - ta )
+
+!-- 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 - ta )
+
+!-- 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
+
+
+
+!------------------------------------------------------------------------------!
+!
+! 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, pair, pet )
+
+    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 ) ::  pair  !< Air pressure                (hPa)
+
+!-- Output arguments:
+    REAL(wp), INTENT ( OUT ) ::  pet  !< PET                         (°C)
+
+!-- Internal variables:
+    REAL(wp) ::  acl        !< clothing area                        (m²)
+    REAL(wp) ::  adu        !< Du Bois area                         (m²)
+    REAL(wp) ::  aeff       !< effective area                       (m²)
+    REAL(wp) ::  ere        !< energy ballance                      (W)
+    REAL(wp) ::  erel       !< latent energy ballance               (W)
+    REAL(wp) ::  esw        !< Energy-loss through sweat evap.      (W)
+    REAL(wp) ::  facl       !< Surface area extension through clothing (factor)
+    REAL(wp) ::  feff       !< Surface modification by posture      (factor)
+    REAL(wp) ::  rdcl       !< Diffusion resistence of clothing     (factor)
+    REAL(wp) ::  rdsk       !< Diffusion resistence of skin         (factor)
+    REAL(wp) ::  rtv
+    REAL(wp) ::  vpts       !< Sat. vapor pressure over skin        (hPa)
+    REAL(wp) ::  tsk        !< Skin temperature                     (°C)
+    REAL(wp) ::  tcl        !< Clothing temperature                 (°C)
+    REAL(wp) ::  wetsk      !< Fraction of wet skin                 (factor)
+
+!-- Variables:
+    REAL(wp) :: int_heat    !< 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) :: clo         !< 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
+
+!-- Configuration, keep standard parameters!
+    age   = 35._wp
+    mbody = 75._wp
+    ht    =  1.75_wp
+    work  = 80._wp
+    eta   =  0._wp
+    clo   =  0.9_wp
+    fcl   =  1.15_wp
+
+!-- Call subfunctions
+    CALL in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,    &
+            vpa, work )
+
+    CALL heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, &
+            int_heat,mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,      &
+            vpts, wetsk )
+
+    CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair, rdcl,&
+            rdsk, rtv, ta, tcl, tsk, pet, vpts, wetsk )
+
+
+ END SUBROUTINE calculate_pet_static
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate internal energy ballance
+!------------------------------------------------------------------------------!
+ SUBROUTINE in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
+    vpa, work )
+
+!-- Input arguments:
+    REAL(wp), INTENT( IN )  ::  pair      !< air pressure             (hPa)
+    REAL(wp), INTENT( IN )  ::  ta        !< air temperature          (°C)
+    REAL(wp), INTENT( IN )  ::  vpa       !< vapor pressure           (hPa)
+    REAL(wp), INTENT( IN )  ::  age       !< Persons age              (a)
+    REAL(wp), INTENT( IN )  ::  mbody     !< Persons body mass        (kg)
+    REAL(wp), INTENT( IN )  ::  ht        !< Persons height           (m)
+    REAL(wp), INTENT( IN )  ::  work      !< Work load                (W)
+    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency      (dimensionless)
+
+!-- Output arguments:
+    REAL(wp), INTENT( OUT ) ::  ere       !< energy ballance          (W)
+    REAL(wp), INTENT( OUT ) ::  erel      !< latent energy ballance   (W)
+    REAL(wp), INTENT( OUT ) ::  int_heat  !< internal heat production (W)
+    REAL(wp), INTENT( OUT ) ::  rtv       !< respiratory volume
+
+!-- Constants:
+!     REAL(wp), PARAMETER :: cair = 1010._wp        !< replaced by c_p
+!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
+
+!-- Internal variables:
+    REAL(wp) ::  eres                     !< Sensible respiratory heat flux (W)
+    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
+
+    int_heat = 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 = c_p * (ta - tex) * rtv
+
+!-- Latent respiration energy
+    vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
+    erel = 0.623_wp * l_v / pair * ( 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, clo, ere, erel, esw, facl, fcl, feff,  &
+        ht, int_heat, mbody, pair, 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 )  ::  int_heat  !< internal heat production (W)
+    REAL(wp), INTENT( IN )  ::  pair   !< 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)
+    REAL(wp), INTENT( IN )  ::  mbody  !< body mass                (kg)
+    REAL(wp), INTENT( IN )  ::  ht     !< height                   (m)
+    REAL(wp), INTENT( IN )  ::  clo    !< clothing insulation      (clo)
+    REAL(wp), INTENT( IN )  ::  fcl    !< factor for surface area increase by clothing
+
+!-- 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)
+
+!-- Cconstants:
+!     REAL(wp), PARAMETER :: cair = 1010._wp      !< replaced by c_p
+    REAL(wp), PARAMETER :: cb   = 3640._wp        !< 
+    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
+    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
+!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
+    REAL(wp), PARAMETER :: food =    0._wp        !< Heat gain by food        (W)
+    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
+    REAL(wp), PARAMETER :: rob  =    1.06_wp      !< 
+
+!-- Internal variables
+    REAL(wp) ::  c(0:10)        !< Core temperature array           (°C)
+    REAL(wp) ::  cbare          !< Convection through bare skin
+    REAL(wp) ::  cclo           !< Convection through clothing
+    REAL(wp) ::  csum           !< Convection in total
+    REAL(wp) ::  di             !< difference between r1 and r2
+    REAL(wp) ::  ed             !< energy transfer by diffusion     (W)
+    REAL(wp) ::  enbal          !< energy ballance                  (W)
+    REAL(wp) ::  enbal2         !< energy ballance (storage, last cycle)
+    REAL(wp) ::  eswdif         !< difference between sweat production and evaporation potential
+    REAL(wp) ::  eswphy         !< sweat created by physiology
+    REAL(wp) ::  eswpot         !< potential sweat evaporation
+    REAL(wp) ::  fec            !< 
+    REAL(wp) ::  hc             !< 
+    REAL(wp) ::  he             !< 
+    REAL(wp) ::  htcl           !< 
+    REAL(wp) ::  r1             !< 
+    REAL(wp) ::  r2             !< 
+    REAL(wp) ::  rbare          !< Radiational loss of bare skin    (W/m²)
+    REAL(wp) ::  rcl            !< 
+    REAL(wp) ::  rclo           !< Radiational loss of clothing     (W/m²)
+    REAL(wp) ::  rclo2          !< Longwave radiation gain or loss  (W/m²)
+    REAL(wp) ::  rsum           !< Radiational loss or gain         (W/m²)
+    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             !< modification step                (K)
+    REAL(wp) ::  y              !< fraction of bare skin
+    INTEGER(iwp) ::  count1     !< running index
+    INTEGER(iwp) ::  count3     !< running index
+    INTEGER(iwp) ::  j          !< running index
+    INTEGER(iwp) ::  i          !< running index
+    LOGICAL ::  skipIncreaseCount   !< iteration control flag
+
+    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 * (pair /po) ** 0.55_wp
+    feff = 0.725_wp                     !< Posture: 0.725 for stading
+    facl = (- 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp     &
+          * (clo ** 3._wp)) / 100._wp
+
+    IF ( facl > 1._wp )   facl = 1._wp
+    rcl = ( clo / 6.45_wp ) / facl
+    IF ( clo >= 2._wp )  y  = 1._wp
+
+    IF ( ( clo > 0.6_wp ) .AND. ( clo < 2._wp ) )  y = ( ht - 0.2_wp ) / ht
+    IF ( ( clo <= 0.6_wp ) .AND. ( clo > 0.3_wp ) ) y = 0.5_wp
+    IF ( ( clo <= 0.3_wp ) .AND. ( clo > 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 * sigma_sb * ( ( tcl + degc_to_k )**4._wp -             &
+            ( tmrt + degc_to_k )** 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 * sigma_sb *                  &
+            ( ( tmrt + degc_to_k )** 4._wp - ( tsk + degc_to_k )** 4._wp )
+          rclo  = feff * acl * emcl * sigma_sb *                               &
+            ( ( tmrt + degc_to_k )** 4._wp - ( tcl + degc_to_k )** 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)  = int_heat + 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  !< no need for sweating
+
+          sw = swm
+          eswphy = - sw * l_v
+          he     = 0.633_wp * hc / ( pair * c_p )
+          fec    = 1._wp / ( 1._wp + 0.92_wp * hc * rcl )
+          eswpot = he * ( vpa - vpts ) * adu * l_v * fec
+          wetsk  = eswphy / eswpot
+
+          IF ( wetsk > 1._wp ) wetsk = 1._wp
+!
+!--       Sweat production > evaporation?
+          eswdif = eswphy - eswpot
+
+          IF ( eswdif <= 0._wp ) esw = eswpot  !< Limit is evaporation
+          IF ( eswdif > 0._wp ) esw = eswphy   !< Limit is sweat production
+          IF ( esw  > 0._wp )   esw = 0._wp    !< Sweat can't be evaporated, no more cooling effect
+
+!--       Diffusion
+          rdsk = 0.79_wp * 10._wp ** 7._wp
+          rdcl = 0._wp
+          ed   = l_v / ( 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 = int_heat + 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 / l_v * 3600._wp * ( -1000._wp )
+       wr = erel / l_v * 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, int_heat, pair,   &
+        rdcl, rdsk, rtv, ta, tcl, tsk, pet, 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 ) ::  int_heat  !< internal heat production (W)
+    REAL(wp), INTENT( IN ) ::  pair  !< 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   !< respiratory volume
+    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 ) ::  pet   !< PET                          (°C)
+
+!-- Cconstants:
+!     REAL(wp), PARAMETER :: cair = 1010._wp        !< replaced by c_p
+    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
+    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
+!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
+    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
+
+!-- Internal variables
+    REAL ( wp ) ::  cbare             !< Convection through bare skin
+    REAL ( wp ) ::  cclo              !< Convection through clothing
+    REAL ( wp ) ::  csum              !< Convection in total
+    REAL ( wp ) ::  ed                !< Diffusion                      (W)
+    REAL ( wp ) ::  enbal             !< Energy ballance                (W)
+    REAL ( wp ) ::  enbal2            !< Energy ballance (last iteration cycle)
+    REAL ( wp ) ::  ere               !< Energy ballance result         (W)
+    REAL ( wp ) ::  erel              !< Latent energy ballance         (W)
+    REAL ( wp ) ::  eres              !< Sensible respiratory heat flux (W)
+    REAL ( wp ) ::  hc                !<
+    REAL ( wp ) ::  rbare             !< Radiational loss of bare skin  (W/m²)
+    REAL ( wp ) ::  rclo              !< Radiational loss of clothing   (W/m²)
+    REAL ( wp ) ::  rsum              !< Radiational loss or gain       (W/m²)
+    REAL ( wp ) ::  tex               !< Temperat. of exhaled air       (°C)
+    REAL ( wp ) ::  vpex              !< Vapor pressure of exhaled air  (hPa)
+    REAL ( wp ) ::  xx                !< Delta PET per iteration        (K)
+
+    INTEGER ( iwp ) ::  count1        !< running index
+    INTEGER ( iwp ) ::  i             !< running index
+
+    pet = 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 * ( pair / po ) ** 0.55_wp
+
+!--       Radiation
+          aeff  = adu * feff
+          rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
+              ( ( pet + degc_to_k ) ** 4._wp - ( tsk + degc_to_k ) ** 4._wp )
+          rclo  = feff * acl * emcl * sigma_sb *                               &
+              ( ( pet + degc_to_k ) ** 4._wp - ( tcl + degc_to_k ) ** 4._wp )
+          rsum  = rbare + rclo
+
+!--       Covection
+          cbare = hc * ( pet - tsk ) * adu * ( 1._wp - facl )
+          cclo  = hc * ( pet - tcl ) * acl
+          csum  = cbare + cclo
+
+!--       Diffusion
+          ed = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( 12._wp -    &
+             vpts )
+
+!--       Respiration
+          tex  = 0.47_wp * pet + 21._wp
+          eres = c_p * ( pet - tex ) * rtv
+          vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
+          erel = 0.623_wp * l_v / pair * ( 12._wp - vpex ) * rtv
+          ere  = eres + erel
+
+!--       Energy ballance
+          enbal = int_heat + 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 )  pet = pet - xx
+          IF ( enbal < 0._wp )  pet = pet + 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_mod

palm/trunk/SOURCE/check_parameters.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3517 2018-11-12 16:27:08Z gronemeier
 ! bugfix: renamed 'w*2pt*' -> 'w*2theta*'
 !                 'w*pt*2' -> 'w*theta*2'
@@ -721 +724 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_check_data_output, biom_check_parameters
+        ONLY:  bio_check_data_output, bio_check_parameters
 
     USE bulk_cloud_model_mod,                                                  &
@@ -1507 +1510 @@
 
 !-- Check the module settings
-    IF ( biometeorology )       CALL biom_check_parameters
+    IF ( biometeorology )       CALL bio_check_parameters
     IF ( bulk_cloud_model )     CALL bcm_check_parameters
     IF ( gust_module_enabled )  CALL gust_check_parameters
@@ -3211 +3214 @@
 
              IF ( unit == 'illegal'  .AND.  biometeorology )  THEN
-                CALL biom_check_data_output( var, unit )
+                CALL bio_check_data_output( var, unit )
              ENDIF
 

palm/trunk/SOURCE/data_output_2d.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3467 2018-10-30 19:05:21Z suehring
 ! Implementation of a new aerosol module salsa.
 ! 
@@ -263 +266 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_data_output_2d
+        ONLY:  bio_data_output_2d
 
     USE bulk_cloud_model_mod,                                                  &
@@ -1324 +1327 @@
                 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 )
+                   CALL bio_data_output_2d( av, do2d(av,if), found, grid,      &
+                                            local_pf, two_d, nzb_do, nzt_do,   &
+                                            fill_value )
                 ENDIF
 

palm/trunk/SOURCE/data_output_3d.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3467 2018-10-30 19:05:21Z suehring
 ! Implementation of a new aerosol module salsa.
 ! 
@@ -236 +239 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_data_output_3d
+        ONLY:  bio_data_output_3d
 
     USE bulk_cloud_model_mod,                                                  &
@@ -790 +793 @@
 
              IF ( .NOT. found  .AND.  biometeorology )  THEN
-                 CALL biom_data_output_3d( av, do3d(av,if), found, local_pf,   &
-                                           nzb_do, nzt_do )
+                 CALL bio_data_output_3d( av, do3d(av,if), found, local_pf,    &
+                                          nzb_do, nzt_do )
              ENDIF
 

palm/trunk/SOURCE/header.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3524 2018-11-14 13:36:44Z raasch
 ! unused variables removed
 ! 
@@ -408 +411 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_header
+        ONLY:  bio_header
 
     USE bulk_cloud_model_mod,                                                  &
@@ -1151 +1154 @@
 !
 !-- Header information from other modules
-    IF ( biometeorology      )  CALL biom_header ( io )
+    IF ( biometeorology      )  CALL bio_header ( io )
     IF ( gust_module_enabled )  CALL gust_header( io )
     IF ( land_surface        )  CALL lsm_header( io )

palm/trunk/SOURCE/init_3d_model.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3524 2018-11-14 13:36:44Z raasch
 ! preprocessor directive added to avoid the compiler to complain about unused
 ! variable
@@ -554 +557 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_init, biom_init_arrays
+        ONLY:  bio_init, bio_init_arrays
 
     USE bulk_cloud_model_mod,                                                  &
@@ -1084 +1087 @@
 !
 !-- Allocate arrays for other modules
-    IF ( biometeorology      )  CALL biom_init_arrays
+    IF ( biometeorology      )  CALL bio_init_arrays
     IF ( bulk_cloud_model    )  CALL bcm_init_arrays
     IF ( gust_module_enabled )  CALL gust_init_arrays
@@ -2464 +2467 @@
     IF ( biometeorology )  THEN
         CALL location_message( 'initializing biometeorology module', .FALSE. )
-        CALL biom_init
+        CALL bio_init
         CALL location_message( 'finished', .TRUE. )
     ENDIF

palm/trunk/SOURCE/multi_agent_system_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3448 2018-10-29 18:14:31Z kanani
 ! Adjustment of biometeorology calls,
 ! implement some agent biometeorology
@@ -342 +345 @@
 
     USE biometeorology_mod,                                                   &
-        ONLY:  biom_calc_ipt, biom_determine_input_at
+        ONLY:  bio_calc_ipt, bio_calculate_mrt_grid, bio_get_thermal_index_input_ij
 
 
@@ -522 +525 @@
 !--    those agents to be deleted after the timestep
        deleted_agents = 0
-!
+
        IF ( biometeorology )  THEN
+!
+!--       Fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
+          CALL bio_calculate_mrt_grid ( .FALSE. )
 !
 !--       Call of human thermal comfort mod (and UV exposure)
@@ -537 +543 @@
 !
 !--             Evaluation of social forces
-!                CALL biom_dynamic(i,j)
+!                CALL bio_dynamic( i, j )
 !
 !--             Determine local meteorological conditions
-                CALL biom_determine_input_at ( .FALSE., i, j, ta, vp, v,      &
-                                               pair, tmrt )
+                CALL bio_get_thermal_index_input_ij ( .FALSE., i, j, ta, vp,  &
+                                                      v, pair, tmrt )
 
                 DO  a = 1, number_of_agents
@@ -547 +553 @@
 !--                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 )
+                   CALL bio_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
 

palm/trunk/SOURCE/netcdf_interface_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3485 2018-11-03 17:09:40Z gronemeier
 ! Write geographic coordinates as global attributes to file.
 ! 
@@ -607 +610 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_define_netcdf_grid
+        ONLY:  bio_define_netcdf_grid
 
     USE chemistry_model_mod,                                                   &
@@ -2032 +2035 @@
 !--                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 )
+                      CALL bio_define_netcdf_grid( do3d(av,i), found,          &
+                                                   grid_x, grid_y, grid_z )
                    ENDIF
 
@@ -3090 +3093 @@
 !--                      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 )
+                            CALL bio_define_netcdf_grid( do2d( av, i), found,  &
+                                                         grid_x, grid_y,       &
+                                                         grid_z )
                          ENDIF
 !

palm/trunk/SOURCE/parin.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3473 2018-10-30 20:50:15Z suehring
 ! Add virtual measurement module
 ! 
@@ -459 +462 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_parin
+        ONLY:  bio_parin
 
     USE bulk_cloud_model_mod,                                                  &
@@ -912 +915 @@
 !
 !--       Check for module namelists and read them
-          CALL biom_parin
+          CALL bio_parin
           CALL lsm_parin
           CALL bcm_parin

palm/trunk/SOURCE/sum_up_3d_data.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3467 2018-10-30 19:05:21Z suehring
 ! Implementation of a new aerosol module salsa.
 ! 
@@ -243 +246 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_3d_data_averaging
+        ONLY:  bio_3d_data_averaging
 
     USE bulk_cloud_model_mod,                                                  &
@@ -530 +533 @@
                 ENDIF
 
-                IF ( biometeorology  .AND.  trimvar(1:5) == 'biom_')  THEN
-                   CALL biom_3d_data_averaging( 'allocate', doav(ii) )
+                IF ( biometeorology  .AND.  trimvar(1:4) == 'bio_')  THEN
+                   CALL bio_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
@@ -1184 +1187 @@
              ENDIF
 
-             IF ( biometeorology  .AND.  trimvar(1:5) == 'biom_' )  THEN
-                   CALL biom_3d_data_averaging( 'sum', doav(ii) )
+             IF ( biometeorology  .AND.  trimvar(1:4) == 'bio_' )  THEN
+                CALL bio_3d_data_averaging( 'sum', doav(ii) )
              ENDIF
 

palm/trunk/SOURCE/time_integration.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Changes related to clean-up of biometeorology (dom_dwd_user)
+! 
+! 3524 2018-11-14 13:36:44Z raasch
 ! unused variables removed
 ! 
@@ -412 +415 @@
 
     USE biometeorology_mod,                                                    &
-        ONLY:  biom_calculate_static_grid, time_biom_results
+        ONLY:  bio_calculate_thermal_index_maps, time_bio_results
 
     USE bulk_cloud_model_mod,                                                  &
@@ -1363 +1366 @@
 !--    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
+          CALL bio_calculate_thermal_index_maps ( .FALSE. )
+          time_bio_results = time_since_reference_point
        ENDIF