Index: palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3525)
+++ palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3569)
@@ -21,5 +21,5 @@
 !
 ! Current revisions:
-! -----------------
+! ------------------
 ! 
 ! 
@@ -27,4 +27,13 @@
 ! -----------------
 ! $Id$
+! Consistently use bio_fill_value everywhere,
+! move allocation and initialization of output variables to bio_check_data_output
+! and bio_3d_data_averaging,
+! dom_dwd_user, Schrempf:
+! - integration of UVEM module part from r3474 (edited)
+! - split UTCI regression into 6 parts
+! - all data_output_3d is now explicity casted to sp
+!
+! 3525 2018-11-14 16:06:14Z kanani
 ! Clean up, renaming from "biom" to "bio", summary of thermal index calculation
 ! into one module (dom_dwd_user)
@@ -49,17 +58,21 @@
 ! Authors:
 ! --------
-! @author Dominik Froehlich <dominik.froehlich@dwd.de>
-! @author Jaroslav Resler <resler@cs.cas.cz>
+! @author Dominik Froehlich <dominik.froehlich@dwd.de>, thermal indices
+! @author Jaroslav Resler <resler@cs.cas.cz>, mean radiant temperature
+! @author Michael Schrempf <schrempf@muk.uni-hannover.de>, uv exposure
 !
 !
 ! Description:
 ! ------------
-!> Human thermal comfort module calculating thermal perception of a sample
+!> Biometeorology module consisting of two parts:
+!> 1.: Human thermal comfort module calculating thermal perception of a sample
 !> human being under the current meteorological conditions.
+!> 2.: Calculation of vitamin-D weighted UV exposure
 !>
 !> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
 !>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
 !> @todo Comments start with capital letter --> "!-- Include..." 
-!> @todo Variable and routine names strictly in lowercase letters and in English
+!> @todo uv_vitd3dose-->new output type necessary (cumulative)
+!> @todo consider upwelling radiation in UV
 !>
 !> @note nothing now
@@ -76,5 +89,5 @@
 
     USE basic_constants_and_equations_mod,                                     &
-       ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb
+       ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
 
     USE control_parameters,                                                    &
@@ -83,4 +96,7 @@
               simulated_time, surface_pressure
 
+    USE date_and_time_mod,                                                                                                        &
+        ONLY:  calc_date_and_time, day_of_year, time_utc
+
     USE grid_variables,                                                        &
        ONLY:  ddx, dx, ddy, dy
@@ -92,4 +108,8 @@
     USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
 
+    USE netcdf_data_input_mod,                                                                                                    &
+        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,                                                          &
+               uvem_irradiance_f, uvem_integration_f, building_obstruction_f
+!
 !-- Import radiation model to obtain input for mean radiant temperature
     USE radiation_model_mod,                                                   &
@@ -106,14 +126,16 @@
     PRIVATE
 
+!
 !-- Declare all global variables within the module (alphabetical order)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt 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)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci       !< UTCI results (Â°C)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet        !< PET results  (Â°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)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (Â°C)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av        !< PET results (aver. input)  (Â°C)
 
 
@@ -123,6 +145,7 @@
     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)
+    REAL ( wp ), PARAMETER ::  bio_fill_value = -9999._wp  !< set module fill value, replace by global fill value as soon as available
+!
 !--
-
     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?
@@ -135,4 +158,5 @@
     LOGICAL ::  average_trigger_pet   = .FALSE.  !< update averaged input on call to bio_pet?
 
+    LOGICAL ::  thermal_comfort = .TRUE.  !< Turn all thermal indices on or off
     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
@@ -142,4 +166,52 @@
     LOGICAL ::  bio_utci_av   = .TRUE.   !< Turn index UTCI (averaged input) on or off
 
+!
+!-- UVEM parameters from here
+!
+!-- Declare all global variables within the module (alphabetical order)
+    INTEGER(iwp) ::  ai                      = 0  !< loop index in azimuth direction
+    INTEGER(iwp) ::  bi                      = 0  !< loop index of bit location within an 8bit-integer (one Byte)
+    INTEGER(iwp) ::  clothing                = 1  !< clothing (0=unclothed, 1=Arms,Hands,Face free, 3=Hand,Face free)
+    INTEGER(iwp) ::  iq                      = 0  !< loop index of irradiance quantity
+    INTEGER(iwp) ::  pobi                    = 0  !< loop index of the position of corresponding byte within ibset byte vektor
+    INTEGER(iwp) ::  obstruction_direct_beam = 0  !< Obstruction information for direct beam    
+    INTEGER(iwp) ::  zi                      = 0  !< loop index in zenith direction
+
+    INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0  !< temporary obstruction information stored with ibset 
+    INTEGER(iwp),    DIMENSION(0:359) ::  obstruction_temp2 = 0  !< restored temporary obstruction information from ibset file
+
+    INTEGER(iwp), DIMENSION(0:35,0:9) ::  obstruction       = 1  !< final 2D obstruction information array
+
+    LOGICAL ::  consider_obstructions = .TRUE.   !< namelist parameter (see documentation)
+    LOGICAL ::  sun_in_south          = .FALSE.  !< namelist parameter (see documentation)
+    LOGICAL ::  turn_to_sun           = .TRUE.   !< namelist parameter (see documentation)
+    LOGICAL ::  uv_exposure           = .FALSE.  !< namelist parameter (see documentation)
+
+    REAL(wp) ::  diffuse_exposure            =   0.0_wp  !< calculated exposure by diffuse radiation
+    REAL(wp) ::  direct_exposure             =   0.0_wp  !< calculated exposure by direct solar beam   
+    REAL(wp) ::  orientation_angle           =   0.0_wp  !< orientation of front/face of the human model    
+    REAL(wp) ::  projection_area_direct_beam =   0.0_wp  !< projection area for direct solar beam
+    REAL(wp) ::  saa                         = 180.0_wp  !< solar azimuth angle
+    REAL(wp) ::  startpos_human              =   0.0_wp  !< start value for azimuth interpolation of human geometry array
+    REAL(wp) ::  startpos_saa_float          =   0.0_wp  !< start value for azimuth interpolation of radiance array
+    REAL(wp) ::  sza                         =  20.0_wp  !< solar zenith angle
+    REAL(wp) ::  xfactor                     =   0.0_wp  !< relative x-position used for interpolation
+    REAL(wp) ::  yfactor                     =   0.0_wp  !< relative y-position used for interpolation
+
+    REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp  !< iradiance values extracted from irradiance lookup table  
+
+    REAL(wp), DIMENSION(0:2,0:90) ::  irradiance_lookup_table      = 0.0_wp  !< irradiance lookup table
+    REAL(wp), DIMENSION(0:35,0:9) ::  integration_array            = 0.0_wp  !< solid angle factors for hemispherical integration
+    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area              = 0.0_wp  !< projection areas of a human (all directions)
+    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area_lookup_table = 0.0_wp  !< human geometry lookup table (projection areas)
+    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_direct_temp  = 0.0_wp  !< temporary projection area for direct solar beam
+    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_temp         = 0.0_wp  !< temporary projection area for all directions
+    REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp  !< radiance extracted from radiance_lookup_table  
+    REAL(wp), DIMENSION(0:71,0:9) ::  radiance_array_temp          = 0.0_wp  !< temporary radiance data
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure     !< result variable for instantaneous vitamin-D weighted exposures
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure_av  !< result variable for summation of vitamin-D weighted exposures
+
+    REAL(wp), DIMENSION(0:35,0:9,0:90) ::  radiance_lookup_table   = 0.0_wp  !< radiance lookup table
 
 !
@@ -150,6 +222,9 @@
     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
+    bio_init, bio_parin, bio_perct, bio_perct_av, bio_pet,                     &
+    bio_pet_av, bio_utci, bio_utci_av, thermal_comfort, time_bio_results,      &
+!
+!-- UVEM PUBLIC variables and methods
+    uvem_calc_exposure, uv_exposure
 
 !
@@ -160,65 +235,60 @@
        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 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 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 bio_check_data_output
-        MODULE PROCEDURE bio_check_data_output
-     END INTERFACE bio_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 bio_check_parameters
        MODULE PROCEDURE bio_check_parameters
     END INTERFACE bio_check_parameters
-
+!
 !-- Data output of 2D quantities
     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 bio_data_output_3d
        MODULE PROCEDURE bio_data_output_3d
     END INTERFACE bio_data_output_3d
-
+!
 !-- Definition of data output quantities
     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 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 bio_header
        MODULE PROCEDURE bio_header
     END INTERFACE bio_header
-
+!
 !-- Initialization actions 
     INTERFACE bio_init
        MODULE PROCEDURE bio_init
     END INTERFACE bio_init
-
-!-- Initialization of arrays
-    INTERFACE bio_init_arrays
-       MODULE PROCEDURE bio_init_arrays
-    END INTERFACE bio_init_arrays
-
+!
 !-- Reading of NAMELIST parameters
     INTERFACE bio_parin
@@ -226,4 +296,9 @@
     END INTERFACE bio_parin
 
+!
+!-- Calculate UV exposure grid
+    INTERFACE uvem_calc_exposure
+       MODULE PROCEDURE uvem_calc_exposure
+    END INTERFACE uvem_calc_exposure
 
  CONTAINS
@@ -235,4 +310,8 @@
 !> Sum up and time-average biom input quantities as well as allocate
 !> the array necessary for storing the average.
+!> There is a considerable difference to the 3d_data_averaging subroutines
+!> used by other modules:
+!> For the thermal indices, the module needs to average the input conditions
+!> not the result!
 !------------------------------------------------------------------------------!
  SUBROUTINE bio_3d_data_averaging( mode, variable )
@@ -260,8 +339,36 @@
           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_perct .AND. .NOT. average_trigger_utci &
-                .AND. .NOT. average_trigger_pet ) THEN
+!
+!--          Averaging, as well as the allocation of the required grids must be
+!            done only once, independent from for how many thermal indices 
+!            averaged output is desired.
+!            Therefore wee need to memorize which index is the one that controls
+!            the averaging (what must be the first thermal index called).
+!            Indices are in unknown order as depending on the input file, 
+!            determine first index to average und update only once
+
+!--          Only proceed here if this was not done for any index before. This
+!            is done only once during the whole model run.
+             IF ( .NOT. average_trigger_perct .AND.                            &
+                  .NOT. average_trigger_utci  .AND.                            &
+                  .NOT. average_trigger_pet ) THEN
+!
+!--             Allocate the required grids
+                IF ( .NOT. ALLOCATED( perct_av ) ) THEN
+                   ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
+                ENDIF
+                perct_av = REAL( bio_fill_value, KIND = wp )
+
+                IF ( .NOT. ALLOCATED( utci_av ) ) THEN
+                   ALLOCATE( utci_av (nys:nyn,nxl:nxr) )
+                ENDIF
+                utci_av = REAL( bio_fill_value, KIND = wp )
+
+                IF ( .NOT. ALLOCATED( pet_av ) ) THEN
+                   ALLOCATE( pet_av (nys:nyn,nxl:nxr) )
+                ENDIF
+                pet_av = REAL( bio_fill_value, KIND = wp )
+!
+!--             Memorize the first index called to control averaging
                 IF ( TRIM( variable ) == 'bio_perct*' ) THEN
                     average_trigger_perct = .TRUE.
@@ -274,40 +381,48 @@
                 ENDIF
              ENDIF
-
-!--          Only continue if updateindex
+!
+!--          Only continue if var is the index, that controls averaging.
+!            Break immediatelly (doing nothing) for the other indices.
              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!
+!
+!--          Now memorize which of the input grids are not averaged by other 
+!            modules. Set averaging switch to .TRUE. in that case.
              IF ( .NOT. ALLOCATED( pt_av ) )  THEN
                 ALLOCATE( pt_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
                 aver_perct = .TRUE.
+                pt_av = 0.0_wp
              ENDIF
-             pt_av = 0.0_wp
 
              IF ( .NOT. ALLOCATED( q_av ) )  THEN
                 ALLOCATE( q_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
                 aver_q = .TRUE.
+                q_av = 0.0_wp
              ENDIF
-             q_av = 0.0_wp
 
              IF ( .NOT. ALLOCATED( u_av ) )  THEN
                 ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 aver_u = .TRUE.
+                u_av = 0.0_wp
              ENDIF
-             u_av = 0.0_wp
 
              IF ( .NOT. ALLOCATED( v_av ) )  THEN
                 ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 aver_v = .TRUE.
+                v_av = 0.0_wp
              ENDIF
-             v_av = 0.0_wp
 
              IF ( .NOT. ALLOCATED( w_av ) )  THEN
                 ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
                 aver_w = .TRUE.
+                w_av = 0.0_wp
              ENDIF
-             w_av = 0.0_wp
+
+          CASE ( 'uvem_vitd3dose*' )
+             IF ( .NOT. ALLOCATED( vitd3_exposure_av ) )  THEN
+                ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             vitd3_exposure_av = 0.0_wp
 
           CASE DEFAULT
@@ -335,6 +450,6 @@
 
           CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
-
-!--          Only continue if updateindex
+!
+!--          Only continue if updateindex, see above
              IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
                 RETURN
@@ -393,6 +508,16 @@
                 ENDDO
              ENDIF
-
-           CASE DEFAULT
+!
+!--       This is a cumulated dose. No mode == 'average' for this quantity.
+          CASE ( 'uvem_vitd3dose*' )
+             IF ( ALLOCATED( vitd3_exposure_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      vitd3_exposure_av(j,i) = vitd3_exposure_av(j,i) + vitd3_exposure(j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE DEFAULT
              CONTINUE
 
@@ -409,12 +534,12 @@
 
           CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
-
-!--          Only continue if update index
-             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
+!
+!--          Only continue if update index, see above
+             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
 
              IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
@@ -422,5 +547,6 @@
                    DO  j = nys, nyn
                       DO  k = nzb, nzt+1
-                         pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                         pt_av(k,j,i) = pt_av(k,j,i) /                         &
+                            REAL( average_count_3d, KIND=wp )
                       ENDDO
                    ENDDO
@@ -432,5 +558,6 @@
                    DO  j = nys, nyn
                       DO  k = nzb, nzt+1
-                         q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                         q_av(k,j,i) = q_av(k,j,i) /                           &
+                            REAL( average_count_3d, KIND=wp )
                       ENDDO
                    ENDDO
@@ -442,5 +569,6 @@
                    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 )
+                         u_av(k,j,i) = u_av(k,j,i) /                           &
+                            REAL( average_count_3d, KIND=wp )
                       ENDDO
                    ENDDO
@@ -452,5 +580,6 @@
                    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 )
+                         v_av(k,j,i) = v_av(k,j,i) /                           &
+                            REAL( average_count_3d, KIND=wp )
                       ENDDO
                    ENDDO
@@ -462,12 +591,17 @@
                    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 )
+                         w_av(k,j,i) = w_av(k,j,i) /                           &
+                            REAL( average_count_3d, KIND=wp )
                       ENDDO
                    ENDDO
                 ENDDO
              ENDIF
-
-!--          Udate thermal indices with derived averages
+!
+!--          Udate all thermal index grids with updated averaged input
              CALL bio_calculate_thermal_index_maps ( .TRUE. )
+
+!
+!--     No averaging for UVEM since we are calculating a dose (only sum is
+!--     calculated and saved to av.nc file)
 
         END SELECT
@@ -485,5 +619,5 @@
 !> Check data output for biometeorology model
 !------------------------------------------------------------------------------!
- SUBROUTINE bio_check_data_output( var, unit )
+ SUBROUTINE bio_check_data_output( var, unit, i, ilen, k )
 
     USE control_parameters,                                                    &
@@ -495,9 +629,74 @@
     CHARACTER (LEN=*) ::  var      !< The variable in question
 
+    INTEGER(iwp) ::  i      !<
+    INTEGER(iwp) ::  ilen   !<   
+    INTEGER(iwp) ::  k      !<
 
     SELECT CASE ( TRIM( var ) )
-     
-       CASE( 'bio_mrt', 'bio_pet*', 'bio_perct*', 'bio_utci*' )
+!
+!-- Allocate a temporary array with the desired output dimensions.
+       CASE ( 'bio_mrt')
           unit = 'degree_C'
+          IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
+             ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
+          ENDIF
+          tmrt_grid = REAL( bio_fill_value, KIND = wp )
+
+       CASE ( 'bio_perct*' )
+          unit = 'degree_C'
+          IF ( .NOT. ALLOCATED( perct ) )  THEN
+             ALLOCATE( perct (nys:nyn,nxl:nxr) )
+          ENDIF
+          perct = REAL( bio_fill_value, KIND = wp )
+
+       CASE ( 'bio_utci*' )
+          unit = 'degree_C'
+          IF ( .NOT. ALLOCATED( utci ) )  THEN
+             ALLOCATE( utci (nys:nyn,nxl:nxr) )
+          ENDIF
+          utci = REAL( bio_fill_value, KIND = wp )
+
+       CASE ( 'bio_pet*' )
+          unit = 'degree_C'
+          IF ( .NOT. ALLOCATED( pet ) )  THEN
+             ALLOCATE( pet (nys:nyn,nxl:nxr) )
+          ENDIF
+          pet = REAL( bio_fill_value, KIND = wp )
+
+       CASE ( 'uvem_vitd3*' )
+          IF (  .NOT.  uv_exposure )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
+                      'res a namelist &uvexposure_par'
+             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+             message_string = 'illegal value for data_output: "' //            &
+                              TRIM( var ) // '" & only 2d-horizontal ' //      &
+                              'cross sections are allowed for this value'
+             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+          ENDIF
+          unit = 'IU/s'
+          IF ( .NOT. ALLOCATED( vitd3_exposure ) )  THEN
+             ALLOCATE( vitd3_exposure(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          vitd3_exposure = 0.0_wp
+
+       CASE ( 'uvem_vitd3dose*' )
+          IF (  .NOT.  uv_exposure )  THEN
+             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
+                      'res  a namelist &uvexposure_par'
+             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
+             message_string = 'illegal value for data_output: "' //            &
+                              TRIM( var ) // '" & only 2d-horizontal ' //      &
+                              'cross sections are allowed for this value'
+             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
+          ENDIF
+          unit = 'IU/av-h'
+          IF ( .NOT. ALLOCATED( vitd3_exposure_av ) )  THEN
+             ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
+          ENDIF
+          vitd3_exposure_av = 0.0_wp
 
        CASE DEFAULT
@@ -506,7 +705,8 @@
     END SELECT
 
-    IF ( unit /= 'illegal' )  THEN
+    IF ( thermal_comfort .AND. unit == 'degree_C' )  THEN
 !
 !--    Futher checks if output belongs to biometeorology
+!      Break if required modules "radiation" and "humidity" are not running.
        IF ( .NOT.  radiation ) THEN
           message_string = 'output of "' // TRIM( var ) // '" require'         &
@@ -541,8 +741,8 @@
     IMPLICIT NONE
 
-
+!
 !--    Disabled as long as radiation model not available
 
-       IF ( .NOT. humidity )  THEN
+       IF ( thermal_comfort  .AND.  .NOT.  humidity )  THEN
           message_string = 'The estimation of thermal comfort requires '    // &
                            'air humidity information, but humidity module ' // &
@@ -562,8 +762,6 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf,           &
-                                      two_d, nzb_do, nzt_do, fill_value )
-
-    USE indices,                                                               &
-       ONLY: nxl, nxl, nxr, nxr, nyn, nyn, nys, nys, nzb, nzt
+                                two_d, nzb_do, nzt_do )
+
 
     USE kinds
@@ -571,5 +769,5 @@
 
     IMPLICIT NONE
-
+!
 !-- Input variables
     CHARACTER (LEN=*), INTENT(IN) ::  variable    !< Char identifier to select var for output
@@ -577,6 +775,5 @@
     INTEGER(iwp), INTENT(IN)      ::  nzb_do      !< Unused. 2D. nz bottom to nz top
     INTEGER(iwp), INTENT(IN)      ::  nzt_do      !< Unused.
-    REAL(wp), INTENT(in)          ::  fill_value  !< Fill value for unassigned locations
-
+!
 !-- Output variables
     CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
@@ -584,7 +781,6 @@
     LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables, horizontal cross sections, must be .TRUE.
     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return
-
+!
 !-- Internal variables
-    CHARACTER (LEN=:), allocatable ::  variable_short  !< Trimmed version of char variable
     INTEGER(iwp) ::  i        !< Running index, x-dir
     INTEGER(iwp) ::  j        !< Running index, y-dir
@@ -593,20 +789,166 @@
 
 
-    variable_short = TRIM( variable )
-    IF ( variable_short(1:4) /= 'bio_' ) THEN
-       found = .FALSE.
-       grid  = 'none'
-    ENDIF
-
     found = .TRUE.
-    local_pf = fill_value
-
-    SELECT CASE ( variable_short )
+    local_pf = bio_fill_value
+
+    SELECT CASE ( TRIM( variable ) )
 
 
         CASE ( 'bio_mrt_xy' )
-            grid = 'zu1'
-            two_d = .FALSE.  !< can be calculated for several levels
-            local_pf = REAL( fill_value, KIND = wp )
+           grid = 'zu1'
+           two_d = .FALSE.  !< can be calculated for several levels
+           local_pf = REAL( bio_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) = perct(j,i)
+                 ENDDO
+              ENDDO
+           ELSE
+              DO  i = nxl, nxr
+                 DO  j = nys, nyn
+                    local_pf(i,j,nzb+1) = perct_av(j,i)
+                 ENDDO
+              ENDDO
+           END IF
+
+
+        CASE ( 'bio_utci*_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) = utci(j,i)
+                 ENDDO
+              ENDDO
+           ELSE
+              DO  i = nxl, nxr
+                 DO  j = nys, nyn
+                    local_pf(i,j,nzb+1) = utci_av(j,i)
+                 ENDDO
+              ENDDO
+           END IF
+
+
+        CASE ( 'bio_pet*_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) = pet(j,i)
+                 ENDDO
+              ENDDO
+           ELSE
+              DO  i = nxl, nxr
+                 DO  j = nys, nyn
+                    local_pf(i,j,nzb+1) = pet_av(j,i)
+                 ENDDO
+              ENDDO
+           END IF
+
+           !
+!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein. 
+!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged. 
+       CASE ( 'uvem_vitd3*_xy' )        ! 2d-array
+          IF ( av == 0 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   local_pf(i,j,nzb+1) = vitd3_exposure(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+       CASE ( 'uvem_vitd3dose*_xy' )        ! 2d-array
+          IF ( av == 1 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   local_pf(i,j,nzb+1) = vitd3_exposure_av(j,i)
+                ENDDO
+             ENDDO
+          ENDIF
+
+          two_d = .TRUE.
+          grid = 'zu1'
+
+
+       CASE DEFAULT
+          found = .FALSE.
+          grid  = 'none'
+
+    END SELECT
+
+
+ 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
+
+!     REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
+
+    found = .TRUE.
+
+    SELECT CASE ( TRIM( variable ) )
+
+        CASE ( 'bio_mrt' )
+            local_pf = REAL( bio_fill_value, KIND = sp )
             DO  l = 1, nmrtbl
                i = mrtbl(ix,l)
@@ -617,147 +959,14 @@
                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
+                     local_pf(i,j,k) = REAL(((human_absorb * mrtinsw(l) +      &
+                                    human_emiss * mrtinlw(l)) /                &
+                                    (human_emiss * sigma_sb)) ** .25_wp -      &
+                                    degc_to_k, kind=sp )
                   ELSE
-                     local_pf(i,j,k) = (human_emiss * mrtinlw(l) /             &
-                                        sigma_sb) ** .25_wp - degc_to_k
+                     local_pf(i,j,k) = REAL((human_emiss * mrtinlw(l) /        &
+                                    sigma_sb) ** .25_wp - degc_to_k, kind=sp )  !< why not (human_emiss * sigma_sb) as above?
                   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) = perct(j,i)
-                 ENDDO
-              ENDDO
-            ELSE
-              DO  i = nxl, nxr
-                 DO  j = nys, nyn
-                    local_pf(i,j,nzb+1) = perct_av(j,i)
-                 ENDDO
-              ENDDO
-            END IF
-
-
-        CASE ( 'bio_utci*_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) = utci_grid(j,i)
-                 ENDDO
-              ENDDO
-            ELSE
-              DO  i = nxl, nxr
-                 DO  j = nys, nyn
-                    local_pf(i,j,nzb+1) = utci_av_grid(j,i)
-                 ENDDO
-              ENDDO
-            END IF
-
-
-        CASE ( 'bio_pet*_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) = pet_grid(j,i)
-                 ENDDO
-              ENDDO
-            ELSE
-              DO  i = nxl, nxr
-                 DO  j = nys, nyn
-                    local_pf(i,j,nzb+1) = pet_av_grid(j,i)
-                 ENDDO
-              ENDDO
-            END IF
-
-
-       CASE DEFAULT
-          found = .FALSE.
-          grid  = 'none'
-
-    END SELECT
-
-
- END SUBROUTINE 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)
+                  local_pf(i,j,k) = REAL(mrt_av_grid(l), kind=sp)
                ENDIF
             ENDDO
@@ -780,8 +989,8 @@
 
     IMPLICIT NONE
-
+!
 !-- Input variables
     CHARACTER (LEN=*), INTENT(IN)  ::  var      !< Name of output variable
-
+!
 !-- Output variables
     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
@@ -790,5 +999,5 @@
 
     LOGICAL, INTENT(OUT)           ::  found    !< Flag if output var is found
-
+!
 !-- Local variables
     LOGICAL      :: is2d  !< Var is 2d?
@@ -810,5 +1019,11 @@
        grid_y = 'y'
        grid_z = 'zu'
-       IF ( is2d ) grid_z = 'zu1'
+       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' ) grid_z = 'zu1'
+    ENDIF
+
+    IF ( is2d  .AND.  var(1:4) == 'uvem' ) THEN
+       grid_x = 'x'
+       grid_y = 'y'
+       grid_z = 'zu1'
     ENDIF
 
@@ -824,8 +1039,8 @@
 
     IMPLICIT NONE
-
+!
 !-- Input variables
     INTEGER(iwp), INTENT(IN) ::  io           !< Unit of the output file
-
+!
 !-- Internal variables
     CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
@@ -857,11 +1072,13 @@
         ONLY: message_string
 
+    USE netcdf_data_input_mod,                                                &
+        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,      &
+               uvem_irradiance_f, uvem_integration_f, building_obstruction_f
+
     IMPLICIT NONE
-
+!
 !-- Internal vriables
     REAL ( wp )  :: height  !< current height in meters
-
-    INTEGER ( iwp )  :: i  !< iteration index
-
+!
 !-- Determine cell level corresponding to 1.1 m above ground level
 !   (gravimetric center of sample human)
@@ -876,67 +1093,15 @@
 
     IF ( .NOT. radiation_interactions ) THEN
-       message_string = 'The mrt calculation requires ' // &
-                        'enabled radiation_interactions but it ' // &
+       message_string = 'The mrt calculation requires ' //                     &
+                        'enabled radiation_interactions but it ' //            &
                         'is disabled!'
        CALL message( 'check_parameters', 'PAHU03', 0, 0, -1, 6, 0 )
     ENDIF
 
+!
+!-- Init UVEM and load lookup tables
+    CALL netcdf_data_input_uvem
+
  END SUBROUTINE bio_init
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Allocate biom arrays and define pointers if required
-!> init_3d_model 1047ff
-!------------------------------------------------------------------------------!
- SUBROUTINE bio_init_arrays
-
-    IMPLICIT NONE
-
-!-- Allocate a temporary array with the desired output dimensions.
-!   Initialization omitted for performance, will be overwritten anyway
-    IF ( .NOT. ALLOCATED( tmrt_grid ) ) THEN
-      ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
-    ENDIF
-
-    IF ( bio_perct ) THEN
-      IF ( .NOT. ALLOCATED( perct ) ) THEN
-        ALLOCATE( perct (nys:nyn,nxl:nxr) )
-      ENDIF
-    ENDIF
-
-    IF ( bio_utci ) THEN
-      IF ( .NOT. ALLOCATED( utci_grid ) ) THEN
-        ALLOCATE( utci_grid (nys:nyn,nxl:nxr) )
-      ENDIF
-    ENDIF
-
-    IF ( bio_pet ) THEN
-      IF ( .NOT. ALLOCATED( pet_grid ) ) THEN
-        ALLOCATE( pet_grid (nys:nyn,nxl:nxr) )
-      ENDIF
-    END IF
-
-    IF ( bio_perct_av ) THEN
-      IF ( .NOT. ALLOCATED( perct_av ) ) THEN
-        ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
-      ENDIF
-    ENDIF
-
-    IF ( bio_utci_av ) THEN
-      IF ( .NOT. ALLOCATED( utci_av_grid ) ) THEN
-        ALLOCATE( utci_av_grid (nys:nyn,nxl:nxr) )
-      ENDIF
-    ENDIF
-
-    IF ( bio_pet_av ) THEN
-      IF ( .NOT. ALLOCATED( pet_av_grid ) ) THEN
-        ALLOCATE( pet_av_grid (nys:nyn,nxl:nxr) )
-      ENDIF
-
-    END IF
-
- END SUBROUTINE bio_init_arrays
 
 
@@ -959,5 +1124,14 @@
                                           bio_perct_av,                        &
                                           bio_utci,                            &
-                                          bio_utci_av
+                                          bio_utci_av,                         &
+                                          thermal_comfort,                     &
+!
+!-- UVEM namelist parameters
+                                          clothing,                            &
+                                          consider_obstructions,               &
+                                          orientation_angle,                   &
+                                          sun_in_south,                        &
+                                          turn_to_sun,                         &
+                                          uv_exposure
 
 
@@ -1003,4 +1177,5 @@
 
     LOGICAL, INTENT(IN)         ::  av    !< use averaged input?
+!
 !-- Internal variables
     INTEGER(iwp)                ::  i     !< Running index, x-dir, radiation coordinates
@@ -1009,8 +1184,8 @@
     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.
+!   fluxes calculated by RTM and assign into 2D grid
+    tmrt_grid = bio_fill_value
     DO  l = 1, nmrtbl
        i = mrtbl(ix,l)
@@ -1042,10 +1217,10 @@
 
     IMPLICIT NONE
-
+!
 !-- Input variables
     LOGICAL,      INTENT ( IN ) ::  average_input  !< Determine averaged input conditions?
     INTEGER(iwp), INTENT ( IN ) ::  i     !< Running index, x-dir
     INTEGER(iwp), INTENT ( IN ) ::  j     !< Running index, y-dir
-
+!
 !-- Output parameters
     REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (Â°C)
@@ -1054,16 +1229,16 @@
     REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
     REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
-
+!
 !-- Internal variables
     INTEGER(iwp)                ::  k     !< Running index, z-dir
-    INTEGER(iwp)                ::  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)                ::  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
+!     INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
 
     REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
 
-
+!
 !-- Determine cell level closest to 1.1m above ground
 !   by making use of truncation due to int cast
@@ -1076,11 +1251,12 @@
        k_wind = k + 1_iwp
     ENDIF
-
+!
 !-- Determine local values:
     IF ( average_input ) THEN
+!
 !--    Calculate ta from Tp assuming dry adiabatic laps rate
        ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
 
-       vp = -999._wp
+       vp = bio_fill_value
        IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
           vp = q_av(k, j, i)
@@ -1091,8 +1267,9 @@
           0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) )
     ELSE
+!
 !-- Calculate ta from Tp assuming dry adiabatic laps rate
        ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
 
-       vp = -999._wp
+       vp = bio_fill_value
        IF ( humidity ) THEN
           vp = q(k, j, i)
@@ -1104,5 +1281,5 @@
 
     ENDIF
-
+!
 !-- Local air pressure
     pair = surface_pressure
@@ -1116,8 +1293,9 @@
        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)
+    tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
     IF ( radiation ) THEN
+!
 !--    Use MRT from RTM precalculated in tmrt_grid
        tmrt = tmrt_grid(j,i)
@@ -1136,58 +1314,67 @@
 
     IMPLICIT NONE
-
+!
 !-- Input attributes
     LOGICAL, INTENT ( IN ) ::  av  !< Calculate based on averaged input conditions?
-
+!
 !-- Internal variables
-    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) ::  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
-
+    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) ::  perct_ij     !< Perceived temperature            (Â°C)
+    REAL(wp) ::  utci_ij      !< Universal thermal climate index  (Â°C)
+    REAL(wp) ::  pet_ij       !< Physiologically equivalent temperature  (Â°C)
+    REAL(wp) ::  tmrt_ij      !< Mean radiant temperature         (Â°C)
+
+!
 !-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
-    CALL bio_calculate_mrt_grid ( av )
-
+    IF ( simulated_time > 0.0_wp ) THEN
+       CALL bio_calculate_mrt_grid ( av )
+    ENDIF
 
     DO i = nxl, nxr
        DO j = nys, nyn
+!
 !--       Determine local input conditions
-          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!
-
+          tmrt_ij = bio_fill_value
+          vp      = bio_fill_value
+!
+!--       Determine input only if 
+          IF ( simulated_time > 0.0_wp ) THEN
+             CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,           &
+                ws, pair, tmrt_ij )
+          END IF
+!
 !--       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
-          ELSE
+          pet_ij   = bio_fill_value   !< set fail value, e.g. valid for within
+          perct_ij = bio_fill_value   !< some obstacle
+          utci_ij  = bio_fill_value
+          IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp ) ) THEN
+!
 !--          Calculate static thermal indices based on local tmrt
-             clo = -999._wp
-             
+             clo = bio_fill_value
+
              IF ( bio_perct ) THEN
+!
 !--          Estimate local perceived temperature
-                CALL calculate_perct_static( ta, vp, ws, tmrt_tmp, pair, clo,  &
-                   perct_tmp )
+                CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, clo,   &
+                   perct_ij )
              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 )
+                CALL calculate_utci_static( ta, vp, ws, tmrt_ij,               &
+                   bio_output_height, utci_ij )
              ENDIF
-             
+
              IF ( bio_pet ) THEN
+!
 !--          Estimate local physiologically equivalent temperature
-                CALL calculate_pet_static( ta, vp, ws, tmrt_tmp, pair, pet_tmp )
+                CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij )
              ENDIF
           END IF
@@ -1195,24 +1382,26 @@
 
           IF ( av ) THEN
+!
 !--          Write results for selected averaged indices
              IF ( bio_perct_av )  THEN
-                perct_av(j, i) = perct_tmp
+                perct_av(j, i) = perct_ij
              END IF
              IF ( bio_utci_av ) THEN
-                utci_av_grid(j, i) = utci_tmp
+                utci_av(j, i) = utci_ij
              END IF
              IF ( bio_pet_av ) THEN
-                pet_av_grid(j, i)  = pet_tmp
+                pet_av(j, i)  = pet_ij
              END IF
           ELSE
+!
 !--          Write result for selected indices
              IF ( bio_perct )  THEN
-                perct(j, i) = perct_tmp
+                perct(j, i) = perct_ij
              END IF
              IF ( bio_utci ) THEN
-                utci_grid(j, i) = utci_tmp
+                utci(j, i) = utci_ij
              END IF
              IF ( bio_pet ) THEN
-                pet_grid(j, i)  = pet_tmp
+                pet(j, i)  = pet_ij
              END IF
           END IF
@@ -1232,5 +1421,5 @@
 
     IMPLICIT NONE
-
+!
 !-- Input parameters
     REAL(wp), INTENT ( IN )  ::  ta   !< Air temperature                  (Â°C)
@@ -1246,5 +1435,5 @@
                                          !  (without metabolism!)         (W)
     INTEGER(iwp), INTENT ( IN ) ::  sex  !< Sex of agent (1 = male, 2 = female)
-
+!
 !-- Both, input and output parameters
     Real(wp), INTENT ( INOUT )  ::  energy_storage    !< Energy storage   (W/mÂ²)
@@ -1253,9 +1442,11 @@
     Real(wp), INTENT ( INOUT )  ::  actlev  !< Individuals activity level
                                             !  per unit surface area      (W/mÂ²)
+!
 !-- Output parameters
     REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (Â°C)
-
+!
 !-- If clo equals the initial value, this is the initial call
     IF ( clo <= -998._wp ) THEN
+!
 !--    Initialize instationary perceived temperature with personalized 
 !      PT as an initial guess, set actlev and clo
@@ -1264,4 +1455,5 @@
           ipt )
     ELSE
+!
 !--    Estimate local instatinoary perceived temperature
        CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,     &
@@ -1287,43 +1479,71 @@
 !> www.utci.org
 !------------------------------------------------------------------------------!
- SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci )
+ SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci_ij )
 
     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)
+    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)
-
+    REAL(wp), INTENT ( OUT ) ::  utci_ij  !< Universal Thermal Climate Index (Â°C)
+
+    REAL(WP) ::  ta           !< air temperature modified by offset (Â°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)
+    REAL(WP) ::  part_ta      !< Air temperature related part of the regression
+    REAL(WP) ::  ta2          !< 2 times ta
+    REAL(WP) ::  ta3          !< 3 times ta
+    REAL(WP) ::  ta4          !< 4 times ta
+    REAL(WP) ::  ta5          !< 5 times ta
+    REAL(WP) ::  ta6          !< 6 times ta
+    REAL(WP) ::  part_va      !< Vapour pressure related part of the regression
+    REAL(WP) ::  va2          !< 2 times va
+    REAL(WP) ::  va3          !< 3 times va
+    REAL(WP) ::  va4          !< 4 times va
+    REAL(WP) ::  va5          !< 5 times va
+    REAL(WP) ::  va6          !< 6 times va
+    REAL(WP) ::  part_d_tmrt  !< Mean radiant temp. related part of the reg.
+    REAL(WP) ::  d_tmrt2      !< 2 times d_tmrt
+    REAL(WP) ::  d_tmrt3      !< 3 times d_tmrt
+    REAL(WP) ::  d_tmrt4      !< 4 times d_tmrt
+    REAL(WP) ::  d_tmrt5      !< 5 times d_tmrt
+    REAL(WP) ::  d_tmrt6      !< 6 times d_tmrt
+    REAL(WP) ::  part_pa      !< Air pressure related part of the regression
+    REAL(WP) ::  pa2          !< 2 times pa
+    REAL(WP) ::  pa3          !< 3 times pa
+    REAL(WP) ::  pa4          !< 4 times pa
+    REAL(WP) ::  pa5          !< 5 times pa
+    REAL(WP) ::  pa6          !< 6 times pa
+    REAL(WP) ::  part_pa2     !< Air pressure^2 related part of the regression
+    REAL(WP) ::  part_pa3     !< Air pressure^3 related part of the regression
+    REAL(WP) ::  part_pa46    !< Air pressure^4-6 related part of the regression
+
+!
 !-- 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
+       utci_ij = bio_fill_value
        RETURN
     ENDIF
-
+!
 !-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
     IF ( ta > 50._wp ) THEN
@@ -1335,5 +1555,5 @@
        ta = -50._wp
     ENDIF
-
+!
 !-- For routine application. For wind speeds and relative
 !   humidity values below 0.5 m/s or 5%, respectively, the
@@ -1342,220 +1562,256 @@
     IF ( va > 17._wp ) va = 17._wp
 
+!
+!-- Pre-calculate multiples of input parameters to save time later
+
+    ta2 = ta  * ta
+    ta3 = ta2 * ta
+    ta4 = ta3 * ta
+    ta5 = ta4 * ta
+    ta6 = ta5 * ta
+
+    va2 = va  * va
+    va3 = va2 * va
+    va4 = va3 * va
+    va5 = va4 * va
+    va6 = va5 * va
+
+    d_tmrt2 = d_tmrt  * d_tmrt
+    d_tmrt3 = d_tmrt2 * d_tmrt
+    d_tmrt4 = d_tmrt3 * d_tmrt
+    d_tmrt5 = d_tmrt4 * d_tmrt
+    d_tmrt6 = d_tmrt5 * d_tmrt
+
+    pa2 = pa  * pa
+    pa3 = pa2 * pa
+    pa4 = pa3 * pa
+    pa5 = pa4 * pa
+    pa6 = pa5 * pa
+
+!
+!-- Pre-calculate parts of the regression equation
+    part_ta = (  6.07562052e-01_wp )       +                                   &
+              ( -2.27712343e-02_wp ) * ta  +                                   &
+              (  8.06470249e-04_wp ) * ta2 +                                   &
+              ( -1.54271372e-04_wp ) * ta3 +                                   &
+              ( -3.24651735e-06_wp ) * ta4 +                                   &
+              (  7.32602852e-08_wp ) * ta5 +                                   &
+              (  1.35959073e-09_wp ) * ta6
+
+    part_va = ( -2.25836520e+00_wp ) * va +                                    &
+        (  8.80326035e-02_wp ) * ta  * va +                                    &
+        (  2.16844454e-03_wp ) * ta2 * va +                                    &
+        ( -1.53347087e-05_wp ) * ta3 * va +                                    &
+        ( -5.72983704e-07_wp ) * ta4 * va +                                    &
+        ( -2.55090145e-09_wp ) * ta5 * va +                                    &
+        ( -7.51269505e-01_wp ) *       va2 +                                   &
+        ( -4.08350271e-03_wp ) * ta  * va2 +                                   &
+        ( -5.21670675e-05_wp ) * ta2 * va2 +                                   &
+        (  1.94544667e-06_wp ) * ta3 * va2 +                                   &
+        (  1.14099531e-08_wp ) * ta4 * va2 +                                   &
+        (  1.58137256e-01_wp ) *       va3 +                                   &
+        ( -6.57263143e-05_wp ) * ta  * va3 +                                   &
+        (  2.22697524e-07_wp ) * ta2 * va3 +                                   &
+        ( -4.16117031e-08_wp ) * ta3 * va3 +                                   &
+        ( -1.27762753e-02_wp ) *       va4 +                                   &
+        (  9.66891875e-06_wp ) * ta  * va4 +                                   &
+        (  2.52785852e-09_wp ) * ta2 * va4 +                                   &
+        (  4.56306672e-04_wp ) *       va5 +                                   &
+        ( -1.74202546e-07_wp ) * ta  * va5 +                                   &
+        ( -5.91491269e-06_wp ) * va6
+
+    part_d_tmrt = (  3.98374029e-01_wp ) *       d_tmrt +                      &
+            (  1.83945314e-04_wp ) * ta  *       d_tmrt +                      &
+            ( -1.73754510e-04_wp ) * ta2 *       d_tmrt +                      &
+            ( -7.60781159e-07_wp ) * ta3 *       d_tmrt +                      &
+            (  3.77830287e-08_wp ) * ta4 *       d_tmrt +                      &
+            (  5.43079673e-10_wp ) * ta5 *       d_tmrt +                      &
+            ( -2.00518269e-02_wp ) *       va  * d_tmrt +                      &
+            (  8.92859837e-04_wp ) * ta  * va  * d_tmrt +                      &
+            (  3.45433048e-06_wp ) * ta2 * va  * d_tmrt +                      &
+            ( -3.77925774e-07_wp ) * ta3 * va  * d_tmrt +                      &
+            ( -1.69699377e-09_wp ) * ta4 * va  * d_tmrt +                      &
+            (  1.69992415e-04_wp ) *       va2 * d_tmrt +                      &
+            ( -4.99204314e-05_wp ) * ta  * va2 * d_tmrt +                      &
+            (  2.47417178e-07_wp ) * ta2 * va2 * d_tmrt +                      &
+            (  1.07596466e-08_wp ) * ta3 * va2 * d_tmrt +                      &
+            (  8.49242932e-05_wp ) *       va3 * d_tmrt +                      &
+            (  1.35191328e-06_wp ) * ta  * va3 * d_tmrt +                      &
+            ( -6.21531254e-09_wp ) * ta2 * va3 * d_tmrt +                      &
+            ( -4.99410301e-06_wp ) * va4 *       d_tmrt +                      &
+            ( -1.89489258e-08_wp ) * ta  * va4 * d_tmrt +                      &
+            (  8.15300114e-08_wp ) *       va5 * d_tmrt +                      &
+            (  7.55043090e-04_wp ) *             d_tmrt2 +                     &
+            ( -5.65095215e-05_wp ) * ta  *       d_tmrt2 +                     &
+            ( -4.52166564e-07_wp ) * ta2 *       d_tmrt2 +                     &
+            (  2.46688878e-08_wp ) * ta3 *       d_tmrt2 +                     &
+            (  2.42674348e-10_wp ) * ta4 *       d_tmrt2 +                     &
+            (  1.54547250e-04_wp ) *       va  * d_tmrt2 +                     &
+            (  5.24110970e-06_wp ) * ta  * va  * d_tmrt2 +                     &
+            ( -8.75874982e-08_wp ) * ta2 * va  * d_tmrt2 +                     &
+            ( -1.50743064e-09_wp ) * ta3 * va  * d_tmrt2 +                     &
+            ( -1.56236307e-05_wp ) *       va2 * d_tmrt2 +                     &
+            ( -1.33895614e-07_wp ) * ta  * va2 * d_tmrt2 +                     &
+            (  2.49709824e-09_wp ) * ta2 * va2 * d_tmrt2 +                     &
+            (  6.51711721e-07_wp ) *       va3 * d_tmrt2 +                     &
+            (  1.94960053e-09_wp ) * ta  * va3 * d_tmrt2 +                     &
+            ( -1.00361113e-08_wp ) *       va4 * d_tmrt2 +                     &
+            ( -1.21206673e-05_wp ) *             d_tmrt3 +                     &
+            ( -2.18203660e-07_wp ) * ta  *       d_tmrt3 +                     &
+            (  7.51269482e-09_wp ) * ta2 *       d_tmrt3 +                     &
+            (  9.79063848e-11_wp ) * ta3 *       d_tmrt3 +                     &
+            (  1.25006734e-06_wp ) *       va  * d_tmrt3 +                     &
+            ( -1.81584736e-09_wp ) * ta  * va  * d_tmrt3 +                     &
+            ( -3.52197671e-10_wp ) * ta2 * va  * d_tmrt3 +                     &
+            ( -3.36514630e-08_wp ) *       va2 * d_tmrt3 +                     &
+            (  1.35908359e-10_wp ) * ta  * va2 * d_tmrt3 +                     &
+            (  4.17032620e-10_wp ) *       va3 * d_tmrt3 +                     &
+            ( -1.30369025e-09_wp ) *             d_tmrt4 +                     &
+            (  4.13908461e-10_wp ) * ta  *       d_tmrt4 +                     &
+            (  9.22652254e-12_wp ) * ta2 *       d_tmrt4 +                     &
+            ( -5.08220384e-09_wp ) *       va  * d_tmrt4 +                     &
+            ( -2.24730961e-11_wp ) * ta  * va  * d_tmrt4 +                     &
+            (  1.17139133e-10_wp ) *       va2 * d_tmrt4 +                     &
+            (  6.62154879e-10_wp ) *             d_tmrt5 +                     &
+            (  4.03863260e-13_wp ) * ta  *       d_tmrt5 +                     &
+            (  1.95087203e-12_wp ) *       va  * d_tmrt5 +                     &
+            ( -4.73602469e-12_wp ) *             d_tmrt6
+
+    part_pa = (  5.12733497e+00_wp ) *                pa +                     &
+       ( -3.12788561e-01_wp ) * ta  *                 pa +                     &
+       ( -1.96701861e-02_wp ) * ta2 *                 pa +                     &
+       (  9.99690870e-04_wp ) * ta3 *                 pa +                     &
+       (  9.51738512e-06_wp ) * ta4 *                 pa +                     &
+       ( -4.66426341e-07_wp ) * ta5 *                 pa +                     &
+       (  5.48050612e-01_wp ) *       va  *           pa +                     &
+       ( -3.30552823e-03_wp ) * ta  * va  *           pa +                     &
+       ( -1.64119440e-03_wp ) * ta2 * va  *           pa +                     &
+       ( -5.16670694e-06_wp ) * ta3 * va  *           pa +                     &
+       (  9.52692432e-07_wp ) * ta4 * va  *           pa +                     &
+       ( -4.29223622e-02_wp ) *       va2 *           pa +                     &
+       (  5.00845667e-03_wp ) * ta  * va2 *           pa +                     &
+       (  1.00601257e-06_wp ) * ta2 * va2 *           pa +                     &
+       ( -1.81748644e-06_wp ) * ta3 * va2 *           pa +                     &
+       ( -1.25813502e-03_wp ) *       va3 *           pa +                     &
+       ( -1.79330391e-04_wp ) * ta  * va3 *           pa +                     &
+       (  2.34994441e-06_wp ) * ta2 * va3 *           pa +                     &
+       (  1.29735808e-04_wp ) *       va4 *           pa +                     &
+       (  1.29064870e-06_wp ) * ta  * va4 *           pa +                     &
+       ( -2.28558686e-06_wp ) *       va5 *           pa +                     &
+       ( -3.69476348e-02_wp ) *             d_tmrt  * pa +                     &
+       (  1.62325322e-03_wp ) * ta  *       d_tmrt  * pa +                     &
+       ( -3.14279680e-05_wp ) * ta2 *       d_tmrt  * pa +                     &
+       (  2.59835559e-06_wp ) * ta3 *       d_tmrt  * pa +                     &
+       ( -4.77136523e-08_wp ) * ta4 *       d_tmrt  * pa +                     &
+       (  8.64203390e-03_wp ) *       va  * d_tmrt  * pa +                     &
+       ( -6.87405181e-04_wp ) * ta  * va  * d_tmrt  * pa +                     &
+       ( -9.13863872e-06_wp ) * ta2 * va  * d_tmrt  * pa +                     &
+       (  5.15916806e-07_wp ) * ta3 * va  * d_tmrt  * pa +                     &
+       ( -3.59217476e-05_wp ) *       va2 * d_tmrt  * pa +                     &
+       (  3.28696511e-05_wp ) * ta  * va2 * d_tmrt  * pa +                     &
+       ( -7.10542454e-07_wp ) * ta2 * va2 * d_tmrt  * pa +                     &
+       ( -1.24382300e-05_wp ) *       va3 * d_tmrt  * pa +                     &
+       ( -7.38584400e-09_wp ) * ta  * va3 * d_tmrt  * pa +                     &
+       (  2.20609296e-07_wp ) *       va4 * d_tmrt  * pa +                     &
+       ( -7.32469180e-04_wp ) *             d_tmrt2 * pa +                     &
+       ( -1.87381964e-05_wp ) * ta  *       d_tmrt2 * pa +                     &
+       (  4.80925239e-06_wp ) * ta2 *       d_tmrt2 * pa +                     &
+       ( -8.75492040e-08_wp ) * ta3 *       d_tmrt2 * pa +                     &
+       (  2.77862930e-05_wp ) *       va  * d_tmrt2 * pa +                     &
+       ( -5.06004592e-06_wp ) * ta  * va  * d_tmrt2 * pa +                     &
+       (  1.14325367e-07_wp ) * ta2 * va  * d_tmrt2 * pa +                     &
+       (  2.53016723e-06_wp ) *       va2 * d_tmrt2 * pa +                     &
+       ( -1.72857035e-08_wp ) * ta  * va2 * d_tmrt2 * pa +                     &
+       ( -3.95079398e-08_wp ) *       va3 * d_tmrt2 * pa +                     &
+       ( -3.59413173e-07_wp ) *             d_tmrt3 * pa +                     &
+       (  7.04388046e-07_wp ) * ta  *       d_tmrt3 * pa +                     &
+       ( -1.89309167e-08_wp ) * ta2 *       d_tmrt3 * pa +                     &
+       ( -4.79768731e-07_wp ) *       va  * d_tmrt3 * pa +                     &
+       (  7.96079978e-09_wp ) * ta  * va  * d_tmrt3 * pa +                     &
+       (  1.62897058e-09_wp ) *       va2 * d_tmrt3 * pa +                     &
+       (  3.94367674e-08_wp ) *             d_tmrt4 * pa +                     &
+       ( -1.18566247e-09_wp ) * ta *        d_tmrt4 * pa +                     &
+       (  3.34678041e-10_wp ) *       va  * d_tmrt4 * pa +                     &
+       ( -1.15606447e-10_wp ) *             d_tmrt5 * pa
+
+    part_pa2 = ( -2.80626406e+00_wp ) *               pa2 +                    &
+       (  5.48712484e-01_wp ) * ta  *                 pa2 +                    &
+       ( -3.99428410e-03_wp ) * ta2 *                 pa2 +                    &
+       ( -9.54009191e-04_wp ) * ta3 *                 pa2 +                    &
+       (  1.93090978e-05_wp ) * ta4 *                 pa2 +                    &
+       ( -3.08806365e-01_wp ) *       va *            pa2 +                    &
+       (  1.16952364e-02_wp ) * ta  * va *            pa2 +                    &
+       (  4.95271903e-04_wp ) * ta2 * va *            pa2 +                    &
+       ( -1.90710882e-05_wp ) * ta3 * va *            pa2 +                    &
+       (  2.10787756e-03_wp ) *       va2 *           pa2 +                    &
+       ( -6.98445738e-04_wp ) * ta  * va2 *           pa2 +                    &
+       (  2.30109073e-05_wp ) * ta2 * va2 *           pa2 +                    &
+       (  4.17856590e-04_wp ) *       va3 *           pa2 +                    &
+       ( -1.27043871e-05_wp ) * ta  * va3 *           pa2 +                    &
+       ( -3.04620472e-06_wp ) *       va4 *           pa2 +                    &
+       (  5.14507424e-02_wp ) *             d_tmrt  * pa2 +                    &
+       ( -4.32510997e-03_wp ) * ta  *       d_tmrt  * pa2 +                    &
+       (  8.99281156e-05_wp ) * ta2 *       d_tmrt  * pa2 +                    &
+       ( -7.14663943e-07_wp ) * ta3 *       d_tmrt  * pa2 +                    &
+       ( -2.66016305e-04_wp ) *       va  * d_tmrt  * pa2 +                    &
+       (  2.63789586e-04_wp ) * ta  * va  * d_tmrt  * pa2 +                    &
+       ( -7.01199003e-06_wp ) * ta2 * va  * d_tmrt  * pa2 +                    &
+       ( -1.06823306e-04_wp ) *       va2 * d_tmrt  * pa2 +                    &
+       (  3.61341136e-06_wp ) * ta  * va2 * d_tmrt  * pa2 +                    &
+       (  2.29748967e-07_wp ) *       va3 * d_tmrt  * pa2 +                    &
+       (  3.04788893e-04_wp ) *             d_tmrt2 * pa2 +                    &
+       ( -6.42070836e-05_wp ) * ta  *       d_tmrt2 * pa2 +                    &
+       (  1.16257971e-06_wp ) * ta2 *       d_tmrt2 * pa2 +                    &
+       (  7.68023384e-06_wp ) *       va  * d_tmrt2 * pa2 +                    &
+       ( -5.47446896e-07_wp ) * ta  * va  * d_tmrt2 * pa2 +                    &
+       ( -3.59937910e-08_wp ) *       va2 * d_tmrt2 * pa2 +                    &
+       ( -4.36497725e-06_wp ) *             d_tmrt3 * pa2 +                    &
+       (  1.68737969e-07_wp ) * ta  *       d_tmrt3 * pa2 +                    &
+       (  2.67489271e-08_wp ) *       va  * d_tmrt3 * pa2 +                    &
+       (  3.23926897e-09_wp ) *             d_tmrt4 * pa2
+
+    part_pa3 = ( -3.53874123e-02_wp ) *               pa3 +                    &
+       ( -2.21201190e-01_wp ) * ta  *                 pa3 +                    &
+       (  1.55126038e-02_wp ) * ta2 *                 pa3 +                    &
+       ( -2.63917279e-04_wp ) * ta3 *                 pa3 +                    &
+       (  4.53433455e-02_wp ) *       va  *           pa3 +                    &
+       ( -4.32943862e-03_wp ) * ta  * va  *           pa3 +                    &
+       (  1.45389826e-04_wp ) * ta2 * va  *           pa3 +                    &
+       (  2.17508610e-04_wp ) *       va2 *           pa3 +                    &
+       ( -6.66724702e-05_wp ) * ta  * va2 *           pa3 +                    &
+       (  3.33217140e-05_wp ) *       va3 *           pa3 +                    &
+       ( -2.26921615e-03_wp ) *             d_tmrt  * pa3 +                    &
+       (  3.80261982e-04_wp ) * ta  *       d_tmrt  * pa3 +                    &
+       ( -5.45314314e-09_wp ) * ta2 *       d_tmrt  * pa3 +                    &
+       ( -7.96355448e-04_wp ) *       va  * d_tmrt  * pa3 +                    &
+       (  2.53458034e-05_wp ) * ta  * va  * d_tmrt  * pa3 +                    &
+       ( -6.31223658e-06_wp ) *       va2 * d_tmrt  * pa3 +                    &
+       (  3.02122035e-04_wp ) *             d_tmrt2 * pa3 +                    &
+       ( -4.77403547e-06_wp ) * ta  *       d_tmrt2 * pa3 +                    &
+       (  1.73825715e-06_wp ) *       va  * d_tmrt2 * pa3 +                    &
+       ( -4.09087898e-07_wp ) *             d_tmrt3 * pa3
+
+    part_pa46 = (  6.14155345e-01_wp ) *              pa4 +                    &
+       ( -6.16755931e-02_wp ) * ta  *                 pa4 +                    &
+       (  1.33374846e-03_wp ) * ta2 *                 pa4 +                    &
+       (  3.55375387e-03_wp ) *       va  *           pa4 +                    &
+       ( -5.13027851e-04_wp ) * ta  * va *            pa4 +                    &
+       (  1.02449757e-04_wp ) *       va2 *           pa4 +                    &
+       ( -1.48526421e-03_wp ) *             d_tmrt  * pa4 +                    &
+       ( -4.11469183e-05_wp ) * ta  *       d_tmrt  * pa4 +                    &
+       ( -6.80434415e-06_wp ) *       va  * d_tmrt  * pa4 +                    &
+       ( -9.77675906e-06_wp ) *             d_tmrt2 * pa4 +                    &
+       (  8.82773108e-02_wp ) *                       pa5 +                    &
+       ( -3.01859306e-03_wp ) * ta  *                 pa5 +                    &
+       (  1.04452989e-03_wp ) *       va  *           pa5 +                    &
+       (  2.47090539e-04_wp ) *             d_tmrt  * pa5 +                    &
+       (  1.48348065e-03_wp ) *                       pa6
+!
 !-- 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 
-
+    utci_ij = ta + part_ta + part_va + part_d_tmrt + part_pa + part_pa2 +      &
+        part_pa3 + part_pa46
+!
 !-- Consider offset in result
-    utci = utci + offset
+    utci_ij = utci_ij + offset
 
  END SUBROUTINE calculate_utci_static
@@ -1570,8 +1826,8 @@
 !> Value of perct is the Perceived Temperature, degree centigrade
 !------------------------------------------------------------------------------!
- SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct )
+ SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct_ij )
 
     IMPLICIT NONE
-
+!
 !-- Type of input of the argument list
     REAL(wp), INTENT ( IN )  :: ta   !< Local air temperature (degC)
@@ -1580,13 +1836,13 @@
     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)
-
+    REAL(wp), INTENT ( OUT ) :: perct_ij  !< 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)
@@ -1614,20 +1870,23 @@
 
     LOGICAL :: sultrieness
-
+!
 !-- Initialise
-    perct = 9999.0_wp
+    perct_ij = bio_fill_value
 
     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
@@ -1636,5 +1895,5 @@
 
        IF ( pmva > 0._wp ) THEN
-
+!
 !--       Case summer clothing insulation: heat load ?
           clo = sclo
@@ -1643,6 +1902,7 @@
           pmv_s = pmva
           IF ( pmva <= 0._wp ) THEN
+!
 !--          Case: comfort achievable by varying clothing insulation 
-!--          Between winter and summer set values
+!            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 )
@@ -1662,5 +1922,5 @@
        ENDIF
     ELSE
-
+!
 !--    First guess: summer clothing insulation: heat load
        clo = sclo
@@ -1669,5 +1929,5 @@
 
        IF ( pmva < 0._wp ) THEN
-
+!
 !--       Case winter clothing insulation: cold stress ?
           clo = wclo
@@ -1677,5 +1937,5 @@
 
           IF ( pmva >= 0._wp ) THEN
-
+!
 !--          Case: comfort achievable by varying clothing insulation 
 !            between winter and summer set values
@@ -1698,10 +1958,11 @@
 
     ENDIF
-
+!
 !-- Determine perceived temperature by regression equation + adjustments
     pmvs = pmva
-    CALL perct_regression ( pmva, clo, perct )
-    ptc = perct
+    CALL perct_regression ( pmva, clo, perct_ij )
+    ptc = perct_ij
     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 )
@@ -1712,13 +1973,14 @@
           pmvs   = -0.11_wp
        ENDIF
-       CALL perct_regression ( pmvs, clo, perct )
+       CALL perct_regression ( pmvs, clo, perct_ij )
     ENDIF
 !     clo_fanger = clo
     clon = clo
-    IF ( clo > 0.5_wp .AND. perct <= 8.73_wp ) THEN
+    IF ( clo > 0.5_wp .AND. perct_ij <= 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 )
+       clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
        clo = clon
     ENDIF
@@ -1727,14 +1989,16 @@
     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 )
+       CALL perct_regression ( pmvs, clo, perct_ij )
        IF ( sult_lim < 99._wp ) THEN
-          IF ( (perct - ptc) > sult_lim ) sultrieness = .TRUE.
+          IF ( (perct_ij - ptc) > sult_lim ) sultrieness = .TRUE.
+!
 !--       Set factor to threshold for sultriness
           IF ( dgtcstd /= 0_iwp ) THEN
-             d_std = ( ( perct - ptc ) - dgtcm ) / dgtcstd
+             d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
           ENDIF
        ENDIF
@@ -1762,13 +2026,15 @@
 
     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 ) )
@@ -1788,5 +2054,5 @@
 
     IMPLICIT NONE
-
+!
 !-- Input variables of argument list:
     REAL(wp), INTENT ( IN )  :: ta     !< Ambient temperature (degC)
@@ -1803,5 +2069,5 @@
     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)
@@ -1814,5 +2080,5 @@
 !           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
@@ -1831,9 +2097,10 @@
     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
+    x_ridder    = bio_fill_value
     pmva        = 0._wp
     clo_lower   = sclo
@@ -1884,4 +2151,5 @@
              y_lower = y_new
           ELSE
+!
 !--          Never get here in x_ridder: singularity in y
              nerr       = -1_iwp
@@ -1895,4 +2163,5 @@
           ENDIF
        ENDDO
+!
 !--    x_ridder exceed maximum iterations
        nerr       = -2_iwp
@@ -1904,4 +2173,5 @@
        x_ridder = clo_upper
     ELSE
+!
 !--    x_ridder not bracketed by u_clo and o_clo
        nerr = -3_iwp
@@ -1919,5 +2189,5 @@
 !> individual and reference environment.
 !------------------------------------------------------------------------------!
- SUBROUTINE perct_regression( pmv, clo, perct )
+ SUBROUTINE perct_regression( pmv, clo, perct_ij )
 
     IMPLICIT NONE
@@ -1926,13 +2196,13 @@
     REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
 
-    REAL(wp), INTENT ( OUT ) ::  perct   !< perct (degC) corresponding to given PMV / clo
+    REAL(wp), INTENT ( OUT ) ::  perct_ij   !< perct (degC) corresponding to given PMV / clo
 
     IF ( pmv <= -0.11_wp ) THEN
-       perct = 5.805_wp + 12.6784_wp * pmv
+       perct_ij = 5.805_wp + 12.6784_wp * pmv
     ELSE
        IF ( pmv >= + 0.01_wp ) THEN
-          perct = 16.826_wp + 6.163_wp * pmv
+          perct_ij = 16.826_wp + 6.163_wp * pmv
        ELSE
-          perct = 21.258_wp - 9.558_wp * clo
+          perct_ij = 21.258_wp - 9.558_wp * clo
        ENDIF
     ENDIF
@@ -1953,5 +2223,5 @@
 
     IMPLICIT NONE
-
+!
 !-- Input variables of argument list:
     REAL(wp), INTENT ( IN ) ::  ta       !< Ambient air temperature (degC)
@@ -1963,5 +2233,5 @@
     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, 
@@ -1969,5 +2239,5 @@
 !            (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)
@@ -1992,12 +2262,12 @@
     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
@@ -2005,14 +2275,14 @@
        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
@@ -2021,5 +2291,5 @@
     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
@@ -2029,19 +2299,19 @@
           * 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 +        &
@@ -2052,11 +2322,11 @@
        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 )
@@ -2073,5 +2343,5 @@
 
     IMPLICIT NONE
-
+!
 !-- Input variables of argument list:
     REAL(wp),     INTENT ( IN )  :: pmva     !< Actual Predicted Mean Vote (PMV) according to Fanger
@@ -2081,5 +2351,5 @@
     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
@@ -2088,5 +2358,5 @@
 !            -3 = rel. humidity set to 5 % or 95 %, respectively
 !            -4 = deltapmv set to avoid pmvs < 0
-
+!
 !-- Internal variable types:
     REAL(wp) ::  pmv          !< 
@@ -2142,5 +2412,5 @@
      +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
@@ -2149,8 +2419,8 @@
        nerr = 0_iwp
     ENDIF
-
+!
 !-- Initialise classic PMV
     pmv  = pmva
-
+!
 !-- Water vapour pressure of air 
     p10  = 0.05_wp * svp_ta
@@ -2161,7 +2431,9 @@
        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
@@ -2169,4 +2441,5 @@
     ENDIF
     IF ( pa > 0._wp ) THEN
+!
 !--    Natural logarithm of pa
        apa = LOG ( pa )
@@ -2174,5 +2447,5 @@
        apa = -5._wp
     ENDIF
-
+!
 !-- Ratio actual water vapour pressure to that of a r.H. of 50 %
     pa_p50   = 0.5_wp * svp_ta
@@ -2184,4 +2457,5 @@
        pa_p50 = 0._wp
     ENDIF
+!
 !-- Square root of wind velocity
     IF ( ws >= 0._wp ) THEN
@@ -2190,12 +2464,14 @@
        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
@@ -2213,5 +2489,5 @@
        ENDIF
     ENDIF
-
+!
 !-- Regression valid for 0. <= pmv <= 6.
     dpmv_1 =                                                                   &
@@ -2239,15 +2515,18 @@
           + 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
@@ -2266,11 +2545,11 @@
 !> perct. 
 !------------------------------------------------------------------------------!
- SUBROUTINE calc_sultr( perct, dperctm, dperctstd, sultr_res )
+ SUBROUTINE calc_sultr( perct_ij, dperctm, dperctstd, sultr_res )
 
     IMPLICIT NONE
-
+!
 !-- Input of the argument list:
-    REAL(wp), INTENT ( IN )  ::  perct      !< Classical perceived temperature: Base is Fanger's PMV
-
+    REAL(wp), INTENT ( IN )  ::  perct_ij      !< 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 
@@ -2279,5 +2558,5 @@
 !             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
@@ -2285,5 +2564,5 @@
     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
@@ -2292,9 +2571,9 @@
     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
@@ -2302,15 +2581,18 @@
     dperctstd = 999999._wp
 
-    IF ( perct < 16.826_wp .OR. perct > 56._wp ) THEN
-!--    Unallowed classical PMV/perct
+    IF ( perct_ij < 16.826_wp .OR. perct_ij > 56._wp ) THEN
+!
+!--    Unallowed value of classical perct!
        RETURN
     ENDIF
-
+!
 !-- Mean deviation dependent on perct
-    dperctm = dperctka + dperctkb * perct + dperctkc * perct**2._wp + dperctkd * perct**3._wp
-
+    dperctm = dperctka + dperctkb * perct_ij + dperctkc * perct_ij**2._wp +    &
+       dperctkd * perct_ij**3._wp
+!
 !-- Mean deviation plus its standard deviation
-    dperctstd = dperctsa + dperctsb * perct + dperctsc * perct**2._wp + dperctsd * perct**3._wp
-
+    dperctstd = dperctsa + dperctsb * perct_ij + dperctsc * perct_ij**2._wp +  &
+       dperctsd * perct_ij**3._wp
+!
 !-- Value of the FUNCTION
     sultr_res = dperctm + faktor * dperctstd
@@ -2330,5 +2612,5 @@
 
     IMPLICIT NONE
-
+!
 !-- Type of input arguments
     REAL(wp), INTENT ( IN ) ::  pmva   !< Fanger's classical predicted mean vote
@@ -2336,10 +2618,10 @@
     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)
@@ -2352,9 +2634,9 @@
     REAL(wp) ::  sqrt_ws
     INTEGER(iwp) ::  i
-    INTEGER(iwp) ::  j
+!     INTEGER(iwp) ::  j
     INTEGER(iwp) ::  i_bin
-
+!
 !-- Coefficient of the 3 regression lines
-    !     1:const   2:*pmvc  3:*ta      4:*sqrt_ws  5:*dtmrt
+!     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 /)
@@ -2363,5 +2645,5 @@
     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
@@ -2384,10 +2666,11 @@
 
     DO i = 1, 3
+!
 !--    Regression constant for given meteorological variables
-       reg_a(i) = coeff(i, 1) + coeff(i, 3) * ta + coeff(i, 4) *             &
+       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
@@ -2408,4 +2691,5 @@
        ENDIF
     ENDDO
+!
 !-- Adjust to operative temperature scaled according 
 !   to classical PMV (Fanger)
@@ -2434,5 +2718,5 @@
     INTEGER(iwp)  ::  i, i_bin
 
-!                             range_1     range_2     range_3
+!                                   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/
@@ -2467,25 +2751,25 @@
 !> for perct < 10 (degC)
 !------------------------------------------------------------------------------!
- REAL(wp) FUNCTION ireq_neutral( perct, ireq_minimal, nerr )
+ REAL(wp) FUNCTION ireq_neutral( perct_ij, ireq_minimal, nerr )
 
     IMPLICIT NONE
-
+!
 !-- Type declaration of arguments
-    REAL(wp),     INTENT ( IN )  ::  perct
+    REAL(wp),     INTENT ( IN )  ::  perct_ij
     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
-
+    top02 = 1.8788_wp + 0.9296_wp * perct_ij
+!
 !-- 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
@@ -2495,5 +2779,5 @@
        ireq_neutral = 0.5_wp
     ENDIF
-
+!
 !-- Minimal required clothing insulation: maximal acceptable body cooling
     ireq_minimal = 1.26_wp - 0.0588_wp * top02
@@ -2523,5 +2807,5 @@
 
     height = height_cm * 100._wp
-
+!
 !-- According to Gehan-George, for children
     IF ( age < 19_iwp ) THEN
@@ -2533,5 +2817,5 @@
        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.
@@ -2598,9 +2882,9 @@
 
  SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo,             &
-     ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                         &
+     ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                          &
      ipt )
 
     IMPLICIT NONE
-
+!
 !-- Input parameters
     REAL(wp), INTENT(in) ::  age        !< Persons age          (years)
@@ -2615,5 +2899,5 @@
     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
@@ -2622,5 +2906,5 @@
     REAL(wp), INTENT(out) ::  t_clothing
     REAL(wp), INTENT(out) ::  ipt
-
+!
 !-- Internal variables
     REAL(wp), PARAMETER :: eps = 0.0005_wp
@@ -2644,5 +2928,5 @@
     REAL(wp) ::  top
     REAL(wp) ::  a_surf
-    REAL(wp) ::  acti
+!     REAL(wp) ::  acti
     INTEGER(iwp) ::  ncount
     INTEGER(iwp) ::  nerr_cold
@@ -2653,11 +2937,12 @@
     storage = 0._wp
     CALL persdat ( age, weight, height, sex, work, a_surf, actlev )
-
+!
 !-- Initialise
-    t_clothing = -999.0_wp
-    ipt        = -999.0_wp
+    t_clothing = bio_fill_value
+    ipt        = bio_fill_value
     nerr       = 0_wp
     ncount     = 0_wp
     sultrieness    = .FALSE.
+!
 !-- Tresholds: clothing insulation (account for model inaccuracies)
 !-- Summer clothing
@@ -2665,11 +2950,11 @@
 !-- 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
+       t_clothing = bio_fill_value  ! force initial run
        CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
           t_clothing, storage, dt, pmva )
@@ -2677,12 +2962,13 @@
 
        IF ( pmva > 0._wp ) THEN
-
+!
 !--       Case summer clothing insulation: heat load ?            
           clo = sclo
-          t_clothing = -999.0_wp  ! force initial run
+          t_clothing = bio_fill_value  ! 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
@@ -2695,5 +2981,5 @@
           ELSE IF ( pmva > 0.06_wp ) THEN
              clo = 0.5_wp
-             t_clothing = -999.0_wp
+             t_clothing = bio_fill_value
              CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
                 t_clothing, storage, dt, pmva )
@@ -2701,5 +2987,5 @@
        ELSE IF ( pmva < -0.11_wp ) THEN
           clo = 1.75_wp
-          t_clothing = -999.0_wp
+          t_clothing = bio_fill_value
           CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
              t_clothing, storage, dt, pmva )
@@ -2707,8 +2993,8 @@
 
     ELSE
-
+!
 !--    First guess: summer clothing insulation: heat load
        clo = sclo
-       t_clothing = -999.0_wp
+       t_clothing = bio_fill_value
        CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
           t_clothing, storage, dt, pmva )
@@ -2716,8 +3002,8 @@
 
        IF ( pmva < 0._wp ) THEN
-
+!
 !--       Case winter clothing insulation: cold stress ?
           clo = wclo
-          t_clothing = -999.0_wp
+          t_clothing = bio_fill_value
           CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
              t_clothing, storage, dt, pmva )
@@ -2725,5 +3011,5 @@
 
           IF ( pmva >= 0._wp ) THEN
-
+!
 !--          Case: comfort achievable by varying clothing insulation 
 !            between winter and summer set values
@@ -2736,5 +3022,5 @@
           ELSE IF ( pmva < -0.11_wp ) THEN
              clo = 1.75_wp
-             t_clothing = -999.0_wp
+             t_clothing = bio_fill_value
              CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
                 t_clothing, storage, dt, pmva )
@@ -2742,5 +3028,5 @@
        ELSE IF ( pmva > 0.06_wp ) THEN
           clo = 0.5_wp
-          t_clothing = -999.0_wp
+          t_clothing = bio_fill_value
           CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
              t_clothing, storage, dt, pmva )
@@ -2748,5 +3034,5 @@
 
     ENDIF
-
+!
 !-- Determine perceived temperature by regression equation + adjustments
     pmvs = pmva
@@ -2754,4 +3040,5 @@
     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 )
@@ -2767,4 +3054,5 @@
     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 
@@ -2777,4 +3065,5 @@
     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 )
@@ -2808,5 +3097,5 @@
 
     IMPLICIT NONE
-
+!
 !-- Type of input of the argument list
     REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature             (Â°C)
@@ -2819,12 +3108,12 @@
     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
@@ -2841,19 +3130,20 @@
 
     LOGICAL ::  sultrieness
-
+!
 !-- Initialise
-    ipt = -999.0_wp
+    ipt = bio_fill_value
 
     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 )
@@ -2866,5 +3156,5 @@
        CALL perct_regression ( pmvs, clo, ipt )
     ENDIF
-
+!
 !-- Consider sultriness if appropriate
     ptc = ipt
@@ -2873,4 +3163,5 @@
     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 )
@@ -2898,5 +3189,5 @@
 
     IMPLICIT NONE
-
+!
 !--  Input argument types
     REAL(wp), INTENT ( IN )  ::  ta       !< Air temperature          (Â°C)
@@ -2909,5 +3200,5 @@
     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)
@@ -2915,5 +3206,5 @@
     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
@@ -2928,5 +3219,5 @@
     REAL(wp) ::  gc           !< preliminary result storage
     REAL(wp) ::  t_clothing   !< clothing temperature                   (Â°C)
-    REAL(wp) ::  hr           !< radiational heat resistence
+!     REAL(wp) ::  hr           !< radiational heat resistence
     REAL(wp) ::  clo          !< clothing insulation index              (clo)
     REAL(wp) ::  ws           !< wind speed                             (m/s)
@@ -2946,12 +3237,12 @@
     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
@@ -2959,11 +3250,11 @@
        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
@@ -2972,8 +3263,9 @@
     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
+    niter = INT( dt * 10._wp, KIND=iwp )
+    IF ( niter < 1 ) niter = 1_iwp
     adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
     IF ( adjustrate >= 1._wp ) adjustrate = 1._wp
@@ -2982,5 +3274,5 @@
 
     IF ( t_cloth <= -998.0_wp ) THEN  ! If initial run
-       niter = 3_wp
+       niter = 3_iwp
        adjustrate = 1._wp
        adjustrate_cloth = 1._wp
@@ -2993,35 +3285,35 @@
           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
@@ -3041,8 +3333,8 @@
 !------------------------------------------------------------------------------!
 
- SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet )
+ SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet_ij )
 
     IMPLICIT NONE
-
+!
 !-- Input arguments:
     REAL(wp), INTENT( IN ) ::  ta    !< Air temperature             (Â°C)
@@ -3051,8 +3343,8 @@
     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)
-
+    REAL(wp), INTENT ( OUT ) ::  pet_ij  !< PET                     (Â°C)
+!
 !-- Internal variables:
     REAL(wp) ::  acl        !< clothing area                        (mÂ²)
@@ -3071,8 +3363,8 @@
     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)
@@ -3085,5 +3377,5 @@
 !     INTEGER(iwp) :: pos     !< Posture: 1 = standing, 2 = sitting
 !     INTEGER(iwp) :: sex     !< Sex: 1 = male, 2 = female
-
+!
 !-- Configuration, keep standard parameters!
     age   = 35._wp
@@ -3094,5 +3386,5 @@
     clo   =  0.9_wp
     fcl   =  1.15_wp
-
+!
 !-- Call subfunctions
     CALL in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,    &
@@ -3103,6 +3395,6 @@
             vpts, wetsk )
 
-    CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair, rdcl,&
-            rdsk, rtv, ta, tcl, tsk, pet, vpts, wetsk )
+    CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,      &
+            rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
 
 
@@ -3117,5 +3409,5 @@
  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)
@@ -3127,5 +3419,5 @@
     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)
@@ -3133,9 +3425,9 @@
     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)
@@ -3151,14 +3443,14 @@
 
     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
@@ -3176,5 +3468,5 @@
         vpts, wetsk )
 
-
+!
 !-- Input arguments:
     REAL(wp), INTENT( IN )  ::  ere    !< Energy ballance          (W)
@@ -3190,5 +3482,5 @@
     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Â²)
@@ -3204,5 +3496,5 @@
     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
@@ -3214,5 +3506,5 @@
     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)
@@ -3239,5 +3531,5 @@
     REAL(wp) ::  rsum           !< Radiational loss or gain         (W/mÂ²)
     REAL(wp) ::  sw             !< 
-    REAL(wp) ::  swf            !< 
+!     REAL(wp) ::  swf            !< female factor, currently unused
     REAL(wp) ::  swm            !< 
     REAL(wp) ::  tbody          !< 
@@ -3279,5 +3571,5 @@
 
     di   = r2 - r1
-
+!
 !-- Skin temperatur
     DO j = 1, 7
@@ -3295,5 +3587,5 @@
           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
@@ -3303,10 +3595,10 @@
             ( ( 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
@@ -3337,5 +3629,5 @@
              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)
@@ -3360,5 +3652,5 @@
           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
@@ -3366,5 +3658,5 @@
           ed   = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa -     &
              vpts )
-
+!
 !--       Max vb
           vb1 = 34._wp - tsk
@@ -3374,8 +3666,8 @@
           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
@@ -3429,5 +3721,5 @@
        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
@@ -3448,6 +3740,6 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,   &
-        rdcl, rdsk, rtv, ta, tcl, tsk, pet, vpts, wetsk )
-
+        rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
+!
 !-- Input arguments:
     REAL(wp), INTENT( IN ) ::  acl   !< clothing surface area        (mÂ²)
@@ -3466,16 +3758,14 @@
     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)
-
+    REAL(wp), INTENT( OUT ) ::  aeff     !< effective surface area       (mÂ²)
+    REAL(wp), INTENT( OUT ) ::  pet_ij   !< 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
@@ -3499,5 +3789,5 @@
     INTEGER ( iwp ) ::  i             !< running index
 
-    pet = ta
+    pet_ij = ta
     enbal2 = 0._wp
 
@@ -3506,39 +3796,40 @@
           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 )
+              ( ( pet_ij + 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 )
+              ( ( pet_ij + 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
+          cbare = hc * ( pet_ij - tsk ) * adu * ( 1._wp - facl )
+          cclo  = hc * ( pet_ij - 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
+          tex  = 0.47_wp * pet_ij + 21._wp
+          eres = c_p * ( pet_ij - 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
+          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 )  pet = pet - xx
-          IF ( enbal < 0._wp )  pet = pet + xx
+!           IF ( count1 == 3_iwp )  xx = 0.001_wp
+          IF ( enbal > 0._wp )  pet_ij = pet_ij - xx
+          IF ( enbal < 0._wp )  pet_ij = pet_ij + xx
           IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT
           IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT
@@ -3549,4 +3840,232 @@
  END SUBROUTINE pet_iteration
 
+!
+!-- UVEM specific subroutines
+
+!---------------------------------------------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Module-specific routine for new module
+!---------------------------------------------------------------------------------------------------------------------!
+ SUBROUTINE uvem_solar_position
+    
+    USE date_and_time_mod,                                                                                            &
+       ONLY:  calc_date_and_time, day_of_year, time_utc 
+    
+    USE control_parameters,                                                                                           &
+       ONLY:  latitude, longitude   
+
+    IMPLICIT NONE
+    
+   
+    REAL(wp) ::  alpha       = 0.0_wp   !< solar azimuth angle in radiant    
+    REAL(wp) ::  doy_r       = 0.0_wp   !< real format of day_of_year           
+    REAL(wp) ::  declination = 0.0_wp   !< declination
+    REAL(wp) ::  dtor        = 0.0_wp   !< factor to convert degree to radiant
+    REAL(wp) ::  js          = 0.0_wp   !< parameter for solar position calculation
+    REAL(wp) ::  lat         = 52.39_wp !< latitude
+    REAL(wp) ::  lon         = 9.7_wp   !< longitude        
+    REAL(wp) ::  thetar      = 0.0_wp   !< angle for solar zenith angle calculation
+    REAL(wp) ::  thetasr     = 0.0_wp   !< angle for solar azimuth angle calculation    
+    REAL(wp) ::  zgl         = 0.0_wp   !< calculated exposure by direct beam   
+    REAL(wp) ::  woz         = 0.0_wp   !< calculated exposure by diffuse radiation
+    REAL(wp) ::  wsp         = 0.0_wp   !< calculated exposure by direct beam   
+   
+
+    CALL calc_date_and_time
+    doy_r = real(day_of_year)   
+    dtor = pi / 180.0_wp
+    lat = latitude
+    lon = longitude
+!
+!-- calculation of js, necessary for calculation of equation of time (zgl) :
+    js=  72.0_wp * ( doy_r + ( time_utc / 86400.0_wp ) ) / 73.0_wp 
+!
+!-- calculation of equation of time (zgl):
+    zgl = 0.0066_wp + 7.3525_wp * cos( ( js + 85.9_wp ) * dtor ) + 9.9359_wp *                                        &
+    cos( ( 2.0_wp * js + 108.9_wp ) * dtor ) + 0.3387_wp * cos( ( 3 * js + 105.2_wp ) * dtor )
+!
+!-- calculation of apparent solar time woz:
+    woz = ( ( time_utc / 3600.0_wp ) - ( 4.0_wp * ( 15.0_wp - lon ) ) / 60.0_wp ) + ( zgl / 60.0_wp )
+!
+!-- calculation of hour angle (wsp):
+    wsp = ( woz - 12.0_wp ) * 15.0_wp
+!
+!-- calculation of declination:
+    declination = 0.3948_wp - 23.2559_wp * cos( ( js + 9.1_wp ) * dtor ) -                                            &
+    0.3915_wp * cos( ( 2.0_wp * js + 5.4_wp ) * dtor ) - 0.1764_wp * cos( ( 3.0_wp * js + 26.0_wp ) * dtor )
+!
+!-- calculation of solar zenith angle
+    thetar  = acos( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) *                                &
+    cos( lat * dtor ) * cos( declination * dtor ) )
+    thetasr = asin( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) *                                & 
+    cos( lat * dtor ) * cos( declination * dtor ) )
+    sza = thetar / dtor
+!
+!-- calculation of solar azimuth angle
+    IF (woz .LE. 12.0_wp) alpha = pi - acos( ( sin(thetasr) * sin( lat * dtor ) -                                     &
+    sin( declination * dtor ) ) / ( cos(thetasr) * cos( lat * dtor ) ) )    
+    IF (woz .GT. 12.0_wp) alpha = pi + acos( ( sin(thetasr) * sin( lat * dtor ) -                                     &
+    sin( declination * dtor ) ) / ( cos(thetasr) * cos( lat * dtor ) ) )    
+    saa = alpha / dtor
+
+ END SUBROUTINE uvem_solar_position
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Module-specific routine for new module
+!---------------------------------------------------------------------------------------------------------------------!
+ SUBROUTINE uvem_calc_exposure
+
+    USE indices,                                                                                                      &
+        ONLY:  nxlg, nxrg, nyng, nysg, nys, nyn, nxl, nxr 
+    
+    
+    IMPLICIT NONE   
+    
+    INTEGER(iwp) ::  i     !< loop index in x direction
+    INTEGER(iwp) ::  j     !< loop index in y direction
+    INTEGER(iwp) ::  szai  !< loop index for different sza values
+
+    CALL uvem_solar_position
+     
+    IF (sza  .GE.  90) THEN
+       vitd3_exposure(:,:) = 0.0_wp
+    ELSE
+       
+       DO  ai = 0, 35
+          DO  zi = 0, 9
+                projection_area_lookup_table(ai,zi) = uvem_projarea_f%var(clothing,zi,ai)
+          ENDDO
+       ENDDO
+       DO  ai = 0, 35
+          DO  zi = 0, 9
+                integration_array(ai,zi) = uvem_integration_f%var(zi,ai)
+          ENDDO
+       ENDDO
+       DO  ai = 0, 2
+          DO  zi = 0, 90
+                irradiance_lookup_table(ai,zi) = uvem_irradiance_f%var(zi,ai)
+          ENDDO
+       ENDDO
+       DO  ai = 0, 35
+          DO  zi = 0, 9
+             DO  szai = 0, 90
+                radiance_lookup_table(ai,zi,szai) = uvem_radiance_f%var(szai,zi,ai)
+             ENDDO
+          ENDDO
+       ENDDO
+        
+        
+        
+!--    rotate 3D-Model human to desired direction  ----------------------------- 
+       projection_area_temp( 0:35,:) = projection_area_lookup_table
+       projection_area_temp(36:71,:) = projection_area_lookup_table               
+       IF (  .NOT.  turn_to_sun ) startpos_human = orientation_angle / 10.0_wp
+       IF (       turn_to_sun ) startpos_human = saa / 10.0_wp        
+       DO  ai = 0, 35
+          xfactor = ( startpos_human ) - INT( startpos_human )
+          DO  zi = 0, 9
+             projection_area(ai,zi) = ( projection_area_temp( 36 - INT( startpos_human ) - 1 + ai , zi) *             &
+                                      ( xfactor ) )                                                                   &
+                                      +( projection_area_temp( 36 - INT( startpos_human ) + ai , zi) *                &
+                                      ( 1.0_wp - xfactor ) )
+          ENDDO
+       ENDDO           
+!             
+!           
+!--    interpolate to accurate Solar Zenith Angle  ------------------          
+       DO  ai = 0, 35
+          xfactor = (sza)-INT(sza)
+          DO  zi = 0, 9
+             radiance_array(ai,zi) = ( radiance_lookup_table(ai, zi, INT(sza) ) * ( 1.0_wp - xfactor) ) +             &
+             ( radiance_lookup_table(ai,zi,INT(sza) + 1) * xfactor )
+          ENDDO
+       ENDDO
+       DO  iq = 0, 2
+          irradiance(iq) = ( irradiance_lookup_table(iq, INT(sza) ) * ( 1.0_wp - xfactor)) +                          &
+          (irradiance_lookup_table(iq, INT(sza) + 1) * xfactor )
+       ENDDO   
+!         
+!--    interpolate to accurate Solar Azimuth Angle ------------------
+       IF ( sun_in_south )  THEN
+          startpos_saa_float = 180.0_wp / 10.0_wp
+       ELSE 
+          startpos_saa_float = saa / 10.0_wp
+       ENDIF
+       radiance_array_temp( 0:35,:) = radiance_array
+       radiance_array_temp(36:71,:) = radiance_array
+       xfactor = (startpos_saa_float) - INT(startpos_saa_float)
+       DO  ai = 0, 35
+          DO  zi = 0, 9
+             radiance_array(ai,zi) = ( radiance_array_temp( 36 - INT( startpos_saa_float ) - 1 + ai , zi ) *          &
+                                     ( xfactor ) )                                                                    &
+                                     + ( radiance_array_temp( 36 - INT( startpos_saa_float ) + ai , zi )              &
+                                     * ( 1.0_wp - xfactor ) )
+          ENDDO
+       ENDDO 
+!        
+!      
+!--    calculate Projectionarea for direct beam -----------------------------'
+       projection_area_direct_temp( 0:35,:) = projection_area
+       projection_area_direct_temp(36:71,:) = projection_area
+       yfactor = ( sza / 10.0_wp ) - INT( sza / 10.0_wp )
+       xfactor = ( startpos_saa_float ) - INT( startpos_saa_float )
+       projection_area_direct_beam = ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)  ) * &
+                                     ( 1.0_wp - xfactor ) * ( 1.0_wp - yfactor ) ) +                                  &
+                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)  ) * &
+                                     (          xfactor ) * ( 1.0_wp - yfactor ) ) +                                  &
+                                     ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)+1) * &
+                                     ( 1.0_wp - xfactor ) * (          yfactor ) ) +                                  &
+                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)+1) * &
+                                     (          xfactor ) * (          yfactor ) )
+!                                                
+!                                                
+!                                                
+       DO  i = nxl, nxr    !nxlg, nxrg
+          DO  j = nys, nyn    !nysg, nyng 
+!                   
+! !--        extract obstruction from IBSET-Integer_Array ------------------'
+             IF (consider_obstructions ) THEN
+                obstruction_temp1 = building_obstruction_f%var_3d(:,j,i)
+                IF (obstruction_temp1(0)  .NE.  9) THEN 
+                   DO  pobi = 0, 44 
+                      DO  bi = 0, 7 
+                         IF ( btest( obstruction_temp1(pobi), bi )  .EQV.  .TRUE.) THEN
+                            obstruction_temp2( ( pobi * 8 ) + bi ) = 1
+                         ELSE
+                            obstruction_temp2( ( pobi * 8 ) + bi ) = 0
+                         ENDIF
+                      ENDDO
+                   ENDDO        
+                   DO  zi = 0, 9                                          
+                      obstruction(:,zi) = obstruction_temp2( zi * 36 :( zi * 36) + 35 )
+                   ENDDO
+                ELSE 
+                   obstruction(:,:) = 0
+                ENDIF
+             ENDIF
+!              
+! !--        calculated human exposure ------------------'  
+             diffuse_exposure = SUM( radiance_array * projection_area * integration_array * obstruction )      
+         
+             obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) ) 
+             IF (sza  .GE.  89.99_wp) THEN
+                sza = 89.99999_wp
+             ENDIF
+!             
+!--          calculate direct normal irradiance (direct beam) ------------------'
+             direct_exposure = ( irradiance(1) / cos( pi * sza / 180.0_wp ) ) * &
+             projection_area_direct_beam * obstruction_direct_beam  
+               
+             vitd3_exposure(j,i) = ( diffuse_exposure + direct_exposure ) / 1000.0_wp * 70.97_wp 
+             ! unit = international units vitamin D per second              
+          ENDDO
+       ENDDO
+    ENDIF
+
+ END SUBROUTINE uvem_calc_exposure
 
  END MODULE biometeorology_mod
