Changeset 4540 for palm/trunk


Ignore:
Timestamp:
May 18, 2020 3:23:29 PM (15 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/biometeorology_mod.f90

    r4535 r4540  
    11!> @file biometeorology_mod.f90
    2 !--------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of PALM-4U.
    44!
    5 ! PALM-4U is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM-4U is free software: You can redistribute it and/or modify it under the terms of the GNU
     6! General Public License as published by the Free Software Foundation, either version 3 of the
     7! License, or (at your option) any later version.
     8!
     9! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
     10! the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 2018-2020 Deutscher Wetterdienst (DWD)
    1817! Copyright 2018-2020 Institute of Computer Science, Academy of Sciences, Prague
    1918! Copyright 2018-2020 Leibniz Universitaet Hannover
    20 !--------------------------------------------------------------------------------!
     19!--------------------------------------------------------------------------------------------------!
    2120!
    2221! Current revisions:
     
    2726! -----------------
    2827! $Id$
     28! file re-formatted to follow the PALM coding standard
     29!
     30! 4535 2020-05-15 12:07:23Z raasch
    2931! bugfix for restart data format query
    30 ! 
     32!
    3133! 4517 2020-05-03 14:29:30Z raasch
    3234! added restart with MPI-IO for reading local arrays
    33 ! 
     35!
    3436! 4495 2020-04-13 20:11:20Z raasch
    3537! restart data handling with MPI-IO added
    36 ! 
     38!
    3739! 4493 2020-04-10 09:49:43Z pavelkrc
    3840! Revise bad formatting
    39 ! 
     41!
    4042! 4286 2019-10-30 16:01:14Z resler
    4143! implement new palm_date_time_mod
    42 ! 
     44!
    4345! 4223 2019-09-10 09:20:47Z gronemeier
    4446! Corrected "Former revisions" section
    45 ! 
     47!
    4648! 4168 2019-08-16 13:50:17Z suehring
    4749! Replace function get_topography_top_index by topo_top_ind
    48 ! 
     50!
    4951! 4144 2019-08-06 09:11:47Z raasch
    5052! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
    51 ! 
     53!
    5254! 4127 2019-07-30 14:47:10Z suehring
    5355! Output for bio_mrt added (merge from branch resler)
    54 ! 
     56!
    5557! 4126 2019-07-30 11:09:11Z gronemeier
    5658! renamed vitd3_exposure_av into vitd3_dose,
    5759! renamed uvem_calc_exposure into bio_calculate_uv_exposure
    58 ! 
     60!
    5961! 3885 2019-04-11 11:29:34Z kanani
    60 ! Changes related to global restructuring of location messages and introduction
    61 ! of additional debug messages
    62 ! 
     62! Changes related to global restructuring of location messages and introduction of additional debug
     63! messages
     64!
    6365! 3753 2019-02-19 14:48:54Z dom_dwd_user
    64 ! - Added automatic setting of mrt_nlevels in case it was not part of
    65 ! radiation_parameters namelist (or set to 0 accidentially).
     66! - Added automatic setting of mrt_nlevels in case it was not part of radiation_parameters namelist
     67!   (or set to 0 accidentially).
    6668! - Minor speed improvoemnts in perceived temperature calculations.
    6769! - Perceived temperature regression arrays now declared as PARAMETERs.
    68 ! 
     70!
    6971! 3750 2019-02-19 07:29:39Z dom_dwd_user
    7072! - Added addittional safety meassures to bio_calculate_thermal_index_maps.
    7173! - Replaced several REAL (un-)equality comparisons.
    72 ! 
     74!
    7375! 3742 2019-02-14 11:25:22Z dom_dwd_user
    74 ! - Allocation of the input _av grids was moved to the "sum" section of
    75 ! bio_3d_data_averaging to make sure averaging is only done once!
    76 ! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to
    77 ! time_integration to make sure averaged input is updated before calculating.
    78 ! 
     76! - Allocation of the input _av grids was moved to the "sum" section of bio_3d_data_averaging to
     77!   make sure averaging is only done once!
     78! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to time_integration to
     79!   make sure averaged input is updated before calculating.
     80!
    7981! 3740 2019-02-13 12:35:12Z dom_dwd_user
    80 ! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after
    81 ! 'bio_<index>' in the output section of the p3d file.
    82 ! 
     82! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after 'bio_<index>' in the
     83!   output section of the p3d file.
     84!
    8385! 3739 2019-02-13 08:05:17Z dom_dwd_user
    84 ! - Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices
    85 ! set as output quantities.
     86! - Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices set as output
     87!   quantities.
    8688! - Renamed flags "bio_<index>" to "do_calculate_<index>" for better readability
    8789! - Removed everything related to "time_bio_results" as this is never used.
    88 ! - Moved humidity warning to check_data_output
     90! - Moved humidity warning to check_data_output.
    8991! - Fixed bug in mrt calculation introduced with my commit yesterday.
    90 ! 
     92!
    9193! 3735 2019-02-12 09:52:40Z dom_dwd_user
    92 ! - Fixed auto-setting of thermal index calculation flags by output
    93 as originally proposed by resler.
     94! - Fixed auto-setting of thermal index calculation flags by output as originally proposed by
     95 resler.
    9496! - removed bio_pet and outher configuration variables.
    9597! - Updated namelist.
    96 ! 
     98!
    9799! 3711 2019-01-31 13:44:26Z knoop
    98100! Introduced interface routine bio_init_checks + small error message changes
    99 ! 
     101!
    100102! 3693 2019-01-23 15:20:53Z dom_dwd_user
    101 ! Added usage of time_averaged mean radiant temperature, together with calculation,
    102 ! grid and restart routines. General cleanup and commenting.
    103 ! 
     103! Added usage of time_averaged mean radiant temperature, together with calculation, grid and restart
     104! routines. General cleanup and commenting.
     105!
    104106! 3685 2019-01-21 01:02:11Z knoop
    105107! Some interface calls moved to module_interface + cleanup
    106 ! 
     108!
    107109! 3650 2019-01-04 13:01:33Z kanani
    108110! Bugfixes and additions for enabling restarts with biometeorology
    109 ! 
     111!
    110112! 3448 2018-10-29 18:14:31Z kanani
    111113! Initial revision
    112114!
    113 ! 
     115!
    114116!
    115117! Authors:
     
    123125! ------------
    124126!> Biometeorology module consisting of two parts:
    125 !> 1.: Human thermal comfort module calculating thermal perception of a sample
    126 !> human being under the current meteorological conditions.
     127!> 1.: Human thermal comfort module calculating thermal perception of a sample human being under the
     128!> current meteorological conditions.
    127129!> 2.: Calculation of vitamin-D weighted UV exposure
    128130!>
    129131!> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
    130132!>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
    131 !> @todo Comments start with capital letter --> "!-- Include..." 
     133!> @todo Comments start with capital letter --> "!-- Include..."
    132134!> @todo uv_vitd3dose-->new output type necessary (cumulative)
    133135!> @todo consider upwelling radiation in UV
     
    136138!>
    137139!> @bug  no known bugs by now
    138 !------------------------------------------------------------------------------!
     140!--------------------------------------------------------------------------------------------------!
    139141 MODULE biometeorology_mod
    140142
    141     USE arrays_3d,                                                             &
     143    USE arrays_3d,                                                                                 &
    142144        ONLY:  pt, p, u, v, w, q
    143145
    144     USE averaging,                                                             &
     146    USE averaging,                                                                                 &
    145147        ONLY:  pt_av, q_av, u_av, v_av, w_av
    146148
    147     USE basic_constants_and_equations_mod,                                     &
    148         ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
    149 
    150     USE control_parameters,                                                    &
    151         ONLY:  average_count_3d, biometeorology,                               &
    152                debug_output,                                                   &
    153                dz, dz_stretch_factor,                                          &
    154                dz_stretch_level, humidity, initializing_actions, nz_do3d,      &
     149    USE basic_constants_and_equations_mod,                                                         &
     150        ONLY: degc_to_k, c_p, l_v, magnus, pi, sigma_sb
     151
     152    USE control_parameters,                                                                        &
     153        ONLY:  average_count_3d, biometeorology,                                                   &
     154               debug_output,                                                                       &
     155               dz, dz_stretch_factor,                                                              &
     156               dz_stretch_level, humidity, initializing_actions, nz_do3d,                          &
    155157               restart_data_format_output, surface_pressure
    156158
    157     USE grid_variables,                                                        &
     159    USE grid_variables,                                                                            &
    158160        ONLY:  ddx, dx, ddy, dy
    159161
    160     USE indices,                                                               &
    161         ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,   &
    162                nysg, nyng, topo_top_ind
     162    USE indices,                                                                                   &
     163        ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg, nysg, nyng,           &
     164               topo_top_ind
    163165
    164166    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
    165167
    166     USE netcdf_data_input_mod,                                                 &
    167         ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,       &
    168                uvem_irradiance_f, uvem_integration_f, building_obstruction_f
    169 
    170     USE palm_date_time_mod,                                                    &
     168    USE netcdf_data_input_mod,                                                                     &
     169        ONLY: building_obstruction_f, netcdf_data_input_uvem, uvem_integration_f,                  &
     170              uvem_irradiance_f, uvem_projarea_f, uvem_radiance_f
     171
     172    USE palm_date_time_mod,                                                                        &
    171173        ONLY:  get_date_time
    172174!
    173175!-- Import radiation model to obtain input for mean radiant temperature
    174     USE radiation_model_mod,                                                   &
    175         ONLY:  ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                    &
    176                mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation,                     &
    177                radiation_interactions, rad_sw_in,                              &
    178                rad_sw_out, rad_lw_in, rad_lw_out
     176    USE radiation_model_mod,                                                                       &
     177        ONLY:  id, ix, iy, iz, mrt_include_sw, mrt_nlevels,                                        &
     178               mrtbl, mrtinlw, mrtinsw, nmrtbl, radiation,                                         &
     179               rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out, radiation_interactions
    179180
    180181    USE restart_data_mpi_io_mod,                                                                   &
     
    184185    IMPLICIT NONE
    185186
    186     PRIVATE
    187 
    188187!
    189188!-- Declare all global variables within the module (alphabetical order)
    190     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (degree_C)
    191     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct      !< PT results   (degree_C)
    192     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci       !< UTCI results (degree_C)
    193     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet        !< PET results  (degree_C)
    194 !
    195 !-- Grids for averaged thermal indices
    196     REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
    197     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_av_grid  !< tmrt results (degree_C)
    198     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (degree_C)
    199     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (degree_C)
    200     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av        !< PET results (aver. input)  (degree_C)
    201 
    202 
    203     INTEGER( iwp ) ::  bio_cell_level     !< cell level biom calculates for
    204     REAL ( wp )    ::  bio_output_height  !< height output is calculated in m
    205     REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
    206     REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
    207     REAL ( wp ), PARAMETER ::  bio_fill_value = -9999._wp  !< set module fill value, replace by global fill value as soon as available
    208 !
    209 !--
    210     LOGICAL ::  thermal_comfort  = .FALSE.  !< Enables or disables the entire thermal comfort part
    211     LOGICAL ::  do_average_theta = .FALSE.  !< switch: do theta averaging in this module? (if .FALSE. this is done globally)
    212     LOGICAL ::  do_average_q     = .FALSE.  !< switch: do e averaging in this module?
    213     LOGICAL ::  do_average_u     = .FALSE.  !< switch: do u averaging in this module?
    214     LOGICAL ::  do_average_v     = .FALSE.  !< switch: do v averaging in this module?
    215     LOGICAL ::  do_average_w     = .FALSE.  !< switch: do w averaging in this module?
    216     LOGICAL ::  do_average_mrt   = .FALSE.  !< switch: do mrt averaging in this module?
     189    INTEGER(iwp) ::  bio_cell_level     !< cell level biom calculates for
     190
     191    LOGICAL ::  thermal_comfort        = .FALSE.  !< Enables or disables the entire thermal comfort part
     192    LOGICAL ::  do_average_theta       = .FALSE.  !< switch: do theta averaging in this module? (if .FALSE. this is done globally)
     193    LOGICAL ::  do_average_q           = .FALSE.  !< switch: do e averaging in this module?
     194    LOGICAL ::  do_average_u           = .FALSE.  !< switch: do u averaging in this module?
     195    LOGICAL ::  do_average_v           = .FALSE.  !< switch: do v averaging in this module?
     196    LOGICAL ::  do_average_w           = .FALSE.  !< switch: do w averaging in this module?
     197    LOGICAL ::  do_average_mrt         = .FALSE.  !< switch: do mrt averaging in this module?
    217198    LOGICAL ::  average_trigger_perct  = .FALSE.  !< update averaged input on call to bio_perct?
    218199    LOGICAL ::  average_trigger_utci   = .FALSE.  !< update averaged input on call to bio_utci?
     
    227208    LOGICAL ::  do_calculate_mrt2d     = .FALSE.  !< Turn index MRT 2D (averaged or inst) on or off
    228209
     210    REAL(wp)    ::  bio_output_height  !< height output is calculated in m
     211
     212    REAL(wp), PARAMETER ::  bio_fill_value = -9999.0_wp  !< set module fill value, replace by global fill value as soon as available
     213    REAL(wp), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
     214    REAL(wp), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
     215
     216    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
     217
     218    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct      !< PT results   (degree_C)
     219    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet        !< PET results  (degree_C)
     220    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (degree_C)
     221    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci       !< UTCI results (degree_C)
     222!
     223!-- Grids for averaged thermal indices
     224    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (degree_C)
     225    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av        !< PET results (aver. input)  (degree_C)
     226    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_av_grid  !< tmrt results (degree_C)
     227    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (degree_C)
     228
    229229!
    230230!-- UVEM parameters from here
    231 !
    232 !-- Declare all global variables within the module (alphabetical order)
    233     INTEGER(iwp) ::  bio_nmrtbl
    234231    INTEGER(iwp) ::  ai                      = 0  !< loop index in azimuth direction
    235232    INTEGER(iwp) ::  bi                      = 0  !< loop index of bit location within an 8bit-integer (one Byte)
     233    INTEGER(iwp) ::  bio_nmrtbl
    236234    INTEGER(iwp) ::  clothing                = 1  !< clothing (0=unclothed, 1=Arms,Hands,Face free, 3=Hand,Face free)
    237235    INTEGER(iwp) ::  iq                      = 0  !< loop index of irradiance quantity
    238236    INTEGER(iwp) ::  pobi                    = 0  !< loop index of the position of corresponding byte within ibset byte vektor
    239     INTEGER(iwp) ::  obstruction_direct_beam = 0  !< Obstruction information for direct beam   
     237    INTEGER(iwp) ::  obstruction_direct_beam = 0  !< Obstruction information for direct beam
    240238    INTEGER(iwp) ::  zi                      = 0  !< loop index in zenith direction
    241239
    242     INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0  !< temporary obstruction information stored with ibset 
     240    INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0  !< temporary obstruction information stored with ibset
    243241    INTEGER(iwp),    DIMENSION(0:359) ::  obstruction_temp2 = 0  !< restored temporary obstruction information from ibset file
    244242
     
    251249
    252250    REAL(wp) ::  diffuse_exposure            =   0.0_wp  !< calculated exposure by diffuse radiation
    253     REAL(wp) ::  direct_exposure             =   0.0_wp  !< calculated exposure by direct solar beam   
    254     REAL(wp) ::  orientation_angle           =   0.0_wp  !< orientation of front/face of the human model   
     251    REAL(wp) ::  direct_exposure             =   0.0_wp  !< calculated exposure by direct solar beam
     252    REAL(wp) ::  orientation_angle           =   0.0_wp  !< orientation of front/face of the human model
    255253    REAL(wp) ::  projection_area_direct_beam =   0.0_wp  !< projection area for direct solar beam
    256254    REAL(wp) ::  saa                         = 180.0_wp  !< solar azimuth angle
     
    261259    REAL(wp) ::  yfactor                     =   0.0_wp  !< relative y-position used for interpolation
    262260
    263     REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp  !< iradiance values extracted from irradiance lookup table 
     261    REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp  !< iradiance values extracted from irradiance lookup table
    264262
    265263    REAL(wp), DIMENSION(0:2,0:90) ::  irradiance_lookup_table      = 0.0_wp  !< irradiance lookup table
     
    269267    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_direct_temp  = 0.0_wp  !< temporary projection area for direct solar beam
    270268    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_temp         = 0.0_wp  !< temporary projection area for all directions
    271     REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp  !< radiance extracted from radiance_lookup_table 
     269    REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp  !< radiance extracted from radiance_lookup_table
    272270    REAL(wp), DIMENSION(0:71,0:9) ::  radiance_array_temp          = 0.0_wp  !< temporary radiance data
    273271
     
    277275    REAL(wp), DIMENSION(0:35,0:9,0:90) ::  radiance_lookup_table   = 0.0_wp  !< radiance lookup table
    278276
     277
     278    PRIVATE
     279
    279280!
    280281!-- INTERFACES that must be available to other modules (alphabetical order)
    281 
    282     PUBLIC bio_3d_data_averaging, bio_check_data_output,                       &
    283     bio_calculate_mrt_grid, bio_calculate_thermal_index_maps, bio_calc_ipt,    &
    284     bio_check_parameters, bio_data_output_3d, bio_data_output_2d,              &
    285     bio_define_netcdf_grid, bio_get_thermal_index_input_ij, bio_header,        &
    286     bio_init, bio_init_checks, bio_parin, thermal_comfort,                     &
    287     bio_nmrtbl, bio_wrd_local, bio_rrd_local, bio_wrd_global, bio_rrd_global
     282    PUBLIC bio_3d_data_averaging,  bio_calculate_mrt_grid, bio_calculate_thermal_index_maps,       &
     283           bio_calc_ipt, bio_check_data_output, bio_check_parameters,                              &
     284           bio_data_output_2d, bio_data_output_3d,  bio_define_netcdf_grid,                        &
     285           bio_get_thermal_index_input_ij, bio_header, bio_init, bio_init_checks, bio_nmrtbl,      &
     286           bio_parin, bio_rrd_global, bio_rrd_local, bio_wrd_global, bio_wrd_local, thermal_comfort
    288287!
    289288!-- UVEM PUBLIC variables and methods
     
    348347    END INTERFACE bio_header
    349348!
    350 !-- Initialization actions 
     349!-- Initialization actions
    351350    INTERFACE bio_init
    352351       MODULE PROCEDURE bio_init
     
    393392
    394393
    395 !------------------------------------------------------------------------------!
     394!--------------------------------------------------------------------------------------------------!
    396395! Description:
    397396! ------------
    398 !> Sum up and time-average biom input quantities as well as allocate
    399 !> the array necessary for storing the average.
    400 !> There is a considerable difference to the 3d_data_averaging subroutines
    401 !> used by other modules:
    402 !> For the thermal indices, the module needs to average the input conditions
    403 !> not the result!
    404 !------------------------------------------------------------------------------!
     397!> Sum up and time-average biom input quantities as well as allocate the array necessary for storing
     398!> the average.
     399!> There is a considerable difference to the 3d_data_averaging subroutines used by other modules:
     400!> For the thermal indices, the module needs to average the input conditions, not the result!
     401!--------------------------------------------------------------------------------------------------!
    405402 SUBROUTINE bio_3d_data_averaging( mode, variable )
    406403
    407404    IMPLICIT NONE
    408405
    409     CHARACTER (LEN=*) ::  mode     !< averaging mode: allocate, sum, or average
     406    CHARACTER (LEN=*) ::  mode     !< Averaging mode: allocate, sum, or average
    410407    CHARACTER (LEN=*) ::  variable !< The variable in question
    411408
     
    425422                ENDIF
    426423                mrt_av_grid = 0.0_wp
    427                 do_average_mrt = .FALSE.  !< overwrite if that was enabled somehow
     424                do_average_mrt = .FALSE.  !< Overwrite if that was enabled somehow
    428425
    429426
     
    431428
    432429!
    433 !--          Averaging, as well as the allocation of the required grids must be
    434 !--          done only once, independent from for how many thermal indices
    435 !--          averaged output is desired.
    436 !--          Therefore wee need to memorize which index is the one that controls
    437 !--          the averaging (what must be the first thermal index called).
    438 !--          Indices are in unknown order as depending on the input file,
    439 !--          determine first index to average und update only once
    440 !
    441 !--          Only proceed here if this was not done for any index before. This
    442 !--          is done only once during the whole model run.
    443              IF ( .NOT. average_trigger_perct  .AND.                           &
    444                   .NOT. average_trigger_utci   .AND.                           &
    445                   .NOT. average_trigger_pet    .AND.                           &
     430!--          Averaging, as well as the allocation of the required grids must be done only once,
     431!--          independent of for how many thermal indices averaged output is desired.
     432!--          Therefore we need to memorize which index is the one that controls the averaging
     433!--          (what must be the first thermal index called).
     434!--          Indices are in unknown order as depending on the input file, determine first index to
     435!--          average und update only once.
     436!
     437!--          Only proceed here if this was not done for any index before. This is done only once
     438!--          during the whole model run.
     439             IF ( .NOT. average_trigger_perct  .AND.                                               &
     440                  .NOT. average_trigger_utci   .AND.                                               &
     441                  .NOT. average_trigger_pet    .AND.                                               &
    446442                  .NOT. average_trigger_mrt )  THEN
    447443!
     
    461457             ENDIF
    462458!
    463 !--          Allocation of the input _av grids was moved to the "sum" section to
    464 !--          make sure averaging is only done once!
    465 
     459!--          Allocation of the input _av grids was moved to the "sum" section to make sure averaging
     460!--          is only done once!
    466461
    467462          CASE ( 'uvem_vitd3dose*' )
     463
    468464             IF ( .NOT. ALLOCATED( vitd3_dose ) )  THEN
    469465                ALLOCATE( vitd3_dose(nysg:nyng,nxlg:nxrg) )
     
    482478          CASE ( 'bio_mrt' )
    483479!
    484 !--          Consider the case 'bio_mrt' is called after some thermal index. In
    485 !--          that case do_average_mrt will be .TRUE. leading to a double-
    486 !--          averaging.
     480!--          Consider the case 'bio_mrt' is called after some thermal index. In that case
     481!            do_average_mrt will be .TRUE. leading to a double-averaging.
    487482             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
    488483
    489484                IF ( mrt_include_sw )  THEN
    490                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
    491                       ( ( human_absorb * mrtinsw(:) +                          &
    492                       mrtinlw(:) ) /                                           &
    493                       ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
     485                   mrt_av_grid(:) = mrt_av_grid(:)  +                                              &
     486                                    ( ( human_absorb * mrtinsw(:) +                                &
     487                                    mrtinlw(:) ) / ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
    494488                ELSE
    495                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
    496                       ( mrtinlw(:) /                                           &
    497                       ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
     489                   mrt_av_grid(:) = mrt_av_grid(:) +                                               &
     490                                    ( mrtinlw(:) / ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
    498491                ENDIF
    499492             ENDIF
     
    501494          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
    502495!
    503 !--          Only continue if the current index is the one to trigger the input
    504 !--          averaging, see above
    505              IF ( average_trigger_perct  .AND.  TRIM( variable ) /=            &
    506                 'bio_perct*')  RETURN
    507              IF ( average_trigger_utci   .AND.  TRIM( variable ) /=            &
    508                 'bio_utci*')   RETURN
    509              IF ( average_trigger_pet    .AND.  TRIM( variable ) /=            &
    510                 'bio_pet*')    RETURN
    511              IF ( average_trigger_mrt    .AND.  TRIM( variable ) /=            &
    512                 'bio_mrt*')    RETURN
    513 !
    514 !--          Now memorize which of the input grids are not averaged by other
    515 !--          modules. Set averaging switch to .TRUE. and allocate the respective
    516 !--          grid in that case.
     496!--          Only continue if the current index is the one to trigger the input averaging, see
     497!--          above.
     498             IF ( average_trigger_perct  .AND.  TRIM( variable ) /= 'bio_perct*')    RETURN
     499             IF ( average_trigger_utci   .AND.  TRIM( variable ) /= 'bio_utci*' )    RETURN
     500             IF ( average_trigger_pet    .AND.  TRIM( variable ) /= 'bio_pet*'  )    RETURN
     501             IF ( average_trigger_mrt    .AND.  TRIM( variable ) /= 'bio_mrt*'  )    RETURN
     502!
     503!--          Now memorize which of the input grids are not averaged by other modules. Set averaging
     504!--          switch to .TRUE. and allocate the respective grid in that case.
    517505             IF ( .NOT. ALLOCATED( pt_av ) )  THEN  !< if not averaged by other module
    518506                ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    546534
    547535!
    548 !-- u_av, v_av and w_av are always allocated
     536!--          u_av, v_av and w_av are always allocated
    549537             IF ( .NOT. ALLOCATED( u_av ) )  THEN
    550538                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    600588
    601589                IF ( mrt_include_sw )  THEN
    602                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
    603                       ( ( human_absorb * mrtinsw(:) +                          &
    604                       mrtinlw(:) ) /                                           &
    605                       ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
    606                 ELSE
    607                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
    608                       ( mrtinlw(:) /                                           &
    609                       ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
     590                   mrt_av_grid(:) = mrt_av_grid(:) +                                               &
     591                                    ( ( human_absorb * mrtinsw(:) +                                &
     592                                    mrtinlw(:) ) /                                                 &
     593                                      ( human_emiss  * sigma_sb ) )**0.25_wp - degc_to_k
     594               ELSE
     595                   mrt_av_grid(:) = mrt_av_grid(:) +                                               &
     596                                    ( mrtinlw(:) /                                                 &
     597                                    ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k
    610598                ENDIF
    611599             ENDIF
     
    632620          CASE ( 'bio_mrt' )
    633621!
    634 !--          Consider the case 'bio_mrt' is called after some thermal index. In
    635 !--          that case do_average_mrt will be .TRUE. leading to a double-
    636 !--          averaging.
     622!--          Consider the case 'bio_mrt' is called after some thermal index. In that case
     623!--          do_average_mrt will be .TRUE. leading to a double-averaging.
    637624             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
    638625                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
     
    642629!
    643630!--          Only continue if update index, see above
    644              IF ( average_trigger_perct  .AND.                                 &
     631             IF ( average_trigger_perct  .AND.                                                     &
    645632                TRIM( variable ) /= 'bio_perct*' )  RETURN
    646              IF ( average_trigger_utci  .AND.                                  &
     633             IF ( average_trigger_utci  .AND.                                                      &
    647634                TRIM( variable ) /= 'bio_utci*' )  RETURN
    648              IF ( average_trigger_pet   .AND.                                  &
     635             IF ( average_trigger_pet   .AND.                                                      &
    649636                TRIM( variable ) /= 'bio_pet*' )  RETURN
    650              IF ( average_trigger_mrt   .AND.                                  &
     637             IF ( average_trigger_mrt   .AND.                                                      &
    651638                TRIM( variable ) /= 'bio_mrt*' )  RETURN
    652639
     
    655642                   DO  j = nys, nyn
    656643                      DO  k = nzb, nzt+1
    657                          pt_av(k,j,i) = pt_av(k,j,i) /                         &
    658                             REAL( average_count_3d, KIND=wp )
     644                         pt_av(k,j,i) = pt_av(k,j,i) /  REAL( average_count_3d, KIND = wp )
    659645                      ENDDO
    660646                   ENDDO
     
    666652                   DO  j = nys, nyn
    667653                      DO  k = nzb, nzt+1
    668                          q_av(k,j,i) = q_av(k,j,i) /                           &
    669                             REAL( average_count_3d, KIND=wp )
     654                         q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    670655                      ENDDO
    671656                   ENDDO
     
    677662                   DO  j = nysg, nyng
    678663                      DO  k = nzb, nzt+1
    679                          u_av(k,j,i) = u_av(k,j,i) /                           &
    680                             REAL( average_count_3d, KIND=wp )
     664                         u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    681665                      ENDDO
    682666                   ENDDO
     
    688672                   DO  j = nysg, nyng
    689673                      DO  k = nzb, nzt+1
    690                          v_av(k,j,i) = v_av(k,j,i) /                           &
    691                             REAL( average_count_3d, KIND=wp )
     674                         v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    692675                      ENDDO
    693676                   ENDDO
     
    699682                   DO  j = nysg, nyng
    700683                      DO  k = nzb, nzt+1
    701                          w_av(k,j,i) = w_av(k,j,i) /                           &
    702                             REAL( average_count_3d, KIND=wp )
     684                         w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    703685                      ENDDO
    704686                   ENDDO
     
    707689
    708690             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
    709                 mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d,      &
    710                 KIND=wp )
     691                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND = wp )
    711692             ENDIF
    712693
    713694!
    714 !--     No averaging for UVEM since we are calculating a dose (only sum is
    715 !--     calculated and saved to av.nc file)
    716 
     695!--     No averaging for UVEM SINce we are calculating a dose (only sum is calculated and saved to
     696!--     av.nc file)
    717697        END SELECT
    718698
     
    724704
    725705
    726 !------------------------------------------------------------------------------!
     706!--------------------------------------------------------------------------------------------------!
    727707! Description:
    728708! ------------
    729709!> Check data output for biometeorology model
    730 !------------------------------------------------------------------------------!
     710!--------------------------------------------------------------------------------------------------!
    731711 SUBROUTINE bio_check_data_output( var, unit, i, j, ilen, k )
    732712
    733     USE control_parameters,                                                    &
     713    USE control_parameters,                                                                        &
    734714        ONLY: data_output, message_string
    735715
     
    747727!
    748728!--    Allocate a temporary array with the desired output dimensions.
    749 !--    Arrays for time-averaged thermal indices are also allocated here because
    750 !--    they are not running through the standard averaging procedure in
    751 !--    bio_3d_data_averaging as the values of the averaged thermal indices are
    752 !--    derived in a single step based on priorly averaged arrays (see
     729!--    Arrays for time-averaged thermal indices are also allocated here because they are not running
     730!--    through the standard averaging procedure in bio_3d_data_averaging as the values of the
     731!--    averaged thermal indices are derived in a SINgle step based on priorly averaged arrays (see
    753732!--    bio_calculate_thermal_index_maps).
    754733       CASE ( 'bio_mrt', 'bio_mrt*' )
     
    772751                perct = REAL( bio_fill_value, KIND = wp )
    773752             ENDIF
    774           ELSE                              !< if averaged input
     753          ELSE                               !< if averaged input
    775754             do_calculate_perct_av = .TRUE.
    776755             IF ( .NOT. ALLOCATED( perct_av ) )  THEN
     
    816795
    817796       CASE ( 'uvem_vitd3*' )
    818 !           IF (  .NOT. uv_exposure )  THEN
    819 !              message_string = 'output of "' // TRIM( var ) // '" requi' //     &
     797!           IF ( .NOT. uv_exposure )  THEN
     798!              message_string = 'output of "' // TRIM( var ) // '" requi' //                       &
    820799!                       'res a namelist &uvexposure_par'
    821800!              CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
    822801!           ENDIF
    823802          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    824              message_string = 'illegal value for data_output: "' //            &
    825                               TRIM( var ) // '" & only 2d-horizontal ' //      &
     803             message_string = 'illegal value for data_output: "' //                                &
     804                              TRIM( var ) // '" & only 2d-horizontal ' //                          &
    826805                              'cross sections are allowed for this value'
    827806             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
     
    840819!           ENDIF
    841820          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    842              message_string = 'illegal value for data_output: "' //            &
    843                               TRIM( var ) // '" & only 2d-horizontal ' //      &
     821             message_string = 'illegal value for data_output: "' //                                &
     822                              TRIM( var ) // '" & only 2d-horizontal ' //                          &
    844823                              'cross sections are allowed for this value'
    845824             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
     
    857836
    858837!
    859 !--    Further checks if thermal comfort output is desired.
     838!-- Further checks if thermal comfort output is desired.
    860839    IF ( thermal_comfort  .AND.  unit == 'degree_C' )  THEN
    861840!
    862 !--    Break if required modules "radiation" is not avalable.
     841!--    Break if required modules "radiation" is not available.
    863842       IF ( .NOT.  radiation )  THEN
    864           message_string = 'output of "' // TRIM( var ) // '" require'         &
    865                            // 's radiation = .TRUE.'
     843          message_string = 'output of "' // TRIM( var ) // '" require' // 's radiation = .TRUE.'
    866844          CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
    867845          unit = 'illegal'
    868846       ENDIF
    869847!
    870 !--    All "thermal_comfort" outputs except from 'bio_mrt' will also need
    871 !--    humidity input. Check also for that.
     848!--    All "thermal_comfort" outputs except from 'bio_mrt' will also need  humidity input. Check
     849!--    also for that.
    872850       IF ( TRIM( var ) /= 'bio_mrt' )  THEN
    873851          IF ( .NOT.  humidity )  THEN
    874              message_string = 'The estimation of thermal comfort '    //       &
    875                               'requires air humidity information, but ' //     &
     852             message_string = 'The estimation of thermal comfort '    //                           &
     853                              'requires air humidity information, but ' //                         &
    876854                              'humidity module is disabled!'
    877855             CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
     
    885863 END SUBROUTINE bio_check_data_output
    886864
    887 !------------------------------------------------------------------------------!
     865!--------------------------------------------------------------------------------------------------!
    888866! Description:
    889867! ------------
    890868!> Check parameters routine for biom module
    891869!> Currently unused but might come in handy for future checks?
    892 !------------------------------------------------------------------------------!
     870!--------------------------------------------------------------------------------------------------!
    893871 SUBROUTINE bio_check_parameters
    894872
     
    897875
    898876
    899 
    900877 END SUBROUTINE bio_check_parameters
    901878
    902879
    903 !------------------------------------------------------------------------------!
     880!--------------------------------------------------------------------------------------------------!
    904881! Description:
    905882! ------------
    906883!> Subroutine defining 2D output variables
    907884!> data_output_2d 1188ff
    908 !------------------------------------------------------------------------------!
    909  SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf,           &
    910                                 two_d, nzb_do, nzt_do )
     885!--------------------------------------------------------------------------------------------------!
     886 SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf, two_d, nzb_do, nzt_do)
    911887
    912888
     
    924900!-- Output variables
    925901    CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
     902
    926903    LOGICAL, INTENT(OUT)           ::  found  !< Output found?
    927904    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables,
     
    950927              j = mrtbl(iy,l)
    951928              k = mrtbl(iz,l)
    952               IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.           &
     929              IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.                               &
    953930                 j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
    954931              IF ( av == 0 )  THEN
    955932                 IF ( mrt_include_sw )  THEN
    956                     local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +          &
    957                                     mrtinlw(l) ) /                             &
    958                                     ( human_emiss * sigma_sb ) )**.25_wp -     &
    959                                     degc_to_k
     933                    local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +                              &
     934                                      mrtinlw(l) ) /                                               &
     935                                      ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k
    960936                 ELSE
    961                     local_pf(i,j,k) = ( mrtinlw(l)  /                          &
    962                                     ( human_emiss * sigma_sb ) )**.25_wp -     &
    963                                     degc_to_k
     937                    local_pf(i,j,k) = ( mrtinlw(l) /                                               &
     938                                      ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k
    964939                 ENDIF
    965940              ELSE
     
    10401015
    10411016!
    1042 !--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
    1043 !--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
     1017!--    Before data is transfered to local_pf, transfer is in 2D dummy variable and exchange ghost
     1018!--    points therein. However, at this point this is only required for instantaneous arrays,
     1019!--    time-averaged quantities are already exchanged.
    10441020       CASE ( 'uvem_vitd3*_xy' )        ! 2d-array
    10451021          IF ( av == 0 )  THEN
     
    10771053
    10781054
    1079 !------------------------------------------------------------------------------!
     1055!--------------------------------------------------------------------------------------------------!
    10801056! Description:
    10811057! ------------
    10821058!> Subroutine defining 3D output variables (dummy, always 2d!)
    10831059!> data_output_3d 709ff
    1084 !------------------------------------------------------------------------------!
     1060!--------------------------------------------------------------------------------------------------!
    10851061 SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
    10861062
     
    10941070!-- Input variables
    10951071    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
     1072
    10961073    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
    10971074    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
     
    11001077!-- Output variables
    11011078    LOGICAL, INTENT(OUT) ::  found   !< Output found?
     1079
    11021080    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
    11031081!
     
    11081086    INTEGER(iwp) ::  k    !< Running index, z-dir
    11091087
    1110 !     REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
     1088!    REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
    11111089
    11121090    found = .TRUE.
     
    11201098               j = mrtbl(iy,l)
    11211099               k = mrtbl(iz,l)
    1122                IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.          &
     1100               IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.                              &
    11231101                  j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
    11241102               IF ( av == 0 )  THEN
    11251103                  IF ( mrt_include_sw )  THEN
    1126                      local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +   &
    1127                                     mrtinlw(l) ) /                             &
    1128                                     ( human_emiss * sigma_sb ) )**.25_wp -     &
    1129                                     degc_to_k, KIND = sp )
     1104                     local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +                       &
     1105                                       mrtinlw(l) ) /                                              &
     1106                                       ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k,          &
     1107                                        KIND = sp )
    11301108                  ELSE
    1131                      local_pf(i,j,k) = REAL( ( mrtinlw(l) /                    &
    1132                                     ( human_emiss * sigma_sb ) )**.25_wp -     &
    1133                                     degc_to_k, KIND = sp )
     1109                     local_pf(i,j,k) = REAL( ( mrtinlw(l) /                                        &
     1110                                       ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k,          &
     1111                                        KIND = sp )
    11341112                  ENDIF
    11351113               ELSE
     
    11451123 END SUBROUTINE bio_data_output_3d
    11461124
    1147 !------------------------------------------------------------------------------!
     1125!--------------------------------------------------------------------------------------------------!
    11481126! Description:
    11491127! ------------
     
    11511129!> It is called out from subroutine netcdf_interface_mod.
    11521130!> netcdf_interface_mod 918ff
    1153 !------------------------------------------------------------------------------!
     1131!--------------------------------------------------------------------------------------------------!
    11541132 SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
    11551133
     
    11671145!
    11681146!-- Local variables
     1147    INTEGER(iwp) :: l     !< Length of the var array
     1148
    11691149    LOGICAL      :: is2d  !< Var is 2d?
    1170 
    1171     INTEGER(iwp) :: l     !< Length of the var array
    1172 
    11731150
    11741151    found  = .FALSE.
     
    11961173 END SUBROUTINE bio_define_netcdf_grid
    11971174
    1198 !------------------------------------------------------------------------------!
     1175!--------------------------------------------------------------------------------------------------!
    11991176! Description:
    12001177! ------------
    12011178!> Header output for biom module
    12021179!> header 982
    1203 !------------------------------------------------------------------------------!
     1180!--------------------------------------------------------------------------------------------------!
    12041181 SUBROUTINE bio_header( io )
    12051182
     
    12191196    WRITE( io, 3 )  TRIM( ACHAR( bio_cell_level ) )
    12201197
    1221 1   FORMAT (//' Human thermal comfort module information:'/                    &
     11981   FORMAT (//' Human thermal comfort module information:'/                                        &
    12221199              ' ------------------------------'/)
    122312002   FORMAT ('    --> All indices calculated for a height of (m): ', A )
     
    12271204
    12281205
    1229 !------------------------------------------------------------------------------!
     1206!--------------------------------------------------------------------------------------------------!
    12301207! Description:
    12311208! ------------
    12321209!> Initialization of the HTCM
    12331210!> init_3d_model 1987ff
    1234 !------------------------------------------------------------------------------!
     1211!--------------------------------------------------------------------------------------------------!
    12351212 SUBROUTINE bio_init
    12361213
    1237     USE netcdf_data_input_mod,                                                 &
     1214    USE netcdf_data_input_mod,                                                                     &
    12381215        ONLY:  netcdf_data_input_uvem
    12391216
     
    12451222    IF ( debug_output )  CALL debug_message( 'bio_init', 'start' )
    12461223!
    1247 !-- Determine cell level corresponding to 1.1 m above ground level
    1248 !   (gravimetric center of sample human)
     1224!-- Determine cell level corresponding to 1.1 m above ground level (gravimetric center of sample
     1225!-- human)
    12491226
    12501227    bio_cell_level = 0_iwp
     
    12681245
    12691246
    1270 !------------------------------------------------------------------------------!
     1247!--------------------------------------------------------------------------------------------------!
    12711248! Description:
    12721249! ------------
    12731250!> Checks done after the Initialization
    1274 !------------------------------------------------------------------------------!
     1251!--------------------------------------------------------------------------------------------------!
    12751252 SUBROUTINE bio_init_checks
    12761253
    1277     USE control_parameters,                                                    &
     1254    USE control_parameters,                                                                        &
    12781255        ONLY: message_string
    12791256
    12801257    IF ( (.NOT. radiation_interactions) .AND. ( thermal_comfort ) )  THEN
    1281        message_string = 'The mrt calculation requires ' //                     &
    1282                         'enabled radiation_interactions but it ' //            &
     1258       message_string = 'The mrt calculation requires ' //                                         &
     1259                        'enabled radiation_interactions but it ' //                                &
    12831260                        'is disabled!'
    12841261       CALL message( 'bio_init_checks', 'PAHU03', 1, 2, 0, 6, 0 )
     
    12891266
    12901267
    1291 !------------------------------------------------------------------------------!
     1268!--------------------------------------------------------------------------------------------------!
    12921269! Description:
    12931270! ------------
    12941271!> Parin for &biometeorology_parameters for reading biomet parameters
    1295 !------------------------------------------------------------------------------!
     1272!--------------------------------------------------------------------------------------------------!
    12961273 SUBROUTINE bio_parin
    12971274
     
    13001277!
    13011278!-- Internal variables
    1302     CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file 
    1303 
    1304     NAMELIST /biometeorology_parameters/  thermal_comfort,                     &
    1305                                           clothing,                            &
    1306                                           consider_obstructions,               &
    1307                                           orientation_angle,                   &
    1308                                           sun_in_south,                        &
    1309                                           turn_to_sun,                         &
     1279    CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file
     1280
     1281    NAMELIST /biometeorology_parameters/  clothing,                                                &
     1282                                          consider_obstructions,                                   &
     1283                                          orientation_angle,                                       &
     1284                                          sun_in_south,                                            &
     1285                                          thermal_comfort,                                         &
     1286                                          turn_to_sun,                                             &
    13101287                                          uv_exposure
    13111288
     
    13421319 END SUBROUTINE bio_parin
    13431320
    1344 !------------------------------------------------------------------------------!
     1321!--------------------------------------------------------------------------------------------------!
    13451322! Description:
    13461323! ------------
    13471324!> Read module-specific global restart data (Fortran binary format).
    1348 !------------------------------------------------------------------------------!
     1325!--------------------------------------------------------------------------------------------------!
    13491326 SUBROUTINE bio_rrd_global_ftn( found )
    13501327
    1351     USE control_parameters,                                                    &
     1328    USE control_parameters,                                                                        &
    13521329        ONLY:  length, restart_string
    13531330
     
    13631340
    13641341!
    1365 !--    read control flags to determine if input grids need to be averaged
     1342!--    Read control flags to determine if input grids need to be averaged.
    13661343       CASE ( 'do_average_theta' )
    13671344          READ ( 13 )  do_average_theta
     
    13831360
    13841361!
    1385 !--    read control flags to determine which thermal index needs to trigger averaging
     1362!--    Read control flags to determine which thermal index needs to trigger averaging.
    13861363       CASE ( 'average_trigger_perct' )
    13871364          READ ( 13 )  average_trigger_perct
     
    14071384
    14081385
    1409 !------------------------------------------------------------------------------!
     1386!--------------------------------------------------------------------------------------------------!
    14101387! Description:
    14111388! ------------
    14121389!> Read module-specific global restart data (MPI-IO).
    1413 !------------------------------------------------------------------------------!
     1390!--------------------------------------------------------------------------------------------------!
    14141391 SUBROUTINE bio_rrd_global_mpi
    14151392
     
    14331410
    14341411
    1435 !------------------------------------------------------------------------------!
     1412!--------------------------------------------------------------------------------------------------!
    14361413! Description:
    14371414! ------------
    14381415!> Read module-specific local restart data arrays (Fortran binary format).
    1439 !------------------------------------------------------------------------------!
     1416!--------------------------------------------------------------------------------------------------!
    14401417 SUBROUTINE bio_rrd_local_ftn( found )
    14411418
    14421419
    1443     USE control_parameters,                                                    &
     1420    USE control_parameters,                                                                        &
    14441421        ONLY:  length, restart_string
    14451422
     
    14481425
    14491426
    1450     LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
     1427    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .TRUE., no = .FALSE.
    14511428
    14521429    found = .TRUE.
     
    14871464
    14881465
    1489 !------------------------------------------------------------------------------!
     1466!--------------------------------------------------------------------------------------------------!
    14901467! Description:
    14911468! ------------
    14921469!> Write global restart data for the biometeorology module.
    1493 !------------------------------------------------------------------------------!
     1470!--------------------------------------------------------------------------------------------------!
    14941471 SUBROUTINE bio_wrd_global
    14951472
     
    15191496       WRITE ( 14 )  average_trigger_mrt
    15201497
    1521     ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
     1498    ELSEIF ( TRIM( restart_data_format_output(1:3) ) == 'mpi' )  THEN
    15221499
    15231500       CALL wrd_mpi_io( 'do_average_theta', do_average_theta )
     
    15311508       CALL wrd_mpi_io( 'average_trigger_pet', average_trigger_pet )
    15321509       CALL wrd_mpi_io( 'average_trigger_mrt', average_trigger_mrt )
     1510
    15331511    ENDIF
    15341512
     
    15361514
    15371515
    1538 !------------------------------------------------------------------------------!
     1516!--------------------------------------------------------------------------------------------------!
    15391517! Description:
    15401518! ------------
    15411519!> Write local restart data for the biometeorology module.
    1542 !------------------------------------------------------------------------------!
     1520!--------------------------------------------------------------------------------------------------!
    15431521 SUBROUTINE bio_wrd_local
    15441522
     
    15581536       ENDIF
    15591537
    1560     ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
     1538    ELSEIF ( TRIM( restart_data_format_output(1:3) ) == 'mpi' )  THEN
    15611539
    15621540!
     
    15731551 END SUBROUTINE bio_wrd_local
    15741552
    1575 !------------------------------------------------------------------------------!
     1553!--------------------------------------------------------------------------------------------------!
    15761554! Description:
    15771555! ------------
    15781556!> Calculate biometeorology MRT for all 2D grid
    1579 !------------------------------------------------------------------------------!
     1557!--------------------------------------------------------------------------------------------------!
    15801558 SUBROUTINE bio_calculate_mrt_grid ( av )
    15811559
     
    16051583       IF ( ALLOCATED( mrt_av_grid ) )  THEN
    16061584!
    1607 !-- Iterate over the radiation grid (radiation coordinates) and fill the
    1608 !-- tmrt_av_grid (x, y coordinates) where appropriate:
    1609 !-- tmrt_av_grid is written for all i / j if level (k) matches output height.
     1585!--       Iterate over the radiation grid (radiation coordinates) and fill the tmrt_av_grid
     1586!--       (x, y coordinates) where appropriate: tmrt_av_grid is written for all i / j if level (k)
     1587!--       matches output height.
    16101588          DO  l = 1, nmrtbl
    16111589             i = mrtbl(ix,l)
     
    16141592             IF ( k - topo_top_ind(j,i,0) == bio_cell_level + 1_iwp)  THEN
    16151593!
    1616 !-- Averaging was done before, so we can just copy the result here
     1594!--             Averaging was done before, so we can just copy the result here.
    16171595                tmrt_av_grid(j,i) = mrt_av_grid(l)
    16181596
     
    16251603    ELSE
    16261604!
    1627 !-- Calculate biometeorology MRT from local radiation fluxes calculated by RTM and assign
    1628 !-- into 2D grid. Depending on selected output quantities, tmrt_grid might not have been
    1629 !-- allocated in bio_check_data_output yet.
     1605!--    Calculate biometeorology MRT from local radiation fluxes calculated by RTM and assign into 2D
     1606!--    grid. Depending on selected output quantities, tmrt_grid might not have been allocated in
     1607!--   bio_check_data_output yet.
    16301608       IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
    16311609          ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
     
    16391617          IF ( k - topo_top_ind(j,i,0) == bio_cell_level + 1_iwp)  THEN
    16401618             IF ( mrt_include_sw )  THEN
    1641                  tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +              &
    1642                                   mrtinlw(l) )  /                              &
    1643                                   ( human_emiss * sigma_sb ) )**.25_wp -       &
    1644                                   degc_to_k
     1619                tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +                                   &
     1620                                 mrtinlw(l) )  /                                                   &
     1621                                 ( human_emiss * sigma_sb ) )**0.25_wp -                           &
     1622                                 degc_to_k
    16451623             ELSE
    1646                  tmrt_grid(j,i) = ( mrtinlw(l)  /                              &
    1647                                   ( human_emiss * sigma_sb ) )**.25_wp -       &
    1648                                   degc_to_k
     1624                tmrt_grid(j,i) = ( mrtinlw(l)  /                                                   &
     1625                                 ( human_emiss * sigma_sb ) )**0.25_wp -                           &
     1626                                 degc_to_k
    16491627             ENDIF
    16501628          ENDIF
     
    16521630    ENDIF
    16531631
    1654 END SUBROUTINE bio_calculate_mrt_grid
    1655 
    1656 
    1657 !------------------------------------------------------------------------------!
     1632 END SUBROUTINE bio_calculate_mrt_grid
     1633
     1634
     1635!--------------------------------------------------------------------------------------------------!
    16581636! Description:
    16591637! ------------
    16601638!> Calculate static thermal indices for 2D grid point i, j
    1661 !------------------------------------------------------------------------------!
    1662  SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws,   &
    1663                                             pair, tmrt )
     1639!--------------------------------------------------------------------------------------------------!
     1640 SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws, pair, tmrt )
    16641641
    16651642    IMPLICIT NONE
     
    16711648!
    16721649!-- Output parameters
     1650    REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
     1651    REAL(wp), INTENT ( OUT )    ::  ta    !< Air temperature                 (degree_C)
    16731652    REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (degree_C)
    1674     REAL(wp), INTENT ( OUT )    ::  ta    !< Air temperature                 (degree_C)
    16751653    REAL(wp), INTENT ( OUT )    ::  vp    !< Vapour pressure                 (hPa)
    16761654    REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
    1677     REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
    16781655!
    16791656!-- Internal variables
     
    16841661
    16851662!
    1686 !-- Determine cell level closest to 1.1m above ground
    1687 !   by making use of truncation due to int cast
     1663!-- Determine cell level closest to 1.1m above ground by making use of truncation due to int cast.
    16881664    k = INT( topo_top_ind(j,i,0) + bio_cell_level )  !< Vertical cell center closest to 1.1m
    16891665
     
    17011677       ta = bio_fill_value
    17021678       IF ( ALLOCATED( pt_av ) )  THEN
    1703           ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k
     1679          ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + 0.5_wp ) ) - degc_to_k
    17041680       ENDIF
    17051681
     
    17101686
    17111687       ws = bio_fill_value
    1712        IF ( ALLOCATED( u_av )  .AND.  ALLOCATED( v_av )  .AND.                 &
     1688       IF ( ALLOCATED( u_av )  .AND.  ALLOCATED( v_av )  .AND.                                     &
    17131689          ALLOCATED( w_av ) )  THEN
    1714              ws = ( 0.5_wp * ABS( u_av(k_wind,j,i) + u_av(k_wind,j,i+1) )  +   &
    1715              0.5_wp * ABS( v_av(k_wind,j,i) + v_av(k_wind,j+1,i) )  +          &
    1716              0.5_wp * ABS( w_av(k_wind,j,i) + w_av(k_wind+1,j,i) ) )
     1690             ws = ( 0.5_wp * ABS( u_av(k_wind,j,i) + u_av(k_wind,j,i+1) ) +                        &
     1691                    0.5_wp * ABS( v_av(k_wind,j,i) + v_av(k_wind,j+1,i) ) +                        &
     1692                    0.5_wp * ABS( w_av(k_wind,j,i) + w_av(k_wind+1,j,i) ) )
    17171693       ENDIF
    17181694    ELSE
    17191695!
    1720 !-- Calculate ta from Tp assuming dry adiabatic laps rate
    1721        ta = pt(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
     1696!--    Calculate ta from Tp assuming dry adiabatic laps rate
     1697       ta = pt(k,j,i) - ( 0.0098_wp * dz(1) * (  k + 0.5_wp ) ) - degc_to_k
    17221698
    17231699       vp = bio_fill_value
     
    17261702       ENDIF
    17271703
    1728        ws = ( 0.5_wp * ABS( u(k_wind,j,i) + u(k_wind,j,i+1) )  +               &
    1729           0.5_wp * ABS( v(k_wind,j,i) + v(k_wind,j+1,i) )  +                   &
    1730           0.5_wp * ABS( w(k_wind,j,i) + w(k_wind+1,j,i) ) )
     1704       ws = ( 0.5_wp * ABS( u(k_wind,j,i) + u(k_wind,j,i+1) )  +                                   &
     1705              0.5_wp * ABS( v(k_wind,j,i) + v(k_wind,j+1,i) )  +                                   &
     1706              0.5_wp * ABS( w(k_wind,j,i) + w(k_wind+1,j,i) ) )
    17311707
    17321708    ENDIF
     
    17351711    pair = surface_pressure
    17361712!
    1737 !-- Calculate water vapour pressure at saturation and convert to hPa
    1738 !-- The magnus formula is limited to temperatures up to 333.15 K to
    1739 !   avoid negative values of vp_sat
    1740     IF ( vp > -998._wp )  THEN
     1713!-- Calculate water vapour pressure at saturation and convert to hPa.
     1714!-- The magnus formula is limited to temperatures up to 333.15 K to avoid negative values of vp_sat.
     1715    IF ( vp > -998.0_wp )  THEN
    17411716       vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
    17421717       vp  = vp * pair / ( vp + 0.622_wp )
     
    17441719    ENDIF
    17451720!
    1746 !-- local mtr value at [i,j]
     1721!-- Local mtr value at [i,j]
    17471722    tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
    17481723    IF ( .NOT. average_input )  THEN
     
    17571732
    17581733
    1759 !------------------------------------------------------------------------------!
     1734!--------------------------------------------------------------------------------------------------!
    17601735! Description:
    17611736! ------------
    1762 !> Calculate static thermal indices for any point within a 2D grid
    1763 !> time_integration.f90: 1065ff
    1764 !------------------------------------------------------------------------------!
     1737!> Calculate static thermal indices for any point within a 2D grid time_integration.f90: 1065ff
     1738!--------------------------------------------------------------------------------------------------!
    17651739 SUBROUTINE bio_calculate_thermal_index_maps( av )
    17661740
     
    17741748
    17751749    REAL(wp) ::  clo          !< Clothing index                (no dimension)
     1750    REAL(wp) ::  pair         !< Air pressure                     (hPa)
     1751    REAL(wp) ::  perct_ij     !< Perceived temperature            (degree_C)
     1752    REAL(wp) ::  pet_ij       !< Physiologically equivalent temperature  (degree_C)
    17761753    REAL(wp) ::  ta           !< Air temperature                  (degree_C)
     1754    REAL(wp) ::  tmrt_ij      !< Mean radiant temperature         (degree_C)
     1755    REAL(wp) ::  utci_ij      !< Universal thermal climate index  (degree_C)
    17771756    REAL(wp) ::  vp           !< Vapour pressure                  (hPa)
    17781757    REAL(wp) ::  ws           !< Wind speed    (local level)      (m/s)
    1779     REAL(wp) ::  pair         !< Air pressure                     (hPa)
    1780     REAL(wp) ::  perct_ij     !< Perceived temperature            (degree_C)
    1781     REAL(wp) ::  utci_ij      !< Universal thermal climate index  (degree_C)
    1782     REAL(wp) ::  pet_ij       !< Physiologically equivalent temperature  (degree_C)
    1783     REAL(wp) ::  tmrt_ij      !< Mean radiant temperature         (degree_C)
    1784 
    1785 !
    1786 !-- Check if some thermal index is desired. Don't do anything if, e.g. only
    1787 !-- bio_mrt is desired.
    1788     IF ( do_calculate_perct  .OR.  do_calculate_perct_av  .OR.                 &
    1789        do_calculate_utci  .OR.  do_calculate_utci_av  .OR.                     &
    1790        do_calculate_pet  .OR.  do_calculate_pet_av  .OR.                       &
    1791        do_calculate_mrt2d )  THEN
     1758
     1759!
     1760!-- Check if some thermal index is desired. Don't do anything if, e.g. only bio_mrt is desired.
     1761    IF ( do_calculate_perct    .OR.  do_calculate_perct_av  .OR.  do_calculate_utci    .OR.        &
     1762         do_calculate_utci_av  .OR.  do_calculate_pet       .OR.  do_calculate_pet_av  .OR.        &
     1763         do_calculate_mrt2d )  THEN
    17921764
    17931765!
     
    18031775!
    18041776!--          Determine local meteorological conditions
    1805              CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,           &
    1806                                                    ws, pair, tmrt_ij )
     1777             CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp, ws, pair, tmrt_ij )
    18071778!
    18081779!--          Only proceed if input is available
     
    18101781             perct_ij = bio_fill_value   !< within some obstacle
    18111782             utci_ij  = bio_fill_value
    1812              IF ( .NOT. ( tmrt_ij <= -998._wp  .OR.  vp <= -998._wp  .OR.      &
    1813                 ws <= -998._wp  .OR.  ta <= -998._wp ) )  THEN
     1783             IF ( .NOT. ( tmrt_ij <= -998.0_wp  .OR.  vp <= -998.0_wp  .OR.   ws <= -998.0_wp  .OR.&
     1784                  ta <= -998.0_wp ) )  THEN
    18141785!
    18151786!--             Calculate static thermal indices based on local tmrt
     
    18191790!
    18201791!--                Estimate local perceived temperature
    1821                    CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair,     &
    1822                       clo, perct_ij )
     1792                   CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, clo, perct_ij )
    18231793                ENDIF
    18241794
     
    18261796!
    18271797!--                Estimate local universal thermal climate index
    1828                    CALL calculate_utci_static( ta, vp, ws, tmrt_ij,            &
    1829                       bio_output_height, utci_ij )
     1798                   CALL calculate_utci_static( ta, vp, ws, tmrt_ij, bio_output_height, utci_ij )
    18301799                ENDIF
    18311800
     
    18331802!
    18341803!--                Estimate local physiologically equivalent temperature
    1835                    CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair,       &
    1836                       pet_ij )
     1804                   CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij )
    18371805                ENDIF
    18381806             ENDIF
     
    18711839 END SUBROUTINE bio_calculate_thermal_index_maps
    18721840
    1873 !------------------------------------------------------------------------------!
     1841!--------------------------------------------------------------------------------------------------!
    18741842! Description:
    18751843! ------------
    18761844!> Calculate dynamic thermal indices (currently only iPT, but expandable)
    1877 !------------------------------------------------------------------------------!
    1878  SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,          &
    1879     t_clo, clo, actlev, age, weight, height, work, sex, ipt )
     1845!--------------------------------------------------------------------------------------------------!
     1846 SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage, t_clo, clo, actlev, age,     &
     1847                          weight, height, work, sex, ipt )
    18801848
    18811849    IMPLICIT NONE
    18821850!
    18831851!-- Input parameters
    1884     REAL(wp), INTENT ( IN )  ::  ta   !< Air temperature                  (degree_C)
    1885     REAL(wp), INTENT ( IN )  ::  vp   !< Vapour pressure                  (hPa)
    1886     REAL(wp), INTENT ( IN )  ::  ws   !< Wind speed    (local level)      (m/s)
    1887     REAL(wp), INTENT ( IN )  ::  pair !< Air pressure                     (hPa)
    1888     REAL(wp), INTENT ( IN )  ::  tmrt !< Mean radiant temperature         (degree_C)
    1889     REAL(wp), INTENT ( IN )  ::  dt   !< Time past since last calculation (s)
    1890     REAL(wp), INTENT ( IN )  ::  age  !< Age of agent                     (y)
    1891     REAL(wp), INTENT ( IN )  ::  weight  !< Weight of agent               (Kg)
    1892     REAL(wp), INTENT ( IN )  ::  height  !< Height of agent               (m)
    1893     REAL(wp), INTENT ( IN )  ::  work    !< Mechanical workload of agent
    1894                                          !  (without metabolism!)         (W)
    18951852    INTEGER(iwp), INTENT ( IN ) ::  sex  !< Sex of agent (1 = male, 2 = female)
     1853
     1854    REAL(wp), INTENT ( IN )  ::  age     !< Age of agent                     (y)
     1855    REAL(wp), INTENT ( IN )  ::  dt      !< Time past SINce last calculation (s)
     1856    REAL(wp), INTENT ( IN )  ::  height  !< Height of agent                  (m)
     1857    REAL(wp), INTENT ( IN )  ::  pair    !< Air pressure                     (hPa)
     1858    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature                  (degree_C)
     1859    REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature         (degree_C)
     1860    REAL(wp), INTENT ( IN )  ::  vp      !< Vapour pressure                  (hPa)
     1861    REAL(wp), INTENT ( IN )  ::  weight  !< Weight of agent                  (Kg)
     1862    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical workload of agent  (without metabolism!) (W)
     1863    REAL(wp), INTENT ( IN )  ::  ws      !< Wind speed  (local level)        (m/s)
     1864
    18961865!
    18971866!-- Both, input and output parameters
     1867    Real(wp), INTENT ( INOUT )  ::  actlev            !< Individuals activity level
     1868                                                      !< per unit surface area      (W/m²)
     1869    Real(wp), INTENT ( INOUT )  ::  clo               !< Current clothing in sulation
    18981870    Real(wp), INTENT ( INOUT )  ::  energy_storage    !< Energy storage   (W/m²)
    1899     Real(wp), INTENT ( INOUT )  ::  t_clo   !< Clothing temperature       (degree_C)
    1900     Real(wp), INTENT ( INOUT )  ::  clo     !< Current clothing in sulation
    1901     Real(wp), INTENT ( INOUT )  ::  actlev  !< Individuals activity level
    1902                                             !  per unit surface area      (W/m²)
     1871    Real(wp), INTENT ( INOUT )  ::  t_clo             !< Clothing temperature       (degree_C)
    19031872!
    19041873!-- Output parameters
    19051874    REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (degree_C)
    19061875!
    1907 !-- return immediatelly if nothing to do!
     1876!-- Return immediatelly if nothing to do!
    19081877    IF ( .NOT. thermal_comfort )  THEN
    19091878        RETURN
     
    19111880!
    19121881!-- If clo equals the initial value, this is the initial call
    1913     IF ( clo <= -998._wp )  THEN
    1914 !
    1915 !--    Initialize instationary perceived temperature with personalized
    1916 !      PT as an initial guess, set actlev and clo
    1917        CALL ipt_init( age, weight, height, sex, work, actlev, clo,            &
    1918           ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,                   &
    1919           ipt )
     1882    IF ( clo <= -998.0_wp )  THEN
     1883!
     1884!--    Initialize instationary perceived temperature with personalized PT as an initial guess, set
     1885!--    actlev and clo
     1886       CALL ipt_init( age, weight, height, sex, work, actlev, clo, ta, vp, ws, tmrt, pair, dt,     &
     1887                      energy_storage, t_clo, ipt )
    19201888    ELSE
    19211889!
    19221890!--    Estimate local instatinoary perceived temperature
    1923        CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,     &
    1924           clo, actlev, work, ipt )
     1891       CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo, clo, actlev, work, ipt )
    19251892    ENDIF
    19261893
     
    19291896
    19301897
    1931 !------------------------------------------------------------------------------!
     1898!--------------------------------------------------------------------------------------------------!
    19321899! Description:
    19331900! ------------
     
    19351902!> computed by a 6th order approximation
    19361903!>
    1937 !> UTCI regression equation after
    1938 !> Bröde P, Fiala D, Blazejczyk K, Holmér I, Jendritzky G, Kampmann B, Tinz B,
    1939 !> Havenith G (2012) Deriving the operational procedure for the Universal Thermal
    1940 !> Climate Index (UTCI). International Journal of Biometeorology 56 (3):481-494.
    1941 !> doi:10.1007/s00484-011-0454-1
     1904!> UTCI regression equation according to
     1905!> Bröde P, Fiala D, Blazejczyk K, Holmér I, Jendritzky G, Kampmann B, Tinz B, Havenith G (2012)
     1906!> Deriving the operational procedure for the Universal Thermal Climate Index (UTCI). International
     1907!> Journal of Biometeorology 56 (3):481-494. doi:10.1007/s00484-011-0454-1
    19421908!>
    19431909!> original source available at:
    19441910!> www.utci.org
    1945 !------------------------------------------------------------------------------!
     1911!--------------------------------------------------------------------------------------------------!
    19461912 SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci_ij )
    19471913
     
    19491915!
    19501916!-- Type of input of the argument list
     1917    REAL(WP), INTENT ( IN )  ::  hag      !< Height of wind speed input (m)
    19511918    REAL(WP), INTENT ( IN )  ::  ta_in    !< Local air temperature (degree_C)
     1919    REAL(WP), INTENT ( IN )  ::  tmrt     !< Local mean radiant temperature (degree_C)
    19521920    REAL(WP), INTENT ( IN )  ::  vp       !< Loacl vapour pressure (hPa)
    19531921    REAL(WP), INTENT ( IN )  ::  ws_hag   !< Incident wind speed (m/s)
    1954     REAL(WP), INTENT ( IN )  ::  tmrt     !< Local mean radiant temperature (degree_C)
    1955     REAL(WP), INTENT ( IN )  ::  hag      !< Height of wind speed input (m)
    19561922!
    19571923!-- Type of output of the argument list
    1958     REAL(wp), INTENT ( OUT ) ::  utci_ij  !< Universal Thermal Climate Index (degree_C)
    1959 
     1924    REAL(WP) ::  d_tmrt       !< delta-tmrt               (degree_C)
     1925    REAL(WP) ::  d_tmrt2      !< 2 times d_tmrt
     1926    REAL(WP) ::  d_tmrt3      !< 3 times d_tmrt
     1927    REAL(WP) ::  d_tmrt4      !< 4 times d_tmrt
     1928    REAL(WP) ::  d_tmrt5      !< 5 times d_tmrt
     1929    REAL(WP) ::  d_tmrt6      !< 6 times d_tmrt
     1930    REAL(WP) ::  offset       !< utci deviation by ta cond. exceeded      (degree_C)
     1931    REAL(WP) ::  pa           !< air pressure in kPa      (kPa)
     1932    REAL(WP) ::  pa2          !< 2 times pa
     1933    REAL(WP) ::  pa3          !< 3 times pa
     1934    REAL(WP) ::  pa4          !< 4 times pa
     1935    REAL(WP) ::  pa5          !< 5 times pa
     1936    REAL(WP) ::  pa6          !< 6 times pa
     1937    REAL(WP) ::  part_d_tmrt  !< Mean radiant temp. related part of the reg.
     1938    REAL(WP) ::  part_pa      !< Air pressure related part of the regression
     1939    REAL(WP) ::  part_pa2     !< Air pressure^2 related part of the regression
     1940    REAL(WP) ::  part_pa3     !< Air pressure^3 related part of the regression
     1941    REAL(WP) ::  part_pa46    !< Air pressure^4-6 related part of the regression
     1942    REAL(WP) ::  part_ta      !< Air temperature related part of the regression
     1943    REAL(WP) ::  part_va      !< Vapour pressure related part of the regression
    19601944    REAL(WP) ::  ta           !< air temperature modified by offset (degree_C)
    1961     REAL(WP) ::  pa           !< air pressure in kPa      (kPa)
    1962     REAL(WP) ::  d_tmrt       !< delta-tmrt               (degree_C)
    1963     REAL(WP) ::  va           !< wind speed at 10 m above ground level    (m/s)
    1964     REAL(WP) ::  offset       !< utci deviation by ta cond. exceeded      (degree_C)
    1965     REAL(WP) ::  part_ta      !< Air temperature related part of the regression
    19661945    REAL(WP) ::  ta2          !< 2 times ta
    19671946    REAL(WP) ::  ta3          !< 3 times ta
     
    19691948    REAL(WP) ::  ta5          !< 5 times ta
    19701949    REAL(WP) ::  ta6          !< 6 times ta
    1971     REAL(WP) ::  part_va      !< Vapour pressure related part of the regression
     1950    REAL(WP) ::  va           !< wind speed at 10 m above ground level    (m/s)
    19721951    REAL(WP) ::  va2          !< 2 times va
    19731952    REAL(WP) ::  va3          !< 3 times va
     
    19751954    REAL(WP) ::  va5          !< 5 times va
    19761955    REAL(WP) ::  va6          !< 6 times va
    1977     REAL(WP) ::  part_d_tmrt  !< Mean radiant temp. related part of the reg.
    1978     REAL(WP) ::  d_tmrt2      !< 2 times d_tmrt
    1979     REAL(WP) ::  d_tmrt3      !< 3 times d_tmrt
    1980     REAL(WP) ::  d_tmrt4      !< 4 times d_tmrt
    1981     REAL(WP) ::  d_tmrt5      !< 5 times d_tmrt
    1982     REAL(WP) ::  d_tmrt6      !< 6 times d_tmrt
    1983     REAL(WP) ::  part_pa      !< Air pressure related part of the regression
    1984     REAL(WP) ::  pa2          !< 2 times pa
    1985     REAL(WP) ::  pa3          !< 3 times pa
    1986     REAL(WP) ::  pa4          !< 4 times pa
    1987     REAL(WP) ::  pa5          !< 5 times pa
    1988     REAL(WP) ::  pa6          !< 6 times pa
    1989     REAL(WP) ::  part_pa2     !< Air pressure^2 related part of the regression
    1990     REAL(WP) ::  part_pa3     !< Air pressure^3 related part of the regression
    1991     REAL(WP) ::  part_pa46    !< Air pressure^4-6 related part of the regression
     1956
     1957
     1958    REAL(wp), INTENT ( OUT ) ::  utci_ij  !< Universal Thermal Climate Index (degree_C)
    19921959
    19931960!
    19941961!-- Initialize
    1995     offset = 0._wp
     1962    offset = 0.0_wp
    19961963    ta = ta_in
    19971964    d_tmrt = tmrt - ta_in
     
    20051972!
    20061973!-- Check if input values in range after Broede et al. (2012)
    2007     IF ( ( d_tmrt > 70._wp )  .OR.  ( d_tmrt < -30._wp )  .OR.                 &
    2008        ( vp >= 50._wp ) )  THEN
     1974    IF ( ( d_tmrt > 70.0_wp )  .OR.  ( d_tmrt < -30.0_wp )  .OR.  ( vp >= 50.0_wp ) )  THEN
    20091975       utci_ij = bio_fill_value
    20101976       RETURN
     
    20121978!
    20131979!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
    2014     IF ( ta > 50._wp )  THEN
    2015        offset = ta - 50._wp
    2016        ta = 50._wp
    2017     ENDIF
    2018     IF ( ta < -50._wp )  THEN
    2019        offset = ta + 50._wp
    2020        ta = -50._wp
    2021     ENDIF
    2022 !
    2023 !-- For routine application. For wind speeds and relative
    2024 !-- humidity values below 0.5 m/s or 5%, respectively, the
    2025 !-- user is advised to use the lower bounds for the calculations.
     1980    IF ( ta > 50.0_wp )  THEN
     1981       offset = ta - 50.0_wp
     1982       ta = 50.0_wp
     1983    ENDIF
     1984    IF ( ta < -50.0_wp )  THEN
     1985       offset = ta + 50.0_wp
     1986       ta = -50.0_wp
     1987    ENDIF
     1988!
     1989!-- For routine application. For wind speeds and relative humidity values below 0.5 m/s or 5%,
     1990!-- respectively, the user is advised to use the lower bounds for the calculations.
    20261991    IF ( va < 0.5_wp )  va = 0.5_wp
    2027     IF ( va > 17._wp )  va = 17._wp
     1992    IF ( va > 17.0_wp )  va = 17.0_wp
    20281993
    20291994!
    20301995!-- Pre-calculate multiples of input parameters to save time later
    2031 
    20321996    ta2 = ta  * ta
    20331997    ta3 = ta2 * ta
     
    20562020!
    20572021!-- Pre-calculate parts of the regression equation
    2058     part_ta = (  6.07562052e-01_wp )       +                                   &
    2059               ( -2.27712343e-02_wp ) * ta  +                                   &
    2060               (  8.06470249e-04_wp ) * ta2 +                                   &
    2061               ( -1.54271372e-04_wp ) * ta3 +                                   &
    2062               ( -3.24651735e-06_wp ) * ta4 +                                   &
    2063               (  7.32602852e-08_wp ) * ta5 +                                   &
     2022    part_ta = (  6.07562052e-01_wp )       +                                                       &
     2023              ( -2.27712343e-02_wp ) * ta  +                                                       &
     2024              (  8.06470249e-04_wp ) * ta2 +                                                       &
     2025              ( -1.54271372e-04_wp ) * ta3 +                                                       &
     2026              ( -3.24651735e-06_wp ) * ta4 +                                                       &
     2027              (  7.32602852e-08_wp ) * ta5 +                                                       &
    20642028              (  1.35959073e-09_wp ) * ta6
    20652029
    2066     part_va = ( -2.25836520e+00_wp ) * va +                                    &
    2067         (  8.80326035e-02_wp ) * ta  * va +                                    &
    2068         (  2.16844454e-03_wp ) * ta2 * va +                                    &
    2069         ( -1.53347087e-05_wp ) * ta3 * va +                                    &
    2070         ( -5.72983704e-07_wp ) * ta4 * va +                                    &
    2071         ( -2.55090145e-09_wp ) * ta5 * va +                                    &
    2072         ( -7.51269505e-01_wp ) *       va2 +                                   &
    2073         ( -4.08350271e-03_wp ) * ta  * va2 +                                   &
    2074         ( -5.21670675e-05_wp ) * ta2 * va2 +                                   &
    2075         (  1.94544667e-06_wp ) * ta3 * va2 +                                   &
    2076         (  1.14099531e-08_wp ) * ta4 * va2 +                                   &
    2077         (  1.58137256e-01_wp ) *       va3 +                                   &
    2078         ( -6.57263143e-05_wp ) * ta  * va3 +                                   &
    2079         (  2.22697524e-07_wp ) * ta2 * va3 +                                   &
    2080         ( -4.16117031e-08_wp ) * ta3 * va3 +                                   &
    2081         ( -1.27762753e-02_wp ) *       va4 +                                   &
    2082         (  9.66891875e-06_wp ) * ta  * va4 +                                   &
    2083         (  2.52785852e-09_wp ) * ta2 * va4 +                                   &
    2084         (  4.56306672e-04_wp ) *       va5 +                                   &
    2085         ( -1.74202546e-07_wp ) * ta  * va5 +                                   &
    2086         ( -5.91491269e-06_wp ) * va6
    2087 
    2088     part_d_tmrt = (  3.98374029e-01_wp ) *       d_tmrt +                      &
    2089             (  1.83945314e-04_wp ) * ta  *       d_tmrt +                      &
    2090             ( -1.73754510e-04_wp ) * ta2 *       d_tmrt +                      &
    2091             ( -7.60781159e-07_wp ) * ta3 *       d_tmrt +                      &
    2092             (  3.77830287e-08_wp ) * ta4 *       d_tmrt +                      &
    2093             (  5.43079673e-10_wp ) * ta5 *       d_tmrt +                      &
    2094             ( -2.00518269e-02_wp ) *       va  * d_tmrt +                      &
    2095             (  8.92859837e-04_wp ) * ta  * va  * d_tmrt +                      &
    2096             (  3.45433048e-06_wp ) * ta2 * va  * d_tmrt +                      &
    2097             ( -3.77925774e-07_wp ) * ta3 * va  * d_tmrt +                      &
    2098             ( -1.69699377e-09_wp ) * ta4 * va  * d_tmrt +                      &
    2099             (  1.69992415e-04_wp ) *       va2 * d_tmrt +                      &
    2100             ( -4.99204314e-05_wp ) * ta  * va2 * d_tmrt +                      &
    2101             (  2.47417178e-07_wp ) * ta2 * va2 * d_tmrt +                      &
    2102             (  1.07596466e-08_wp ) * ta3 * va2 * d_tmrt +                      &
    2103             (  8.49242932e-05_wp ) *       va3 * d_tmrt +                      &
    2104             (  1.35191328e-06_wp ) * ta  * va3 * d_tmrt +                      &
    2105             ( -6.21531254e-09_wp ) * ta2 * va3 * d_tmrt +                      &
    2106             ( -4.99410301e-06_wp ) * va4 *       d_tmrt +                      &
    2107             ( -1.89489258e-08_wp ) * ta  * va4 * d_tmrt +                      &
    2108             (  8.15300114e-08_wp ) *       va5 * d_tmrt +                      &
    2109             (  7.55043090e-04_wp ) *             d_tmrt2 +                     &
    2110             ( -5.65095215e-05_wp ) * ta  *       d_tmrt2 +                     &
    2111             ( -4.52166564e-07_wp ) * ta2 *       d_tmrt2 +                     &
    2112             (  2.46688878e-08_wp ) * ta3 *       d_tmrt2 +                     &
    2113             (  2.42674348e-10_wp ) * ta4 *       d_tmrt2 +                     &
    2114             (  1.54547250e-04_wp ) *       va  * d_tmrt2 +                     &
    2115             (  5.24110970e-06_wp ) * ta  * va  * d_tmrt2 +                     &
    2116             ( -8.75874982e-08_wp ) * ta2 * va  * d_tmrt2 +                     &
    2117             ( -1.50743064e-09_wp ) * ta3 * va  * d_tmrt2 +                     &
    2118             ( -1.56236307e-05_wp ) *       va2 * d_tmrt2 +                     &
    2119             ( -1.33895614e-07_wp ) * ta  * va2 * d_tmrt2 +                     &
    2120             (  2.49709824e-09_wp ) * ta2 * va2 * d_tmrt2 +                     &
    2121             (  6.51711721e-07_wp ) *       va3 * d_tmrt2 +                     &
    2122             (  1.94960053e-09_wp ) * ta  * va3 * d_tmrt2 +                     &
    2123             ( -1.00361113e-08_wp ) *       va4 * d_tmrt2 +                     &
    2124             ( -1.21206673e-05_wp ) *             d_tmrt3 +                     &
    2125             ( -2.18203660e-07_wp ) * ta  *       d_tmrt3 +                     &
    2126             (  7.51269482e-09_wp ) * ta2 *       d_tmrt3 +                     &
    2127             (  9.79063848e-11_wp ) * ta3 *       d_tmrt3 +                     &
    2128             (  1.25006734e-06_wp ) *       va  * d_tmrt3 +                     &
    2129             ( -1.81584736e-09_wp ) * ta  * va  * d_tmrt3 +                     &
    2130             ( -3.52197671e-10_wp ) * ta2 * va  * d_tmrt3 +                     &
    2131             ( -3.36514630e-08_wp ) *       va2 * d_tmrt3 +                     &
    2132             (  1.35908359e-10_wp ) * ta  * va2 * d_tmrt3 +                     &
    2133             (  4.17032620e-10_wp ) *       va3 * d_tmrt3 +                     &
    2134             ( -1.30369025e-09_wp ) *             d_tmrt4 +                     &
    2135             (  4.13908461e-10_wp ) * ta  *       d_tmrt4 +                     &
    2136             (  9.22652254e-12_wp ) * ta2 *       d_tmrt4 +                     &
    2137             ( -5.08220384e-09_wp ) *       va  * d_tmrt4 +                     &
    2138             ( -2.24730961e-11_wp ) * ta  * va  * d_tmrt4 +                     &
    2139             (  1.17139133e-10_wp ) *       va2 * d_tmrt4 +                     &
    2140             (  6.62154879e-10_wp ) *             d_tmrt5 +                     &
    2141             (  4.03863260e-13_wp ) * ta  *       d_tmrt5 +                     &
    2142             (  1.95087203e-12_wp ) *       va  * d_tmrt5 +                     &
    2143             ( -4.73602469e-12_wp ) *             d_tmrt6
    2144 
    2145     part_pa = (  5.12733497e+00_wp ) *                pa +                     &
    2146        ( -3.12788561e-01_wp ) * ta  *                 pa +                     &
    2147        ( -1.96701861e-02_wp ) * ta2 *                 pa +                     &
    2148        (  9.99690870e-04_wp ) * ta3 *                 pa +                     &
    2149        (  9.51738512e-06_wp ) * ta4 *                 pa +                     &
    2150        ( -4.66426341e-07_wp ) * ta5 *                 pa +                     &
    2151        (  5.48050612e-01_wp ) *       va  *           pa +                     &
    2152        ( -3.30552823e-03_wp ) * ta  * va  *           pa +                     &
    2153        ( -1.64119440e-03_wp ) * ta2 * va  *           pa +                     &
    2154        ( -5.16670694e-06_wp ) * ta3 * va  *           pa +                     &
    2155        (  9.52692432e-07_wp ) * ta4 * va  *           pa +                     &
    2156        ( -4.29223622e-02_wp ) *       va2 *           pa +                     &
    2157        (  5.00845667e-03_wp ) * ta  * va2 *           pa +                     &
    2158        (  1.00601257e-06_wp ) * ta2 * va2 *           pa +                     &
    2159        ( -1.81748644e-06_wp ) * ta3 * va2 *           pa +                     &
    2160        ( -1.25813502e-03_wp ) *       va3 *           pa +                     &
    2161        ( -1.79330391e-04_wp ) * ta  * va3 *           pa +                     &
    2162        (  2.34994441e-06_wp ) * ta2 * va3 *           pa +                     &
    2163        (  1.29735808e-04_wp ) *       va4 *           pa +                     &
    2164        (  1.29064870e-06_wp ) * ta  * va4 *           pa +                     &
    2165        ( -2.28558686e-06_wp ) *       va5 *           pa +                     &
    2166        ( -3.69476348e-02_wp ) *             d_tmrt  * pa +                     &
    2167        (  1.62325322e-03_wp ) * ta  *       d_tmrt  * pa +                     &
    2168        ( -3.14279680e-05_wp ) * ta2 *       d_tmrt  * pa +                     &
    2169        (  2.59835559e-06_wp ) * ta3 *       d_tmrt  * pa +                     &
    2170        ( -4.77136523e-08_wp ) * ta4 *       d_tmrt  * pa +                     &
    2171        (  8.64203390e-03_wp ) *       va  * d_tmrt  * pa +                     &
    2172        ( -6.87405181e-04_wp ) * ta  * va  * d_tmrt  * pa +                     &
    2173        ( -9.13863872e-06_wp ) * ta2 * va  * d_tmrt  * pa +                     &
    2174        (  5.15916806e-07_wp ) * ta3 * va  * d_tmrt  * pa +                     &
    2175        ( -3.59217476e-05_wp ) *       va2 * d_tmrt  * pa +                     &
    2176        (  3.28696511e-05_wp ) * ta  * va2 * d_tmrt  * pa +                     &
    2177        ( -7.10542454e-07_wp ) * ta2 * va2 * d_tmrt  * pa +                     &
    2178        ( -1.24382300e-05_wp ) *       va3 * d_tmrt  * pa +                     &
    2179        ( -7.38584400e-09_wp ) * ta  * va3 * d_tmrt  * pa +                     &
    2180        (  2.20609296e-07_wp ) *       va4 * d_tmrt  * pa +                     &
    2181        ( -7.32469180e-04_wp ) *             d_tmrt2 * pa +                     &
    2182        ( -1.87381964e-05_wp ) * ta  *       d_tmrt2 * pa +                     &
    2183        (  4.80925239e-06_wp ) * ta2 *       d_tmrt2 * pa +                     &
    2184        ( -8.75492040e-08_wp ) * ta3 *       d_tmrt2 * pa +                     &
    2185        (  2.77862930e-05_wp ) *       va  * d_tmrt2 * pa +                     &
    2186        ( -5.06004592e-06_wp ) * ta  * va  * d_tmrt2 * pa +                     &
    2187        (  1.14325367e-07_wp ) * ta2 * va  * d_tmrt2 * pa +                     &
    2188        (  2.53016723e-06_wp ) *       va2 * d_tmrt2 * pa +                     &
    2189        ( -1.72857035e-08_wp ) * ta  * va2 * d_tmrt2 * pa +                     &
    2190        ( -3.95079398e-08_wp ) *       va3 * d_tmrt2 * pa +                     &
    2191        ( -3.59413173e-07_wp ) *             d_tmrt3 * pa +                     &
    2192        (  7.04388046e-07_wp ) * ta  *       d_tmrt3 * pa +                     &
    2193        ( -1.89309167e-08_wp ) * ta2 *       d_tmrt3 * pa +                     &
    2194        ( -4.79768731e-07_wp ) *       va  * d_tmrt3 * pa +                     &
    2195        (  7.96079978e-09_wp ) * ta  * va  * d_tmrt3 * pa +                     &
    2196        (  1.62897058e-09_wp ) *       va2 * d_tmrt3 * pa +                     &
    2197        (  3.94367674e-08_wp ) *             d_tmrt4 * pa +                     &
    2198        ( -1.18566247e-09_wp ) * ta *        d_tmrt4 * pa +                     &
    2199        (  3.34678041e-10_wp ) *       va  * d_tmrt4 * pa +                     &
    2200        ( -1.15606447e-10_wp ) *             d_tmrt5 * pa
    2201 
    2202     part_pa2 = ( -2.80626406e+00_wp ) *               pa2 +                    &
    2203        (  5.48712484e-01_wp ) * ta  *                 pa2 +                    &
    2204        ( -3.99428410e-03_wp ) * ta2 *                 pa2 +                    &
    2205        ( -9.54009191e-04_wp ) * ta3 *                 pa2 +                    &
    2206        (  1.93090978e-05_wp ) * ta4 *                 pa2 +                    &
    2207        ( -3.08806365e-01_wp ) *       va *            pa2 +                    &
    2208        (  1.16952364e-02_wp ) * ta  * va *            pa2 +                    &
    2209        (  4.95271903e-04_wp ) * ta2 * va *            pa2 +                    &
    2210        ( -1.90710882e-05_wp ) * ta3 * va *            pa2 +                    &
    2211        (  2.10787756e-03_wp ) *       va2 *           pa2 +                    &
    2212        ( -6.98445738e-04_wp ) * ta  * va2 *           pa2 +                    &
    2213        (  2.30109073e-05_wp ) * ta2 * va2 *           pa2 +                    &
    2214        (  4.17856590e-04_wp ) *       va3 *           pa2 +                    &
    2215        ( -1.27043871e-05_wp ) * ta  * va3 *           pa2 +                    &
    2216        ( -3.04620472e-06_wp ) *       va4 *           pa2 +                    &
    2217        (  5.14507424e-02_wp ) *             d_tmrt  * pa2 +                    &
    2218        ( -4.32510997e-03_wp ) * ta  *       d_tmrt  * pa2 +                    &
    2219        (  8.99281156e-05_wp ) * ta2 *       d_tmrt  * pa2 +                    &
    2220        ( -7.14663943e-07_wp ) * ta3 *       d_tmrt  * pa2 +                    &
    2221        ( -2.66016305e-04_wp ) *       va  * d_tmrt  * pa2 +                    &
    2222        (  2.63789586e-04_wp ) * ta  * va  * d_tmrt  * pa2 +                    &
    2223        ( -7.01199003e-06_wp ) * ta2 * va  * d_tmrt  * pa2 +                    &
    2224        ( -1.06823306e-04_wp ) *       va2 * d_tmrt  * pa2 +                    &
    2225        (  3.61341136e-06_wp ) * ta  * va2 * d_tmrt  * pa2 +                    &
    2226        (  2.29748967e-07_wp ) *       va3 * d_tmrt  * pa2 +                    &
    2227        (  3.04788893e-04_wp ) *             d_tmrt2 * pa2 +                    &
    2228        ( -6.42070836e-05_wp ) * ta  *       d_tmrt2 * pa2 +                    &
    2229        (  1.16257971e-06_wp ) * ta2 *       d_tmrt2 * pa2 +                    &
    2230        (  7.68023384e-06_wp ) *       va  * d_tmrt2 * pa2 +                    &
    2231        ( -5.47446896e-07_wp ) * ta  * va  * d_tmrt2 * pa2 +                    &
    2232        ( -3.59937910e-08_wp ) *       va2 * d_tmrt2 * pa2 +                    &
    2233        ( -4.36497725e-06_wp ) *             d_tmrt3 * pa2 +                    &
    2234        (  1.68737969e-07_wp ) * ta  *       d_tmrt3 * pa2 +                    &
    2235        (  2.67489271e-08_wp ) *       va  * d_tmrt3 * pa2 +                    &
    2236        (  3.23926897e-09_wp ) *             d_tmrt4 * pa2
    2237 
    2238     part_pa3 = ( -3.53874123e-02_wp ) *               pa3 +                    &
    2239        ( -2.21201190e-01_wp ) * ta  *                 pa3 +                    &
    2240        (  1.55126038e-02_wp ) * ta2 *                 pa3 +                    &
    2241        ( -2.63917279e-04_wp ) * ta3 *                 pa3 +                    &
    2242        (  4.53433455e-02_wp ) *       va  *           pa3 +                    &
    2243        ( -4.32943862e-03_wp ) * ta  * va  *           pa3 +                    &
    2244        (  1.45389826e-04_wp ) * ta2 * va  *           pa3 +                    &
    2245        (  2.17508610e-04_wp ) *       va2 *           pa3 +                    &
    2246        ( -6.66724702e-05_wp ) * ta  * va2 *           pa3 +                    &
    2247        (  3.33217140e-05_wp ) *       va3 *           pa3 +                    &
    2248        ( -2.26921615e-03_wp ) *             d_tmrt  * pa3 +                    &
    2249        (  3.80261982e-04_wp ) * ta  *       d_tmrt  * pa3 +                    &
    2250        ( -5.45314314e-09_wp ) * ta2 *       d_tmrt  * pa3 +                    &
    2251        ( -7.96355448e-04_wp ) *       va  * d_tmrt  * pa3 +                    &
    2252        (  2.53458034e-05_wp ) * ta  * va  * d_tmrt  * pa3 +                    &
    2253        ( -6.31223658e-06_wp ) *       va2 * d_tmrt  * pa3 +                    &
    2254        (  3.02122035e-04_wp ) *             d_tmrt2 * pa3 +                    &
    2255        ( -4.77403547e-06_wp ) * ta  *       d_tmrt2 * pa3 +                    &
    2256        (  1.73825715e-06_wp ) *       va  * d_tmrt2 * pa3 +                    &
    2257        ( -4.09087898e-07_wp ) *             d_tmrt3 * pa3
    2258 
    2259     part_pa46 = (  6.14155345e-01_wp ) *              pa4 +                    &
    2260        ( -6.16755931e-02_wp ) * ta  *                 pa4 +                    &
    2261        (  1.33374846e-03_wp ) * ta2 *                 pa4 +                    &
    2262        (  3.55375387e-03_wp ) *       va  *           pa4 +                    &
    2263        ( -5.13027851e-04_wp ) * ta  * va  *           pa4 +                    &
    2264        (  1.02449757e-04_wp ) *       va2 *           pa4 +                    &
    2265        ( -1.48526421e-03_wp ) *             d_tmrt  * pa4 +                    &
    2266        ( -4.11469183e-05_wp ) * ta  *       d_tmrt  * pa4 +                    &
    2267        ( -6.80434415e-06_wp ) *       va  * d_tmrt  * pa4 +                    &
    2268        ( -9.77675906e-06_wp ) *             d_tmrt2 * pa4 +                    &
    2269        (  8.82773108e-02_wp ) *                       pa5 +                    &
    2270        ( -3.01859306e-03_wp ) * ta  *                 pa5 +                    &
    2271        (  1.04452989e-03_wp ) *       va  *           pa5 +                    &
    2272        (  2.47090539e-04_wp ) *             d_tmrt  * pa5 +                    &
    2273        (  1.48348065e-03_wp ) *                       pa6
    2274 !
    2275 !-- Calculate 6th order polynomial as approximation
    2276     utci_ij = ta + part_ta + part_va + part_d_tmrt + part_pa + part_pa2 +      &
    2277         part_pa3 + part_pa46
     2030    part_va = ( -2.25836520e+00_wp ) *       va  +                                                 &
     2031              (  8.80326035e-02_wp ) * ta  * va  +                                                 &
     2032              (  2.16844454e-03_wp ) * ta2 * va  +                                                 &
     2033              ( -1.53347087e-05_wp ) * ta3 * va  +                                                 &
     2034              ( -5.72983704e-07_wp ) * ta4 * va  +                                                 &
     2035              ( -2.55090145e-09_wp ) * ta5 * va  +                                                 &
     2036              ( -7.51269505e-01_wp ) *       va2 +                                                 &
     2037              ( -4.08350271e-03_wp ) * ta  * va2 +                                                 &
     2038              ( -5.21670675e-05_wp ) * ta2 * va2 +                                                 &
     2039              (  1.94544667e-06_wp ) * ta3 * va2 +                                                 &
     2040              (  1.14099531e-08_wp ) * ta4 * va2 +                                                 &
     2041              (  1.58137256e-01_wp ) *       va3 +                                                 &
     2042              ( -6.57263143e-05_wp ) * ta  * va3 +                                                 &
     2043              (  2.22697524e-07_wp ) * ta2 * va3 +                                                 &
     2044              ( -4.16117031e-08_wp ) * ta3 * va3 +                                                 &
     2045              ( -1.27762753e-02_wp ) *       va4 +                                                 &
     2046              (  9.66891875e-06_wp ) * ta  * va4 +                                                 &
     2047              (  2.52785852e-09_wp ) * ta2 * va4 +                                                 &
     2048              (  4.56306672e-04_wp ) *       va5 +                                                 &
     2049              ( -1.74202546e-07_wp ) * ta  * va5 +                                                 &
     2050              ( -5.91491269e-06_wp ) * va6
     2051
     2052    part_d_tmrt = (  3.98374029e-01_wp ) *             d_tmrt  +                                   &
     2053                  (  1.83945314e-04_wp ) * ta  *       d_tmrt  +                                   &
     2054                  ( -1.73754510e-04_wp ) * ta2 *       d_tmrt  +                                   &
     2055                  ( -7.60781159e-07_wp ) * ta3 *       d_tmrt  +                                   &
     2056                  (  3.77830287e-08_wp ) * ta4 *       d_tmrt  +                                   &
     2057                  (  5.43079673e-10_wp ) * ta5 *       d_tmrt  +                                   &
     2058                  ( -2.00518269e-02_wp ) *       va  * d_tmrt  +                                   &
     2059                  (  8.92859837e-04_wp ) * ta  * va  * d_tmrt  +                                   &
     2060                  (  3.45433048e-06_wp ) * ta2 * va  * d_tmrt  +                                   &
     2061                  ( -3.77925774e-07_wp ) * ta3 * va  * d_tmrt  +                                   &
     2062                  ( -1.69699377e-09_wp ) * ta4 * va  * d_tmrt  +                                   &
     2063                  (  1.69992415e-04_wp ) *       va2 * d_tmrt  +                                   &
     2064                  ( -4.99204314e-05_wp ) * ta  * va2 * d_tmrt  +                                   &
     2065                  (  2.47417178e-07_wp ) * ta2 * va2 * d_tmrt  +                                   &
     2066                  (  1.07596466e-08_wp ) * ta3 * va2 * d_tmrt  +                                   &
     2067                  (  8.49242932e-05_wp ) *       va3 * d_tmrt  +                                   &
     2068                  (  1.35191328e-06_wp ) * ta  * va3 * d_tmrt  +                                   &
     2069                  ( -6.21531254e-09_wp ) * ta2 * va3 * d_tmrt  +                                   &
     2070                  ( -4.99410301e-06_wp ) * va4 *       d_tmrt  +                                   &
     2071                  ( -1.89489258e-08_wp ) * ta  * va4 * d_tmrt  +                                   &
     2072                  (  8.15300114e-08_wp ) *       va5 * d_tmrt  +                                   &
     2073                  (  7.55043090e-04_wp ) *             d_tmrt2 +                                   &
     2074                  ( -5.65095215e-05_wp ) * ta  *       d_tmrt2 +                                   &
     2075                  ( -4.52166564e-07_wp ) * ta2 *       d_tmrt2 +                                   &
     2076                  (  2.46688878e-08_wp ) * ta3 *       d_tmrt2 +                                   &
     2077                  (  2.42674348e-10_wp ) * ta4 *       d_tmrt2 +                                   &
     2078                  (  1.54547250e-04_wp ) *       va  * d_tmrt2 +                                   &
     2079                  (  5.24110970e-06_wp ) * ta  * va  * d_tmrt2 +                                   &
     2080                  ( -8.75874982e-08_wp ) * ta2 * va  * d_tmrt2 +                                   &
     2081                  ( -1.50743064e-09_wp ) * ta3 * va  * d_tmrt2 +                                   &
     2082                  ( -1.56236307e-05_wp ) *       va2 * d_tmrt2 +                                   &
     2083                  ( -1.33895614e-07_wp ) * ta  * va2 * d_tmrt2 +                                   &
     2084                  (  2.49709824e-09_wp ) * ta2 * va2 * d_tmrt2 +                                   &
     2085                  (  6.51711721e-07_wp ) *       va3 * d_tmrt2 +                                   &
     2086                  (  1.94960053e-09_wp ) * ta  * va3 * d_tmrt2 +                                   &
     2087                  ( -1.00361113e-08_wp ) *       va4 * d_tmrt2 +                                   &
     2088                  ( -1.21206673e-05_wp ) *             d_tmrt3 +                                   &
     2089                  ( -2.18203660e-07_wp ) * ta  *       d_tmrt3 +                                   &
     2090                  (  7.51269482e-09_wp ) * ta2 *       d_tmrt3 +                                   &
     2091                  (  9.79063848e-11_wp ) * ta3 *       d_tmrt3 +                                   &
     2092                  (  1.25006734e-06_wp ) *       va  * d_tmrt3 +                                   &
     2093                  ( -1.81584736e-09_wp ) * ta  * va  * d_tmrt3 +                                   &
     2094                  ( -3.52197671e-10_wp ) * ta2 * va  * d_tmrt3 +                                   &
     2095                  ( -3.36514630e-08_wp ) *       va2 * d_tmrt3 +                                   &
     2096                  (  1.35908359e-10_wp ) * ta  * va2 * d_tmrt3 +                                   &
     2097                  (  4.17032620e-10_wp ) *       va3 * d_tmrt3 +                                   &
     2098                  ( -1.30369025e-09_wp ) *             d_tmrt4 +                                   &
     2099                  (  4.13908461e-10_wp ) * ta  *       d_tmrt4 +                                   &
     2100                  (  9.22652254e-12_wp ) * ta2 *       d_tmrt4 +                                   &
     2101                  ( -5.08220384e-09_wp ) *       va  * d_tmrt4 +                                   &
     2102                  ( -2.24730961e-11_wp ) * ta  * va  * d_tmrt4 +                                   &
     2103                  (  1.17139133e-10_wp ) *       va2 * d_tmrt4 +                                   &
     2104                  (  6.62154879e-10_wp ) *             d_tmrt5 +                                   &
     2105                  (  4.03863260e-13_wp ) * ta  *       d_tmrt5 +                                   &
     2106                  (  1.95087203e-12_wp ) *       va  * d_tmrt5 +                                   &
     2107                  ( -4.73602469e-12_wp ) *             d_tmrt6
     2108
     2109    part_pa = (  5.12733497e+00_wp ) *                       pa +                                  &
     2110              ( -3.12788561e-01_wp ) * ta  *                 pa +                                  &
     2111              ( -1.96701861e-02_wp ) * ta2 *                 pa +                                  &
     2112              (  9.99690870e-04_wp ) * ta3 *                 pa +                                  &
     2113              (  9.51738512e-06_wp ) * ta4 *                 pa +                                  &
     2114              ( -4.66426341e-07_wp ) * ta5 *                 pa +                                  &
     2115              (  5.48050612e-01_wp ) *       va  *           pa +                                  &
     2116              ( -3.30552823e-03_wp ) * ta  * va  *           pa +                                  &
     2117              ( -1.64119440e-03_wp ) * ta2 * va  *           pa +                                  &
     2118              ( -5.16670694e-06_wp ) * ta3 * va  *           pa +                                  &
     2119              (  9.52692432e-07_wp ) * ta4 * va  *           pa +                                  &
     2120              ( -4.29223622e-02_wp ) *       va2 *           pa +                                  &
     2121              (  5.00845667e-03_wp ) * ta  * va2 *           pa +                                  &
     2122              (  1.00601257e-06_wp ) * ta2 * va2 *           pa +                                  &
     2123              ( -1.81748644e-06_wp ) * ta3 * va2 *           pa +                                  &
     2124              ( -1.25813502e-03_wp ) *       va3 *           pa +                                  &
     2125              ( -1.79330391e-04_wp ) * ta  * va3 *           pa +                                  &
     2126              (  2.34994441e-06_wp ) * ta2 * va3 *           pa +                                  &
     2127              (  1.29735808e-04_wp ) *       va4 *           pa +                                  &
     2128              (  1.29064870e-06_wp ) * ta  * va4 *           pa +                                  &
     2129              ( -2.28558686e-06_wp ) *       va5 *           pa +                                  &
     2130              ( -3.69476348e-02_wp ) *             d_tmrt  * pa +                                  &
     2131              (  1.62325322e-03_wp ) * ta  *       d_tmrt  * pa +                                  &
     2132              ( -3.14279680e-05_wp ) * ta2 *       d_tmrt  * pa +                                  &
     2133              (  2.59835559e-06_wp ) * ta3 *       d_tmrt  * pa +                                  &
     2134              ( -4.77136523e-08_wp ) * ta4 *       d_tmrt  * pa +                                  &
     2135              (  8.64203390e-03_wp ) *       va  * d_tmrt  * pa +                                  &
     2136              ( -6.87405181e-04_wp ) * ta  * va  * d_tmrt  * pa +                                  &
     2137              ( -9.13863872e-06_wp ) * ta2 * va  * d_tmrt  * pa +                                  &
     2138              (  5.15916806e-07_wp ) * ta3 * va  * d_tmrt  * pa +                                  &
     2139              ( -3.59217476e-05_wp ) *       va2 * d_tmrt  * pa +                                  &
     2140              (  3.28696511e-05_wp ) * ta  * va2 * d_tmrt  * pa +                                  &
     2141              ( -7.10542454e-07_wp ) * ta2 * va2 * d_tmrt  * pa +                                  &
     2142              ( -1.24382300e-05_wp ) *       va3 * d_tmrt  * pa +                                  &
     2143              ( -7.38584400e-09_wp ) * ta  * va3 * d_tmrt  * pa +                                  &
     2144              (  2.20609296e-07_wp ) *       va4 * d_tmrt  * pa +                                  &
     2145              ( -7.32469180e-04_wp ) *             d_tmrt2 * pa +                                  &
     2146              ( -1.87381964e-05_wp ) * ta  *       d_tmrt2 * pa +                                  &
     2147              (  4.80925239e-06_wp ) * ta2 *       d_tmrt2 * pa +                                  &
     2148              ( -8.75492040e-08_wp ) * ta3 *       d_tmrt2 * pa +                                  &
     2149              (  2.77862930e-05_wp ) *       va  * d_tmrt2 * pa +                                  &
     2150              ( -5.06004592e-06_wp ) * ta  * va  * d_tmrt2 * pa +                                  &
     2151              (  1.14325367e-07_wp ) * ta2 * va  * d_tmrt2 * pa +                                  &
     2152              (  2.53016723e-06_wp ) *       va2 * d_tmrt2 * pa +                                  &
     2153              ( -1.72857035e-08_wp ) * ta  * va2 * d_tmrt2 * pa +                                  &
     2154              ( -3.95079398e-08_wp ) *       va3 * d_tmrt2 * pa +                                  &
     2155              ( -3.59413173e-07_wp ) *             d_tmrt3 * pa +                                  &
     2156              (  7.04388046e-07_wp ) * ta  *       d_tmrt3 * pa +                                  &
     2157              ( -1.89309167e-08_wp ) * ta2 *       d_tmrt3 * pa +                                  &
     2158              ( -4.79768731e-07_wp ) *       va  * d_tmrt3 * pa +                                  &
     2159              (  7.96079978e-09_wp ) * ta  * va  * d_tmrt3 * pa +                                  &
     2160              (  1.62897058e-09_wp ) *       va2 * d_tmrt3 * pa +                                  &
     2161              (  3.94367674e-08_wp ) *             d_tmrt4 * pa +                                  &
     2162              ( -1.18566247e-09_wp ) * ta *        d_tmrt4 * pa +                                  &
     2163              (  3.34678041e-10_wp ) *       va  * d_tmrt4 * pa +                                  &
     2164              ( -1.15606447e-10_wp ) *             d_tmrt5 * pa
     2165
     2166    part_pa2 = ( -2.80626406e+00_wp ) *                       pa2 +                                &
     2167               (  5.48712484e-01_wp ) * ta  *                 pa2 +                                &
     2168               ( -3.99428410e-03_wp ) * ta2 *                 pa2 +                                &
     2169               ( -9.54009191e-04_wp ) * ta3 *                 pa2 +                                &
     2170               (  1.93090978e-05_wp ) * ta4 *                 pa2 +                                &
     2171               ( -3.08806365e-01_wp ) *       va *            pa2 +                                &
     2172               (  1.16952364e-02_wp ) * ta  * va *            pa2 +                                &
     2173               (  4.95271903e-04_wp ) * ta2 * va *            pa2 +                                &
     2174               ( -1.90710882e-05_wp ) * ta3 * va *            pa2 +                                &
     2175               (  2.10787756e-03_wp ) *       va2 *           pa2 +                                &
     2176               ( -6.98445738e-04_wp ) * ta  * va2 *           pa2 +                                &
     2177               (  2.30109073e-05_wp ) * ta2 * va2 *           pa2 +                                &
     2178               (  4.17856590e-04_wp ) *       va3 *           pa2 +                                &
     2179               ( -1.27043871e-05_wp ) * ta  * va3 *           pa2 +                                &
     2180               ( -3.04620472e-06_wp ) *       va4 *           pa2 +                                &
     2181               (  5.14507424e-02_wp ) *             d_tmrt  * pa2 +                                &
     2182               ( -4.32510997e-03_wp ) * ta  *       d_tmrt  * pa2 +                                &
     2183               (  8.99281156e-05_wp ) * ta2 *       d_tmrt  * pa2 +                                &
     2184               ( -7.14663943e-07_wp ) * ta3 *       d_tmrt  * pa2 +                                &
     2185               ( -2.66016305e-04_wp ) *       va  * d_tmrt  * pa2 +                                &
     2186               (  2.63789586e-04_wp ) * ta  * va  * d_tmrt  * pa2 +                                &
     2187               ( -7.01199003e-06_wp ) * ta2 * va  * d_tmrt  * pa2 +                                &
     2188               ( -1.06823306e-04_wp ) *       va2 * d_tmrt  * pa2 +                                &
     2189               (  3.61341136e-06_wp ) * ta  * va2 * d_tmrt  * pa2 +                                &
     2190               (  2.29748967e-07_wp ) *       va3 * d_tmrt  * pa2 +                                &
     2191               (  3.04788893e-04_wp ) *             d_tmrt2 * pa2 +                                &
     2192               ( -6.42070836e-05_wp ) * ta  *       d_tmrt2 * pa2 +                                &
     2193               (  1.16257971e-06_wp ) * ta2 *       d_tmrt2 * pa2 +                                &
     2194               (  7.68023384e-06_wp ) *       va  * d_tmrt2 * pa2 +                                &
     2195               ( -5.47446896e-07_wp ) * ta  * va  * d_tmrt2 * pa2 +                                &
     2196               ( -3.59937910e-08_wp ) *       va2 * d_tmrt2 * pa2 +                                &
     2197               ( -4.36497725e-06_wp ) *             d_tmrt3 * pa2 +                                &
     2198               (  1.68737969e-07_wp ) * ta  *       d_tmrt3 * pa2 +                                &
     2199               (  2.67489271e-08_wp ) *       va  * d_tmrt3 * pa2 +                                &
     2200               (  3.23926897e-09_wp ) *             d_tmrt4 * pa2
     2201
     2202    part_pa3 = ( -3.53874123e-02_wp ) *                       pa3 +                                &
     2203               ( -2.21201190e-01_wp ) * ta  *                 pa3 +                                &
     2204               (  1.55126038e-02_wp ) * ta2 *                 pa3 +                                &
     2205               ( -2.63917279e-04_wp ) * ta3 *                 pa3 +                                &
     2206               (  4.53433455e-02_wp ) *       va  *           pa3 +                                &
     2207               ( -4.32943862e-03_wp ) * ta  * va  *           pa3 +                                &
     2208               (  1.45389826e-04_wp ) * ta2 * va  *           pa3 +                                &
     2209               (  2.17508610e-04_wp ) *       va2 *           pa3 +                                &
     2210               ( -6.66724702e-05_wp ) * ta  * va2 *           pa3 +                                &
     2211               (  3.33217140e-05_wp ) *       va3 *           pa3 +                                &
     2212               ( -2.26921615e-03_wp ) *             d_tmrt  * pa3 +                                &
     2213               (  3.80261982e-04_wp ) * ta  *       d_tmrt  * pa3 +                                &
     2214               ( -5.45314314e-09_wp ) * ta2 *       d_tmrt  * pa3 +                                &
     2215               ( -7.96355448e-04_wp ) *       va  * d_tmrt  * pa3 +                                &
     2216               (  2.53458034e-05_wp ) * ta  * va  * d_tmrt  * pa3 +                                &
     2217               ( -6.31223658e-06_wp ) *       va2 * d_tmrt  * pa3 +                                &
     2218               (  3.02122035e-04_wp ) *             d_tmrt2 * pa3 +                                &
     2219               ( -4.77403547e-06_wp ) * ta  *       d_tmrt2 * pa3 +                                &
     2220               (  1.73825715e-06_wp ) *       va  * d_tmrt2 * pa3 +                                &
     2221               ( -4.09087898e-07_wp ) *             d_tmrt3 * pa3
     2222
     2223    part_pa46 = (  6.14155345e-01_wp ) *                       pa4 +                               &
     2224                ( -6.16755931e-02_wp ) * ta  *                 pa4 +                               &
     2225                (  1.33374846e-03_wp ) * ta2 *                 pa4 +                               &
     2226                (  3.55375387e-03_wp ) *       va  *           pa4 +                               &
     2227                ( -5.13027851e-04_wp ) * ta  * va  *           pa4 +                               &
     2228                (  1.02449757e-04_wp ) *       va2 *           pa4 +                               &
     2229                ( -1.48526421e-03_wp ) *             d_tmrt  * pa4 +                               &
     2230                ( -4.11469183e-05_wp ) * ta  *       d_tmrt  * pa4 +                               &
     2231                ( -6.80434415e-06_wp ) *       va  * d_tmrt  * pa4 +                               &
     2232                ( -9.77675906e-06_wp ) *             d_tmrt2 * pa4 +                               &
     2233                (  8.82773108e-02_wp ) *                       pa5 +                               &
     2234                ( -3.01859306e-03_wp ) * ta  *                 pa5 +                               &
     2235                (  1.04452989e-03_wp ) *       va  *           pa5 +                               &
     2236                (  2.47090539e-04_wp ) *             d_tmrt  * pa5 +                               &
     2237                (  1.48348065e-03_wp ) *                       pa6
     2238!
     2239!-- Calculate 6th order polynomial as approximation
     2240    utci_ij = ta + part_ta + part_va + part_d_tmrt + part_pa + part_pa2 + part_pa3 + part_pa46
    22782241!
    22792242!-- Consider offset in result
     
    22852248
    22862249
    2287 !------------------------------------------------------------------------------!
     2250!--------------------------------------------------------------------------------------------------!
    22882251! Description:
    22892252! ------------
    2290 !> calculate_perct_static: Estimation of perceived temperature (PT, degree_C)
     2253!> Calculate_perct_static: Estimation of perceived temperature (PT, degree_C)
    22912254!> Value of perct is the Perceived Temperature, degree centigrade
    2292 !------------------------------------------------------------------------------!
     2255!--------------------------------------------------------------------------------------------------!
    22932256 SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct_ij )
    22942257
     
    22962259!
    22972260!-- Type of input of the argument list
     2261    REAL(wp), INTENT ( IN )  :: pair !< Local barometric air pressure (hPa)
    22982262    REAL(wp), INTENT ( IN )  :: ta   !< Local air temperature (degC)
     2263    REAL(wp), INTENT ( IN )  :: tmrt !< Local mean radiant temperature (degC)
    22992264    REAL(wp), INTENT ( IN )  :: vp   !< Local vapour pressure (hPa)
    2300     REAL(wp), INTENT ( IN )  :: tmrt !< Local mean radiant temperature (degC)
    23012265    REAL(wp), INTENT ( IN )  :: ws   !< Local wind velocitry (m/s)
    2302     REAL(wp), INTENT ( IN )  :: pair !< Local barometric air pressure (hPa)
    23032266!
    23042267!-- Type of output of the argument list
     2268    REAL(wp), INTENT ( OUT ) :: clo       !< Clothing index (dimensionless)
    23052269    REAL(wp), INTENT ( OUT ) :: perct_ij  !< Perceived temperature (degC)
    2306     REAL(wp), INTENT ( OUT ) :: clo       !< Clothing index (dimensionless)
    23072270!
    23082271!-- Parameters for standard "Klima-Michel"
    2309     REAL(wp), PARAMETER :: eta = 0._wp  !< Mechanical work efficiency for walking on flat ground
    2310                                         !< (compare to Fanger (1972) pp 24f)
    2311     REAL(wp), PARAMETER :: actlev = 134.6862_wp  !< Workload by activity per standardized surface (A_Du)
     2272    REAL(wp), PARAMETER :: actlev = 134.6862_wp  !< Workload by activity per standardized surface
     2273                                                 !< (A_Du)
     2274    REAL(wp), PARAMETER :: eta = 0.0_wp          !< Mechanical work efficiency for walking on flat
     2275                                                 !< ground (compare to Fanger (1972) pp 24f)
    23122276!
    23132277!-- Type of program variables
    2314     REAL(wp), PARAMETER :: eps = 0.0005  !< Accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
    2315     REAL(wp) ::  sclo           !< summer clothing insulation
    2316     REAL(wp) ::  wclo           !< winter clothing insulation
    2317     REAL(wp) ::  d_pmv          !< PMV deviation (dimensionless --> PMV)
    2318     REAL(wp) ::  svp_ta         !< saturation vapor pressure    (hPa)
    2319     REAL(wp) ::  sult_lim       !< threshold for sultrieness    (hPa)
    2320     REAL(wp) ::  dgtcm          !< Mean deviation dependent on perct
    2321     REAL(wp) ::  dgtcstd        !< Mean deviation plus its standard deviation
    2322     REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
    2323     REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
    2324     REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
    2325     REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
    2326     REAL(wp) ::  pmva           !< adjusted predicted mean vote
    2327     REAL(wp) ::  ptc            !< perceived temp. for cold conditions (degree_C)
    2328     REAL(wp) ::  d_std          !< factor to threshold for sultriness
    2329     REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
    2330 
    23312278    INTEGER(iwp) :: ncount      !< running index
    23322279    INTEGER(iwp) :: nerr_cold   !< error number (cold conditions)
     
    23342281
    23352282    LOGICAL :: sultrieness
     2283
     2284    REAL(wp), PARAMETER :: eps = 0.0005  !< Accuracy in clothing insulation (clo) for evaluation the
     2285                                         !< root of Fanger's PMV (pmva=0)
     2286
     2287    REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
     2288    REAL(wp) ::  d_pmv          !< PMV deviation (dimensionless --> PMV)
     2289    REAL(wp) ::  dgtcm          !< Mean deviation dependent on perct
     2290    REAL(wp) ::  dgtcstd        !< Mean deviation plus its standard deviation
     2291    REAL(wp) ::  d_std          !< factor to threshold for sultriness
     2292    REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
     2293    REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
     2294    REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
     2295    REAL(wp) ::  pmva           !< adjusted predicted mean vote
     2296    REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
     2297    REAL(wp) ::  ptc            !< perceived temp. for cold conditions (degree_C)
     2298    REAL(wp) ::  sclo           !< summer clothing insulation
     2299    REAL(wp) ::  svp_ta         !< saturation vapor pressure    (hPa)
     2300    REAL(wp) ::  sult_lim       !< threshold for sultrieness    (hPa)
     2301    REAL(wp) ::  wclo           !< winter clothing insulation
     2302
    23362303!
    23372304!-- Initialise
     
    23432310!
    23442311!-- Tresholds: clothing insulation (account for model inaccuracies)
    2345 !
    2346 !-- summer clothing
     2312!-- Summer clothing
    23472313    sclo     = 0.44453_wp
    23482314!
    2349 !-- winter clothing
     2315!-- Winter clothing
    23502316    wclo     = 1.76267_wp
    23512317!
    2352 !-- decision: firstly calculate for winter or summer clothing
    2353     IF ( ta <= 10._wp )  THEN
     2318!-- Eecision: first calculate for winter or summer clothing
     2319    IF ( ta <= 10.0_wp )  THEN
    23542320!
    23552321!--    First guess: winter clothing insulation: cold stress
     
    23582324       pmv_w = pmva
    23592325
    2360        IF ( pmva > 0._wp )  THEN
     2326       IF ( pmva > 0.0_wp )  THEN
    23612327!
    23622328!--       Case summer clothing insulation: heat load ?
     
    23642330          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
    23652331          pmv_s = pmva
    2366           IF ( pmva <= 0._wp )  THEN
    2367 !
    2368 !--          Case: comfort achievable by varying clothing insulation
    2369 !--          Between winter and summer set values
    2370              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,     &
    2371                 pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
     2332          IF ( pmva <= 0.0_wp )  THEN
     2333!
     2334!--          Case: comfort achievable by varying clothing insulation between winter and summer set
     2335!--                values
     2336             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, pmv_s, wclo, pmv_w, eps, &
     2337                              pmva, ncount, clo )
    23722338             IF ( ncount < 0_iwp )  THEN
    23732339                nerr = -1_iwp
     
    23762342          ELSE IF ( pmva > 0.06_wp )  THEN
    23772343             clo = 0.5_wp
    2378              CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
    2379                            pmva )
     2344             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,  pmva )
    23802345          ENDIF
    2381        ELSE IF ( pmva < -0.11_wp )  THEN
     2346       ELSE IF ( pmva < - 0.11_wp )  THEN
    23822347          clo = 1.75_wp
    23832348          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
     
    23902355       pmv_s = pmva
    23912356
    2392        IF ( pmva < 0._wp )  THEN
     2357       IF ( pmva < 0.0_wp )  THEN
    23932358!
    23942359!--       Case winter clothing insulation: cold stress ?
     
    23972362          pmv_w = pmva
    23982363
    2399           IF ( pmva >= 0._wp )  THEN
    2400 !
    2401 !--          Case: comfort achievable by varying clothing insulation
    2402 !--          between winter and summer set values
    2403              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,     &
    2404                                pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
     2364          IF ( pmva >= 0.0_wp )  THEN
     2365!
     2366!--          Case: comfort achievable by varying clothing insulation between winter and summer set
     2367!--                values
     2368             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, pmv_s, wclo, pmv_w, eps, &
     2369                               pmva, ncount, clo )
    24052370             IF ( ncount < 0_iwp )  THEN
    24062371                nerr = -1_iwp
    24072372                RETURN
    24082373             ENDIF
    2409           ELSE IF ( pmva < -0.11_wp )  THEN
     2374          ELSE IF ( pmva < - 0.11_wp )  THEN
    24102375             clo = 1.75_wp
    2411              CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
    2412                            pmva )
     2376             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
    24132377          ENDIF
    24142378       ELSE IF ( pmva > 0.06_wp )  THEN
     
    24232387    CALL perct_regression( pmva, clo, perct_ij )
    24242388    ptc = perct_ij
    2425     IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
     2389    IF ( clo >= 1.75_wp  .AND.  pmva <= - 0.11_wp )  THEN
    24262390!
    24272391!--    Adjust for cold conditions according to Gagge 1986
     
    24292393       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
    24302394       pmvs = pmva - d_pmv
    2431        IF ( pmvs > -0.11_wp )  THEN
    2432           d_pmv  = 0._wp
    2433           pmvs   = -0.11_wp
     2395       IF ( pmvs > - 0.11_wp )  THEN
     2396          d_pmv  = 0.0_wp
     2397          pmvs   = - 0.11_wp
    24342398       ENDIF
    24352399       CALL perct_regression( pmvs, clo, perct_ij )
     
    24392403    IF ( clo > 0.5_wp  .AND.  perct_ij <= 8.73_wp )  THEN
    24402404!
    2441 !--    Required clothing insulation (ireq) is exclusively defined for
    2442 !--    perceived temperatures (perct) less 10 (C) for a
    2443 !--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
     2405!--    Required clothing insulation (ireq) is exclusively defined for perceived temperatures (perct)
     2406!--    less 10 (C) for a reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s.
    24442407       clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
    24452408       clo = clon
     
    24472410    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
    24482411    sultrieness    = .FALSE.
    2449     d_std = -99._wp
     2412    d_std = -99.0_wp
    24502413    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
    24512414!
     
    24552418       pmvs   = pmva + d_pmv
    24562419       CALL perct_regression( pmvs, clo, perct_ij )
    2457        IF ( sult_lim < 99._wp )  THEN
     2420       IF ( sult_lim < 99.0_wp )  THEN
    24582421          IF ( (perct_ij - ptc) > sult_lim )  sultrieness = .TRUE.
    24592422!
     
    24672430 END SUBROUTINE calculate_perct_static
    24682431
    2469 !------------------------------------------------------------------------------!
     2432!--------------------------------------------------------------------------------------------------!
    24702433! Description:
    24712434! ------------
    2472 !> The SUBROUTINE calculates the (saturation) water vapour pressure
    2473 !> (hPa = hecto Pascal) for a given temperature ta (degC).
    2474 !> 'ta' can be the air temperature or the dew point temperature. The first will
    2475 !> result in the current vapor pressure (hPa), the latter will calulate the
    2476 !> saturation vapor pressure (hPa).
    2477 !------------------------------------------------------------------------------!
     2435!> The SUBROUTINE calculates the (saturation) water vapour pressure (hPa = hecto Pascal) for a given
     2436!> temperature ta (degC).
     2437!>'ta' can be the air temperature or the dew point temperature. The first will result in the current
     2438!> vapor pressure (hPa), the latter will calulate the saturation vapor pressure (hPa).
     2439!--------------------------------------------------------------------------------------------------!
    24782440 SUBROUTINE saturation_vapor_pressure( ta, svp_ta )
    24792441
     
    24872449
    24882450
    2489     IF ( ta < 0._wp )  THEN
     2451    IF ( ta < 0.0_wp )  THEN
    24902452!
    24912453!--    ta  < 0 (degC): water vapour pressure over ice
     
    25042466 END SUBROUTINE saturation_vapor_pressure
    25052467
    2506 !------------------------------------------------------------------------------!
     2468!--------------------------------------------------------------------------------------------------!
    25072469! Description:
    25082470! ------------
    2509 !> Find the clothing insulation value clo_res (clo) to make Fanger's Predicted
    2510 !> Mean Vote (PMV) equal comfort (pmva=0) for actual meteorological conditions
    2511 !> (ta,tmrt, vp, ws, pair) and values of individual's activity level
    2512 !------------------------------------------------------------------------------!
    2513  SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo,             &
    2514                        pmv_s, wclo, pmv_w, eps, pmva, nerr,               &
    2515                        clo_res )
     2471!> Find the clothing insulation value clo_res (clo) to make Fanger's Predicted Mean Vote (PMV) equal
     2472!> comfort (pmva=0) for actual meteorological conditions (ta,tmrt, vp, ws, pair) and values of
     2473!> individual's activity level.
     2474!--------------------------------------------------------------------------------------------------!
     2475 SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo, pmv_s, wclo, pmv_w, eps, pmva,  &
     2476                        nerr, clo_res )
    25162477
    25172478    IMPLICIT NONE
    25182479!
    25192480!-- Input variables of argument list:
     2481    REAL(wp), INTENT ( IN )  :: actlev   !< Individuals activity level per unit surface area (W/m2)
     2482    REAL(wp), INTENT ( IN )  :: eps      !< (0.05) accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
     2483    REAL(wp), INTENT ( IN )  :: eta      !< Individuals work efficiency (dimensionless)
     2484    REAL(wp), INTENT ( IN )  :: pair     !< Barometric air pressure (hPa)
     2485    REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
     2486    REAL(wp), INTENT ( IN )  :: pmv_w    !< Fanger's PMV corresponding to wclo
     2487    REAL(wp), INTENT ( IN )  :: sclo     !< Lower threshold of bracketing clothing insulation (clo)
    25202488    REAL(wp), INTENT ( IN )  :: ta       !< Ambient temperature (degC)
    25212489    REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC)
    25222490    REAL(wp), INTENT ( IN )  :: vp       !< Water vapour pressure (hPa)
     2491    REAL(wp), INTENT ( IN )  :: wclo     !< Upper threshold of bracketing clothing insulation (clo)
    25232492    REAL(wp), INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
    2524     REAL(wp), INTENT ( IN )  :: pair     !< Barometric air pressure (hPa)
    2525     REAL(wp), INTENT ( IN )  :: actlev   !< Individuals activity level per unit surface area (W/m2)
    2526     REAL(wp), INTENT ( IN )  :: eta      !< Individuals work efficiency (dimensionless)
    2527     REAL(wp), INTENT ( IN )  :: sclo     !< Lower threshold of bracketing clothing insulation (clo)
    2528     REAL(wp), INTENT ( IN )  :: wclo     !< Upper threshold of bracketing clothing insulation (clo)
    2529     REAL(wp), INTENT ( IN )  :: eps      !< (0.05) accuracy in clothing insulation (clo) for
    2530 !                                          evaluation the root of Fanger's PMV (pmva=0)
    2531     REAL(wp), INTENT ( IN )  :: pmv_w    !< Fanger's PMV corresponding to wclo
    2532     REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
    25332493!
    25342494!-- Output variables of argument list:
    2535     REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
    2536     REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
    25372495    INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
    25382496                                         !< nerr >= 0, o.k., and nerr is the number of iterations for convergence
    25392497                                         !< nerr = -1: error = malfunction of Ridder's convergence method
    2540                                          !< nerr = -2: error = maximum iterations (max_iteration) exceeded 
     2498                                         !< nerr = -2: error = maximum iterations (max_iteration) exceeded
    25412499                                         !< nerr = -3: error = root not bracketed between sclo and wclo
     2500
     2501    REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
     2502    REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
    25422503!
    25432504!-- Type of program variables
    25442505    INTEGER(iwp), PARAMETER  ::  max_iteration = 15_iwp       !< max number of iterations
     2506    INTEGER(iwp) ::  j       !< running index
     2507
    25452508    REAL(wp),     PARAMETER  ::  guess_0       = -1.11e30_wp  !< initial guess
    2546     REAL(wp) ::  x_ridder    !< current guess for clothing insulation   (clo)
     2509
    25472510    REAL(wp) ::  clo_lower   !< lower limit of clothing insulation      (clo)
    25482511    REAL(wp) ::  clo_upper   !< upper limit of clothing insulation      (clo)
     2512    REAL(wp) ::  sroot       !< sqrt of PMV-guess
     2513    REAL(wp) ::  x_average   !< average of x_lower and x_upper          (clo)
    25492514    REAL(wp) ::  x_lower     !< lower guess for clothing insulation     (clo)
     2515    REAL(wp) ::  x_new       !< preliminary result for clothing insulation (clo)
     2516    REAL(wp) ::  x_ridder    !< current guess for clothing insulation   (clo)
    25502517    REAL(wp) ::  x_upper     !< upper guess for clothing insulation     (clo)
    2551     REAL(wp) ::  x_average   !< average of x_lower and x_upper          (clo)
    2552     REAL(wp) ::  x_new       !< preliminary result for clothing insulation (clo)
     2518    REAL(wp) ::  y_average   !< average of y_lower and y_upper
     2519    REAL(wp) ::  y_new       !< preliminary result for pred. mean vote
    25532520    REAL(wp) ::  y_lower     !< predicted mean vote for summer clothing
    25542521    REAL(wp) ::  y_upper     !< predicted mean vote for winter clothing
    2555     REAL(wp) ::  y_average   !< average of y_lower and y_upper
    2556     REAL(wp) ::  y_new       !< preliminary result for pred. mean vote
    2557     REAL(wp) ::  sroot       !< sqrt of PMV-guess
    2558     INTEGER(iwp) ::  j       !< running index
    25592522!
    25602523!-- Initialise
     
    25632526!-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation
    25642527    x_ridder    = bio_fill_value
    2565     pmva        = 0._wp
     2528    pmva        = 0.0_wp
    25662529    clo_lower   = sclo
    25672530    y_lower     = pmv_s
    25682531    clo_upper   = wclo
    25692532    y_upper     = pmv_w
    2570     IF ( ( y_lower > 0._wp  .AND.  y_upper < 0._wp )  .OR.                     &
    2571          ( y_lower < 0._wp  .AND.  y_upper > 0._wp ) )  THEN
     2533    IF ( ( y_lower > 0.0_wp .AND. y_upper < 0.0_wp )  .OR.                                         &
     2534         ( y_lower < 0.0_wp .AND. y_upper > 0.0_wp ) )  THEN
    25722535       x_lower  = clo_lower
    25732536       x_upper  = clo_upper
     
    25762539       DO  j = 1_iwp, max_iteration
    25772540          x_average = 0.5_wp * ( x_lower + x_upper )
    2578           CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta,        &
    2579                         y_average )
     2541          CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta, y_average )
    25802542          sroot = SQRT( y_average**2 - y_lower * y_upper )
    25812543          IF ( ABS( sroot ) < 0.00001_wp )  THEN
     
    25842546             RETURN
    25852547          ENDIF
    2586           x_new = x_average + ( x_average - x_lower ) *                        &
    2587                       ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot )
     2548          x_new = x_average + ( x_average - x_lower ) *                                            &
     2549                  ( SIGN ( 1.0_wp, y_lower - y_upper ) * y_average / sroot )
    25882550          IF ( ABS( x_new - x_ridder ) <= eps )  THEN
    25892551             clo_res = x_ridder
     
    25922554          ENDIF
    25932555          x_ridder = x_new
    2594           CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
    2595                         y_new )
     2556          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta, y_new )
    25962557          IF ( ABS( y_new ) < 0.00001_wp )  THEN
    25972558             clo_res = x_ridder
     
    26122573          ELSE
    26132574!
    2614 !--          Never get here in x_ridder: singularity in y
     2575!--          Never get here in x_ridder: SINgularity in y
    26152576             nerr    = -1_iwp
    26162577             clo_res = x_ridder
     
    26422603 END SUBROUTINE iso_ridder
    26432604
    2644 !------------------------------------------------------------------------------!
     2605!--------------------------------------------------------------------------------------------------!
    26452606! Description:
    26462607! ------------
    2647 !> Regression relations between perceived temperature (perct) and (adjusted)
    2648 !> PMV. The regression presumes the Klima-Michel settings for reference
    2649 !> individual and reference environment.
    2650 !------------------------------------------------------------------------------!
     2608!> Regression relations between perceived temperature (perct) and (adjusted) PMV. The regression
     2609!> presumes the Klima-Michel settings for reference individual and reference environment.
     2610!--------------------------------------------------------------------------------------------------!
    26512611 SUBROUTINE perct_regression( pmv, clo, perct_ij )
    26522612
    26532613    IMPLICIT NONE
    26542614
     2615    REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
    26552616    REAL(wp), INTENT ( IN ) ::  pmv   !< Fangers predicted mean vote (dimensionless)
    2656     REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
    26572617
    26582618    REAL(wp), INTENT ( OUT ) ::  perct_ij   !< perct (degC) corresponding to given PMV / clo
    26592619
    2660     IF ( pmv <= -0.11_wp )  THEN
     2620    IF ( pmv <= - 0.11_wp )  THEN
    26612621       perct_ij = 5.805_wp + 12.6784_wp * pmv
    26622622    ELSE
     
    26702630 END SUBROUTINE perct_regression
    26712631
    2672 !------------------------------------------------------------------------------!
     2632!--------------------------------------------------------------------------------------------------!
    26732633! Description:
    26742634! ------------
    26752635!> FANGER.F90
    26762636!>
    2677 !> SI-VERSION: ACTLEV W m-2, DAMPFDRUCK hPa
    2678 !> Berechnet das aktuelle Predicted Mean Vote nach Fanger
    2679 !>
     2637!> SI-VERSION: ACTLEV W m-2, VAPOUR PRESSURE hPa
     2638!> Calculates the current Predicted Mean Vote according to Fanger.
    26802639!> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s
    2681 !------------------------------------------------------------------------------!
     2640!--------------------------------------------------------------------------------------------------!
    26822641 SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva )
    26832642
     
    26852644!
    26862645!-- Input variables of argument list:
     2646    REAL(wp), INTENT ( IN ) ::  actlev   !< Individuals activity level per unit surface area (W/m2)
     2647    REAL(wp), INTENT ( IN ) ::  eta      !< Individuals mechanical work efficiency (dimensionless)
     2648    REAL(wp), INTENT ( IN ) ::  in_clo   !< Clothing insulation (clo)
     2649    REAL(wp), INTENT ( IN ) ::  in_ws    !< Wind speed (m/s) 1 m above ground
     2650    REAL(wp), INTENT ( IN ) ::  pa       !< Water vapour pressure (hPa)
     2651    REAL(wp), INTENT ( IN ) ::  pair     !< Barometric pressure (hPa) at site
    26872652    REAL(wp), INTENT ( IN ) ::  ta       !< Ambient air temperature (degC)
    26882653    REAL(wp), INTENT ( IN ) ::  tmrt     !< Mean radiant temperature (degC)
    2689     REAL(wp), INTENT ( IN ) ::  pa       !< Water vapour pressure (hPa)
    2690     REAL(wp), INTENT ( IN ) ::  pair     !< Barometric pressure (hPa) at site
    2691     REAL(wp), INTENT ( IN ) ::  in_ws    !< Wind speed (m/s) 1 m above ground
    2692     REAL(wp), INTENT ( IN ) ::  in_clo   !< Clothing insulation (clo)
    2693     REAL(wp), INTENT ( IN ) ::  actlev   !< Individuals activity level per unit surface area (W/m2)
    2694     REAL(wp), INTENT ( IN ) ::  eta      !< Individuals mechanical work efficiency (dimensionless)
     2654
    26952655!
    26962656!-- Output variables of argument list:
    2697     REAL(wp), INTENT ( OUT ) ::  pmva    !< Actual Predicted Mean Vote (PMV, 
    2698                                          !< dimensionless) according to Fanger corresponding to meteorological 
     2657    REAL(wp), INTENT ( OUT ) ::  pmva    !< Actual Predicted Mean Vote (PMV,
     2658                                         !< dimensionless) according to Fanger corresponding to meteorological
    26992659                                         !< (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
    27002660!
    27012661!-- Internal variables
    2702     REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
    2703     REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
     2662    INTEGER(iwp) :: i         !< running index
     2663
    27042664    REAL(wp) ::  activity     !< persons activity  (must stay == actlev, W)
    2705     REAL(wp) ::  t_skin_aver  !< average skin temperature               (degree_C)
    27062665    REAL(wp) ::  bc           !< preliminary result storage
    27072666    REAL(wp) ::  cc           !< preliminary result storage
     2667    REAL(wp) ::  clo          !< clothing insulation index              (clo)
    27082668    REAL(wp) ::  dc           !< preliminary result storage
    27092669    REAL(wp) ::  ec           !< preliminary result storage
     2670    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
    27102671    REAL(wp) ::  gc           !< preliminary result storage
     2672    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
     2673    REAL(wp) ::  hr           !< radiational heat resistence
    27112674    REAL(wp) ::  t_clothing   !< clothing temperature                   (degree_C)
    2712     REAL(wp) ::  hr           !< radiational heat resistence
    2713     REAL(wp) ::  clo          !< clothing insulation index              (clo)
     2675    REAL(wp) ::  t_skin_aver  !< average skin temperature               (degree_C)
    27142676    REAL(wp) ::  ws           !< wind speed                             (m/s)
    2715     REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat 
     2677    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat
    27162678                              !< ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
    27172679    REAL(wp) ::  z2           !< Water vapour diffution through the skin
     
    27202682    REAL(wp) ::  z5           !< Loss of radiational heat
    27212683    REAL(wp) ::  z6           !< Heat loss through forced convection
    2722     INTEGER(iwp) :: i         !< running index
     2684
    27232685!
    27242686!-- Clo must be > 0. to avoid div. by 0!
    27252687    clo = in_clo
    2726     IF ( clo <= 0._wp )  clo = .001_wp
    2727 !
    2728 !-- f_cl = Increase in surface due to clothing
    2729     f_cl = 1._wp + .15_wp * clo
     2688    IF ( clo <= 0.0_wp )  clo = .001_wp
     2689!
     2690!-- f_cl = increase in surface due to clothing
     2691    f_cl = 1.0_wp + 0.15_wp * clo
    27302692!
    27312693!-- Case of free convection (ws < 0.1 m/s ) not considered
    27322694    ws = in_ws
    2733     IF ( ws < .1_wp )  THEN
    2734        ws = .1_wp
     2695    IF ( ws < 0.1_wp )  THEN
     2696       ws = 0.1_wp
    27352697    ENDIF
    27362698!
     
    27382700    heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp )
    27392701!
    2740 !-- Activity = inner heat produktion per standardized surface
    2741     activity = actlev * ( 1._wp - eta )
    2742 !
    2743 !-- T_skin_aver = average skin temperature
    2744     t_skin_aver = 35.7_wp - .0275_wp * activity
     2702!-- Activity = inner heat production per standardized surface
     2703    activity = actlev * ( 1.0_wp - eta )
     2704!
     2705!-- t_skin_aver = average skin temperature
     2706    t_skin_aver = 35.7_wp - 0.0275_wp * activity
    27452707!
    27462708!-- Calculation of constants for evaluation below
    2747     bc = .155_wp * clo * 3.96_wp * 10._wp**( -8 ) * f_cl
     2709    bc = 0.155_wp * clo * 3.96_wp * 10.0_wp**( -8 ) * f_cl
    27482710    cc = f_cl * heat_convection
    2749     ec = .155_wp * clo
    2750     dc = ( 1._wp + ec * cc ) / bc
     2711    ec = 0.155_wp * clo
     2712    dc = ( 1.0_wp + ec * cc ) / bc
    27512713    gc = ( t_skin_aver + bc * ( tmrt + degc_to_k )**4 + ec * cc * ta ) / bc
    27522714!
    2753 !-- Calculation of clothing surface temperature (t_clothing) based on
    2754 !-- Newton-approximation with air temperature as initial guess
     2715!-- Calculation of clothing surface temperature (t_clothing) based on Newton-approximation with air
     2716!-- temperature as initial guess.
    27552717    t_clothing = ta
    27562718    DO  i = 1, 3
    2757        t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing &
    2758           * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc )
     2719       t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing * dc - gc ) /        &
     2720                    ( 4.0_wp * ( t_clothing + degc_to_k )**3 + dc )
    27592721    ENDDO
    27602722!
    2761 !-- Empiric factor for the adaption of the heat ballance equation
    2762 !-- to the psycho-physical scale (Equ. 40 in FANGER)
    2763     z1 = ( .303_wp * EXP( -.036_wp * actlev ) + .0275_wp )
     2723!-- Empiric factor for the adaption of the heat ballance equation to the psycho-physical scale (Equ.
     2724!-- 40 in FANGER)
     2725    z1 = ( 0.303_wp * EXP( - 0.036_wp * actlev ) + 0.0275_wp )
    27642726!
    27652727!-- Water vapour diffution through the skin
    2766     z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
     2728    z2 = 0.31_wp * ( 57.3_wp - 0.07_wp * activity-pa )
    27672729!
    27682730!-- Sweat evaporation from the skin surface
    2769     z3 = .42_wp * ( activity - 58._wp )
     2731    z3 = 0.42_wp * ( activity - 58.0_wp )
    27702732!
    27712733!-- Loss of latent heat through respiration
    2772     z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
    2773       ( 34._wp - ta )
     2734    z4 = 0.0017_wp * actlev * ( 58.7_wp - pa ) + 0.0014_wp * actlev *                              &
     2735         ( 34.0_wp - ta )
    27742736!
    27752737!-- Loss of radiational heat
    2776     z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
    2777        degc_to_k )**4 )
    2778     IF ( ABS( t_clothing - tmrt ) > 0._wp )  THEN
     2738    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt + degc_to_k )**4 )
     2739    IF ( ABS( t_clothing - tmrt ) > 0.0_wp )  THEN
    27792740       hr = z5 / f_cl / ( t_clothing - tmrt )
    27802741    ELSE
    2781        hr = 0._wp
     2742       hr = 0.0_wp
    27822743    ENDIF
    27832744!
     
    27902751 END SUBROUTINE fanger
    27912752
    2792 !------------------------------------------------------------------------------!
     2753!--------------------------------------------------------------------------------------------------!
    27932754! Description:
    27942755! ------------
    2795 !> For pmva > 0 and clo =0.5 the increment (deltapmv) is calculated
    2796 !> that converts pmva into Gagge's et al. (1986) PMV*.
    2797 !------------------------------------------------------------------------------!
     2756!> For pmva > 0 and clo =0.5 the increment (deltapmv) is calculated that converts pmva into Gagge's
     2757!> et al. (1986) PMV*.
     2758!--------------------------------------------------------------------------------------------------!
    27982759 REAL(wp) FUNCTION deltapmv( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
    27992760
     
    28032764!-- Input variables of argument list:
    28042765    REAL(wp),     INTENT ( IN )  :: pmva     !< Actual Predicted Mean Vote (PMV) according to Fanger
     2766    REAL(wp),     INTENT ( IN )  :: svp_ta   !< Saturation water vapour pressure (hPa) at ta
    28052767    REAL(wp),     INTENT ( IN )  :: ta       !< Ambient temperature (degC) at screen level
     2768    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
    28062769    REAL(wp),     INTENT ( IN )  :: vp       !< Water vapour pressure (hPa) at screen level
    2807     REAL(wp),     INTENT ( IN )  :: svp_ta   !< Saturation water vapour pressure (hPa) at ta
    2808     REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
    28092770    REAL(wp),     INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
    28102771
     
    28192780!
    28202781!-- Internal variables:
    2821     REAL(wp) ::  pmv          !< temp storage og predicted mean vote
    2822     REAL(wp) ::  pa_p50       !< ratio actual water vapour pressure to that of relative humidity of 50 %
    2823     REAL(wp) ::  pa           !< vapor pressure (hPa) with hard bounds
     2782    INTEGER(iwp) :: nreg      !<
     2783
    28242784    REAL(wp) ::  apa          !< natural logarithm of pa (with hard lower border)
    28252785    REAL(wp) ::  dapa         !< difference of apa and pa_p50
    2826     REAL(wp) ::  sqvel        !< square root of local wind velocity
     2786    REAL(wp) ::  dpmv_1       !<
     2787    REAL(wp) ::  dpmv_2       !<
    28272788    REAL(wp) ::  dtmrt        !< difference mean radiation to air temperature
     2789    REAL(wp) ::  pa           !< vapor pressure (hPa) with hard bounds
     2790    REAL(wp) ::  pa_p50       !< ratio actual water vapour pressure to that of relative humidity of
     2791                              !< 50 %
     2792    REAL(wp) ::  pmv          !< temp storage og predicted mean vote
     2793    REAL(wp) ::  pmvs         !<
    28282794    REAL(wp) ::  p10          !< lower bound for pa
    28292795    REAL(wp) ::  p95          !< upper bound for pa
    2830     REAL(wp) ::  weight       !<
    2831     REAL(wp) ::  weight2      !<
    2832     REAL(wp) ::  dpmv_1       !<
    2833     REAL(wp) ::  dpmv_2       !<
    2834     REAL(wp) ::  pmvs         !<
    2835     INTEGER(iwp) :: nreg      !<
     2796    REAL(wp) ::  sqvel        !< square root of local wind velocity
     2797    REAL(wp) ::  weight       !<
     2798    REAL(wp) ::  weight2      !<
    28362799
    28372800!
    28382801!-- Regression coefficients:
    2839     REAL(wp), DIMENSION(0:7), PARAMETER ::  bpmv = (/                          &
    2840      -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp,               &
    2841      -0.3551048_wp, -0.4304076_wp, -0.4884961_wp, -0.4897495_wp /)
    2842 
    2843     REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa_p50 = (/                       &
    2844      -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp,               &
    2845       0.4285938_wp,  0.6281256_wp,  0.5067361_wp,  0.3965169_wp /)
    2846 
    2847     REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa = (/                           &
    2848       0.0580284_wp,  0.0836264_wp,  0.1009919_wp,  0.1020777_wp,               &
    2849       0.0898681_wp,  0.0839116_wp,  0.0853258_wp,  0.0866589_wp /)
    2850 
    2851     REAL(wp), DIMENSION(0:7), PARAMETER ::  bapa = (/                          &
    2852      -1.7838788_wp, -2.9306231_wp, -1.6350334_wp,   0.6211547_wp,              &
    2853       3.3918083_wp,  5.5521025_wp,  8.4897418_wp,  16.6265851_wp /)
    2854 
    2855     REAL(wp), DIMENSION(0:7), PARAMETER ::  bdapa = (/                         &
    2856       1.6752720_wp,  2.7379504_wp,  1.2940526_wp,  -1.0985759_wp,              &
    2857      -3.9054732_wp, -6.0403012_wp, -8.9437119_wp, -17.0671201_wp /)
    2858 
    2859     REAL(wp), DIMENSION(0:7), PARAMETER ::  bsqvel = (/                        &
    2860      -0.0315598_wp, -0.0286272_wp, -0.0009228_wp,  0.0483344_wp,               &
    2861       0.0992366_wp,  0.1491379_wp,  0.1951452_wp,  0.2133949_wp /)
    2862 
    2863     REAL(wp), DIMENSION(0:7), PARAMETER ::  bta = (/                           &
    2864       0.0953986_wp,  0.1524760_wp,  0.0564241_wp, -0.0893253_wp,               &
    2865      -0.2398868_wp, -0.3515237_wp, -0.5095144_wp, -0.9469258_wp /)
    2866 
    2867     REAL(wp), DIMENSION(0:7), PARAMETER ::  bdtmrt = (/                        &
    2868      -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp,               &
    2869      -0.0069036_wp, -0.0075844_wp, -0.0079602_wp, -0.0089439_wp /)
    2870 
    2871     REAL(wp), DIMENSION(0:7), PARAMETER ::  aconst = (/                        &
    2872       1.8686215_wp,  3.4260713_wp,   2.0116185_wp,  -0.7777552_wp,             &
    2873      -4.6715853_wp, -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /)
     2802    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpmv = (/                                              &
     2803     - 0.0556602_wp, - 0.1528680_wp, - 0.2336104_wp, - 0.2789387_wp,                               &
     2804     - 0.3551048_wp, - 0.4304076_wp, - 0.4884961_wp, - 0.4897495_wp /)
     2805
     2806    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa_p50 = (/                                           &
     2807     - 0.1607154_wp, - 0.4177296_wp, - 0.4120541_wp, - 0.0886564_wp,                               &
     2808       0.4285938_wp,   0.6281256_wp,   0.5067361_wp,   0.3965169_wp /)
     2809
     2810    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa = (/                                               &
     2811       0.0580284_wp,   0.0836264_wp,   0.1009919_wp,   0.1020777_wp,                               &
     2812       0.0898681_wp,   0.0839116_wp,   0.0853258_wp,   0.0866589_wp /)
     2813
     2814    REAL(wp), DIMENSION(0:7), PARAMETER ::  bapa = (/                                              &
     2815     - 1.7838788_wp, - 2.9306231_wp, - 1.6350334_wp,    0.6211547_wp,                              &
     2816       3.3918083_wp,   5.5521025_wp,   8.4897418_wp,   16.6265851_wp /)
     2817
     2818    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdapa = (/                                             &
     2819       1.6752720_wp,   2.7379504_wp,   1.2940526_wp, -  1.0985759_wp,                              &
     2820     - 3.9054732_wp, - 6.0403012_wp, - 8.9437119_wp, - 17.0671201_wp /)
     2821
     2822    REAL(wp), DIMENSION(0:7), PARAMETER ::  bsqvel = (/                                            &
     2823     - 0.0315598_wp, - 0.0286272_wp, - 0.0009228_wp,   0.0483344_wp,                               &
     2824       0.0992366_wp,   0.1491379_wp,   0.1951452_wp,   0.2133949_wp /)
     2825
     2826    REAL(wp), DIMENSION(0:7), PARAMETER ::  bta = (/                                               &
     2827       0.0953986_wp,   0.1524760_wp,   0.0564241_wp, - 0.0893253_wp,                               &
     2828     - 0.2398868_wp, - 0.3515237_wp, - 0.5095144_wp, - 0.9469258_wp /)
     2829
     2830    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdtmrt = (/                                            &
     2831     - 0.0004672_wp, - 0.0000514_wp, - 0.0018037_wp, - 0.0049440_wp,                               &
     2832     - 0.0069036_wp, - 0.0075844_wp, - 0.0079602_wp, - 0.0089439_wp /)
     2833
     2834    REAL(wp), DIMENSION(0:7), PARAMETER ::  aconst = (/                                            &
     2835       1.8686215_wp,   3.4260713_wp,    2.0116185_wp, -  0.7777552_wp,                             &
     2836     - 4.6715853_wp, - 7.7314281_wp, - 11.7602578_wp, - 23.5934198_wp /)
    28742837
    28752838
     
    28852848    pmv  = pmva
    28862849!
    2887 !-- Water vapour pressure of air 
     2850!-- Water vapour pressure of air
    28882851    p10  = 0.05_wp * svp_ta
    28892852    p95  = 1.00_wp * svp_ta
     
    29022865       ENDIF
    29032866    ENDIF
    2904     IF ( pa > 0._wp )  THEN
     2867    IF ( pa > 0.0_wp )  THEN
    29052868!
    29062869!--    Natural logarithm of pa
    29072870       apa = LOG( pa )
    29082871    ELSE
    2909        apa = -5._wp
     2872       apa = -5.0_wp
    29102873    ENDIF
    29112874!
    29122875!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
    29132876    pa_p50   = 0.5_wp * svp_ta
    2914     IF ( pa_p50 > 0._wp  .AND.  pa > 0._wp )  THEN
     2877    IF ( pa_p50 > 0.0_wp  .AND.  pa > 0.0_wp )  THEN
    29152878       dapa   = apa - LOG( pa_p50 )
    29162879       pa_p50 = pa / pa_p50
    29172880    ELSE
    2918        dapa   = -5._wp
    2919        pa_p50 = 0._wp
     2881       dapa   = -5.0_wp
     2882       pa_p50 = 0.0_wp
    29202883    ENDIF
    29212884!
    29222885!-- Square root of wind velocity
    2923     IF ( ws >= 0._wp )  THEN
     2886    IF ( ws >= 0.0_wp )  THEN
    29242887       sqvel = SQRT( ws )
    29252888    ELSE
    2926        sqvel = 0._wp
     2889       sqvel = 0.0_wp
    29272890    ENDIF
    29282891!
     
    29342897    IF ( nreg < 0_iwp )  THEN
    29352898!
    2936 !--    value of the FUNCTION in the case pmv <= -1
    2937        deltapmv = 0._wp
     2899!--    Value of the FUNCTION in the case pmv <= -1
     2900       deltapmv = 0.0_wp
    29382901       RETURN
    29392902    ENDIF
    2940     weight = MOD ( pmv, 1._wp )
    2941     IF ( weight < 0._wp )  weight = 0._wp
     2903    weight = MOD ( pmv, 1.0_wp )
     2904    IF ( weight < 0.0_wp )  weight = 0.0_wp
    29422905    IF ( nreg > 5_iwp )  THEN
    29432906       nreg  = 5_iwp
    2944        weight   = pmv - 5._wp
    2945        weight2  = pmv - 6._wp
     2907       weight   = pmv - 5.0_wp
     2908       weight2  = pmv - 6.0_wp
    29462909       IF ( weight2 > 0_iwp )  THEN
    29472910          weight = ( weight - weight2 ) / weight
     
    29502913!
    29512914!-- Regression valid for 0. <= pmv <= 6., bounds are checked above
    2952     dpmv_1 =                                                                   &
    2953        + bpa(nreg) * pa                                                        &
    2954        + bpmv(nreg) * pmv                                                      &
    2955        + bapa(nreg) * apa                                                      &
    2956        + bta(nreg) * ta                                                        &
    2957        + bdtmrt(nreg) * dtmrt                                                  &
    2958        + bdapa(nreg) * dapa                                                    &
    2959        + bsqvel(nreg) * sqvel                                                  &
    2960        + bpa_p50(nreg) * pa_p50                                                &
    2961        + aconst(nreg)
    2962 
    2963 !    dpmv_2 = 0._wp
     2915    dpmv_1 =                                                                                       &
     2916             + bpa(nreg)     * pa                                                                  &
     2917             + bpmv(nreg)    * pmv                                                                 &
     2918             + bapa(nreg)    * apa                                                                 &
     2919             + bta(nreg)     * ta                                                                  &
     2920             + bdtmrt(nreg)  * dtmrt                                                               &
     2921             + bdapa(nreg)   * dapa                                                                &
     2922             + bsqvel(nreg)  * sqvel                                                               &
     2923             + bpa_p50(nreg) * pa_p50                                                              &
     2924             + aconst(nreg)
     2925
     2926!    dpmv_2 = 0.0_wp
    29642927!    IF ( nreg < 6_iwp )  THEN  !< nreg is always <= 5, see above
    2965     dpmv_2 =                                                                   &
    2966        + bpa(nreg+1_iwp)     * pa                                              &
    2967        + bpmv(nreg+1_iwp)    * pmv                                             &
    2968        + bapa(nreg+1_iwp)    * apa                                             &
    2969        + bta(nreg+1_iwp)     * ta                                              &
    2970        + bdtmrt(nreg+1_iwp)  * dtmrt                                           &
    2971        + bdapa(nreg+1_iwp)   * dapa                                            &
    2972        + bsqvel(nreg+1_iwp)  * sqvel                                           &
    2973        + bpa_p50(nreg+1_iwp) * pa_p50                                          &
    2974        + aconst(nreg+1_iwp)
     2928    dpmv_2 =                                                                                       &
     2929             + bpa(nreg+1_iwp)     * pa                                                            &
     2930             + bpmv(nreg+1_iwp)    * pmv                                                           &
     2931             + bapa(nreg+1_iwp)    * apa                                                           &
     2932             + bta(nreg+1_iwp)     * ta                                                            &
     2933             + bdtmrt(nreg+1_iwp)  * dtmrt                                                         &
     2934             + bdapa(nreg+1_iwp)   * dapa                                                          &
     2935             + bsqvel(nreg+1_iwp)  * sqvel                                                         &
     2936             + bpa_p50(nreg+1_iwp) * pa_p50                                                        &
     2937             + aconst(nreg+1_iwp)
    29752938!    ENDIF
    29762939!
    29772940!-- Calculate pmv modification
    2978     deltapmv = ( 1._wp - weight ) * dpmv_1 + weight * dpmv_2
     2941    deltapmv = ( 1.0_wp - weight ) * dpmv_1 + weight * dpmv_2
    29792942    pmvs = pmva + deltapmv
    2980     IF ( ( pmvs ) < 0._wp )  THEN
     2943    IF ( ( pmvs ) < 0.0_wp )  THEN
    29812944!
    29822945!--    Prevent negative pmv* due to problems with clothing insulation
     
    29892952!
    29902953!--       Set pmvs to "0" for compliance with summer clothing insulation
    2991           deltapmv = -1._wp * pmva
     2954          deltapmv = -1.0_wp * pmva
    29922955       ENDIF
    29932956    ENDIF
     
    29952958 END FUNCTION deltapmv
    29962959
    2997 !------------------------------------------------------------------------------!
     2960!--------------------------------------------------------------------------------------------------!
    29982961! Description:
    29992962! ------------
    3000 !> The subroutine "calc_sultr" returns a threshold value to perceived
    3001 !> temperature allowing to decide whether the actual perceived temperature
    3002 !> is linked to perecption of sultriness. The threshold values depends
    3003 !> on the Fanger's classical PMV, expressed here as perceived temperature
    3004 !> perct.
    3005 !------------------------------------------------------------------------------!
     2963!> The subroutine "calc_sultr" returns a threshold value to perceived temperature allowing to decide
     2964!> whether the actual perceived temperature is linked to perecption of sultriness. The threshold
     2965!> values depends on the Fanger's classical PMV, expressed here as perceived temperature perct.
     2966!--------------------------------------------------------------------------------------------------!
    30062967 SUBROUTINE calc_sultr( perct_ij, dperctm, dperctstd, sultr_res )
    30072968
     
    30092970!
    30102971!-- Input of the argument list:
    3011     REAL(wp), INTENT ( IN )  ::  perct_ij      !< Classical perceived temperature: Base is Fanger's PMV
     2972    REAL(wp), INTENT ( IN )  ::  perct_ij   !< Classical perceived temperature: Base is Fanger's PMV
    30122973!
    30132974!-- Additional output variables of argument list:
    3014     REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational gt
    3015                                             !< calculated based on Gagge's rational PMV*)
    3016     REAL(wp), INTENT ( OUT ) ::  dperctstd  !< dperctm plus its standard deviation times a factor 
     2975    REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational
     2976                                            !< gt calculated based on Gagge's rational PMV*)
     2977    REAL(wp), INTENT ( OUT ) ::  dperctstd  !< dperctm plus its standard deviation times a factor
    30172978                                            !< determining the significance to perceive sultriness
    30182979    REAL(wp), INTENT ( OUT ) ::  sultr_res
    30192980!
    30202981!-- Types of coefficients mean deviation: third order polynomial
    3021     REAL(wp), PARAMETER ::  dperctka =  7.5776086_wp
    3022     REAL(wp), PARAMETER ::  dperctkb = -0.740603_wp
    3023     REAL(wp), PARAMETER ::  dperctkc =  0.0213324_wp
    3024     REAL(wp), PARAMETER ::  dperctkd = -0.00027797237_wp
     2982    REAL(wp), PARAMETER ::  dperctka =   7.5776086_wp
     2983    REAL(wp), PARAMETER ::  dperctkb = - 0.740603_wp
     2984    REAL(wp), PARAMETER ::  dperctkc =   0.0213324_wp
     2985    REAL(wp), PARAMETER ::  dperctkd = - 0.00027797237_wp
    30252986!
    30262987!-- Types of coefficients mean deviation plus standard deviation
    30272988!-- regression coefficients: third order polynomial
    3028     REAL(wp), PARAMETER ::  dperctsa =  0.0268918_wp
    3029     REAL(wp), PARAMETER ::  dperctsb =  0.0465957_wp
    3030     REAL(wp), PARAMETER ::  dperctsc = -0.00054709752_wp
    3031     REAL(wp), PARAMETER ::  dperctsd =  0.0000063714823_wp
    3032 !
    3033 !-- Factor to mean standard deviation defining SIGNificance for 
     2989    REAL(wp), PARAMETER ::  dperctsa =   0.0268918_wp
     2990    REAL(wp), PARAMETER ::  dperctsb =   0.0465957_wp
     2991    REAL(wp), PARAMETER ::  dperctsc = - 0.00054709752_wp
     2992    REAL(wp), PARAMETER ::  dperctsd =   0.0000063714823_wp
     2993!
     2994!-- Factor to mean standard deviation defining SIGNificance for
    30342995!-- sultriness
    3035     REAL(wp), PARAMETER :: faktor = 1._wp
     2996    REAL(wp), PARAMETER :: faktor = 1.0_wp
    30362997!
    30372998!-- Initialise
    3038     sultr_res = 99._wp
    3039     dperctm   = 0._wp
    3040     dperctstd = 999999._wp
    3041 
    3042     IF ( perct_ij < 16.826_wp  .OR.  perct_ij > 56._wp )  THEN
     2999    sultr_res = 99.0_wp
     3000    dperctm   = 0.0_wp
     3001    dperctstd = 999999.0_wp
     3002
     3003    IF ( perct_ij < 16.826_wp  .OR.  perct_ij > 56.0_wp )  THEN
    30433004!
    30443005!--    Unallowed value of classical perct!
     
    30473008!
    30483009!-- Mean deviation dependent on perct
    3049     dperctm = dperctka + dperctkb * perct_ij + dperctkc * perct_ij**2._wp +    &
    3050        dperctkd * perct_ij**3._wp
     3010    dperctm = dperctka + dperctkb * perct_ij + dperctkc * perct_ij**2.0_wp + dperctkd *            &
     3011              perct_ij**3.0_wp
    30513012!
    30523013!-- Mean deviation plus its standard deviation
    3053     dperctstd = dperctsa + dperctsb * perct_ij + dperctsc * perct_ij**2._wp +  &
    3054        dperctsd * perct_ij**3._wp
     3014    dperctstd = dperctsa + dperctsb * perct_ij + dperctsc * perct_ij**2.0_wp + dperctsd *          &
     3015                perct_ij**3.0_wp
    30553016!
    30563017!-- Value of the FUNCTION
    30573018    sultr_res = dperctm + faktor * dperctstd
    3058     IF ( ABS( sultr_res ) > 99._wp )  sultr_res = +99._wp
     3019    IF ( ABS( sultr_res ) > 99.0_wp )  sultr_res = +99.0_wp
    30593020
    30603021 END SUBROUTINE calc_sultr
    30613022
    3062 !------------------------------------------------------------------------------!
     3023!--------------------------------------------------------------------------------------------------!
    30633024! Description:
    30643025! ------------
    3065 !> Multiple linear regression to calculate an increment delta_cold,
    3066 !> to adjust Fanger's classical PMV (pmva) by Gagge's 2 node model,
    3067 !> applying Fanger's convective heat transfer coefficient, hcf.
    3068 !> Wind velocitiy of the reference environment is 0.10 m/s
    3069 !------------------------------------------------------------------------------!
     3026!> Multiple linear regression to calculate an increment delta_cold, to adjust Fanger's classical PMV
     3027!> (pmva) by Gagge's 2 node model, applying Fanger's convective heat transfer coefficient, hcf.
     3028!> Wind velocitiy of the reference environment is 0.10 m/s
     3029!--------------------------------------------------------------------------------------------------!
    30703030 SUBROUTINE dpmv_cold( pmva, ta, ws, tmrt, nerr, dpmv_cold_res )
    30713031
     
    30753035    REAL(wp), INTENT ( IN ) ::  pmva   !< Fanger's classical predicted mean vote
    30763036    REAL(wp), INTENT ( IN ) ::  ta     !< Air temperature 2 m above ground (degC)
     3037    REAL(wp), INTENT ( IN ) ::  tmrt   !< Mean radiant temperature (degC)
    30773038    REAL(wp), INTENT ( IN ) ::  ws     !< Relative wind velocity 1 m above ground (m/s)
    3078     REAL(wp), INTENT ( IN ) ::  tmrt   !< Mean radiant temperature (degC)
    30793039!
    30803040!-- Type of output argument
    3081     INTEGER(iwp), INTENT ( OUT ) ::  nerr !< Error indicator: 0 = o.k., +1 = denominator for intersection = 0
    3082     REAL(wp),     INTENT ( OUT ) ::  dpmv_cold_res    !< Increment to adjust pmva according to the results of Gagge's
    3083                                                       !< 2 node model depending on the input
     3041    INTEGER(iwp), INTENT ( OUT ) ::  nerr !< Error indicator: 0 = o.k., +1 = denominator for
     3042                                          !< intersection = 0
     3043
     3044    REAL(wp),     INTENT ( OUT ) ::  dpmv_cold_res    !< Increment to adjust pmva according to the
     3045                                                      !< results of Gagge's 2 node model depending on the input
    30843046!
    30853047!-- Type of program variables
     3048    INTEGER(iwp) ::  i          !< running index
     3049    INTEGER(iwp) ::  i_bin      !< result row number
     3050
    30863051    REAL(wp) ::  delta_cold(3)
     3052    REAL(wp) ::  dtmrt          !< delta mean radiant temperature
    30873053    REAL(wp) ::  pmv_cross(2)
    30883054    REAL(wp) ::  reg_a(3)
    30893055    REAL(wp) ::  r_denominator  !< the regression equations denominator
    3090     REAL(wp) ::  dtmrt          !< delta mean radiant temperature
    30913056    REAL(wp) ::  sqrt_ws        !< sqare root of wind speed
    3092     INTEGER(iwp) ::  i          !< running index
    3093     INTEGER(iwp) ::  i_bin      !< result row number
    30943057
    30953058!    REAL(wp) ::  coeff(3,5)  !< unsafe! array is (re-)writable!
     
    31043067!-- Coefficient of the 3 regression lines:
    31053068!      1:const      2:*pmva    3:*ta          4:*sqrt_ws     5:*dtmrt
    3106     REAL(wp), DIMENSION(1:3,1:5), PARAMETER ::  coeff = RESHAPE( (/            &
    3107         0.161_wp,   0.130_wp, -1.125E-03_wp,  1.106E-03_wp, -4.570E-04_wp,     &
    3108         0.795_wp,   0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp,     &
    3109         0.05761_wp, 0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp      &
    3110        /), SHAPE(coeff), order=(/ 2, 1 /) )
     3069    REAL(wp), DIMENSION(1:3,1:5), PARAMETER ::  coeff = RESHAPE( (/                                &
     3070        0.161_wp,   0.130_wp, -1.125E-03_wp,  1.106E-03_wp, -4.570E-04_wp,                         &
     3071        0.795_wp,   0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp,                         &
     3072        0.05761_wp, 0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp                          &
     3073       /), SHAPE( coeff ), order=(/ 2, 1 /)                    )
    31113074!
    31123075!-- Initialise
    31133076    nerr           = 0_iwp
    3114     dpmv_cold_res  = 0._wp
     3077    dpmv_cold_res  = 0.0_wp
    31153078    dtmrt          = tmrt - ta
    31163079    sqrt_ws        = ws
     
    31213084    ENDIF
    31223085
    3123     delta_cold = 0._wp
     3086    delta_cold = 0.0_wp
    31243087    pmv_cross = pmva
    31253088
     
    31273090!-- Determine regression constant for given meteorological conditions
    31283091    DO  i = 1, 3
    3129        reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) *                  &
    3130                   sqrt_ws + coeff(i,5)*dtmrt
     3092       reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) * sqrt_ws + coeff(i,5)*dtmrt
    31313093       delta_cold(i) = reg_a(i) + coeff(i,2) * pmva
    31323094    ENDDO
     
    31523114    ENDDO
    31533115!
    3154 !-- Adjust to operative temperature scaled according
    3155 !-- to classical PMV (Fanger)
     3116!-- Adjust to operative temperature scaled according to classical PMV (Fanger)
    31563117    dpmv_cold_res = delta_cold(i_bin) - dpmv_cold_adj(pmva)
    31573118
    31583119 END SUBROUTINE dpmv_cold
    31593120
    3160 !------------------------------------------------------------------------------!
     3121!--------------------------------------------------------------------------------------------------!
    31613122! Description:
    31623123! ------------
    3163 !> Calculates the summand dpmv_cold_adj adjusting to the operative temperature
    3164 !> scaled according to classical PMV (Fanger) for cold conditions.
    3165 !> Valid for reference environment: v (1m) = 0.10 m/s, dTMRT = 0 K, r.h. = 50 %
    3166 !------------------------------------------------------------------------------!
     3124!> Calculates the summand dpmv_cold_adj adjusting to the operative temperature scaled according to
     3125!> classical PMV (Fanger) for cold conditions. Valid for reference environment: v (1m) = 0.10 m/s,
     3126!> dTMRT = 0 K, r.h. = 50 %
     3127!--------------------------------------------------------------------------------------------------!
    31673128 REAL(wp) FUNCTION dpmv_cold_adj( pmva )
    31683129
    31693130    IMPLICIT NONE
    31703131
    3171     REAL(wp), INTENT ( IN ) ::  pmva        !< (adjusted) Predicted Mean Vote
    3172 
    3173     REAL(wp)      ::  pmv     !< pmv-part of the regression
    31743132    INTEGER(iwp)  ::  i       !< running index
    31753133    INTEGER(iwp)  ::  thr     !< thermal range
     3134
     3135    REAL(wp), INTENT ( IN ) ::  pmva        !< (adjusted) Predicted Mean Vote
     3136
     3137    REAL(wp)  ::  pmv     !< pmv-part of the regression
     3138
    31763139!
    31773140!-- Provide regression coefficients for three thermal ranges:
    3178 !--    slightly cold  cold           very cold
    3179     REAL(wp), DIMENSION(1:3,0:3), PARAMETER ::  coef = RESHAPE( (/             &
    3180        0.0941540_wp, -0.1506620_wp, -0.0871439_wp,                            &
    3181        0.0783162_wp, -1.0612651_wp,  0.1695040_wp,                            &
    3182        0.1350144_wp, -1.0049144_wp, -0.0167627_wp,                            &
    3183        0.1104037_wp, -0.2005277_wp, -0.0003230_wp                              &
    3184        /), SHAPE(coef), order=(/ 1, 2 /) )
     3141!--                                                    slightly cold  cold           very cold
     3142    REAL(wp), DIMENSION(1:3,0:3), PARAMETER ::  coef = RESHAPE( (/                                 &
     3143                                                       0.0941540_wp, -0.1506620_wp, -0.0871439_wp, &
     3144                                                       0.0783162_wp, -1.0612651_wp,  0.1695040_wp, &
     3145                                                       0.1350144_wp, -1.0049144_wp, -0.0167627_wp, &
     3146                                                       0.1104037_wp, -0.2005277_wp, -0.0003230_wp  &
     3147                                                                 /), SHAPE(coef), order=(/ 1, 2 /) )
    31853148!
    31863149!-- Select thermal range
     
    31953158!-- Initialize
    31963159    dpmv_cold_adj = coef(thr,0)
    3197     pmv           = 1._wp
     3160    pmv           = 1.0_wp
    31983161!
    31993162!-- Calculate pmv adjustment (dpmv_cold_adj)
     
    32063169 END FUNCTION dpmv_cold_adj
    32073170
    3208 !------------------------------------------------------------------------------!
     3171!--------------------------------------------------------------------------------------------------!
    32093172! Description:
    32103173! ------------
    3211 !> Based on perceived temperature (perct) as input, ireq_neutral determines
    3212 !> the required clothing insulation (clo) for thermally neutral conditions
    3213 !> (neither body cooling nor body heating). It is related to the Klima-
    3214 !> Michel activity level (134.682 W/m2). IREQ_neutral is only defined
    3215 !> for perct < 10 (degC)
    3216 !------------------------------------------------------------------------------!
     3174!> Based on perceived temperature (perct) as input, ireq_neutral determines the required clothing
     3175!> insulation (clo) for thermally neutral conditions (neither body cooling nor body heating). It is
     3176!> related to the Klima-Michel activity level (134.682 W/m2). IREQ_neutral is only defined for perct
     3177!> < 10 (degC)
     3178!--------------------------------------------------------------------------------------------------!
    32173179 REAL(wp) FUNCTION ireq_neutral( perct_ij, ireq_minimal, nerr )
    32183180
     
    32203182!
    32213183!-- Type declaration of arguments
     3184    INTEGER(iwp), INTENT ( OUT ) ::  nerr
     3185
    32223186    REAL(wp),     INTENT ( IN )  ::  perct_ij
    32233187    REAL(wp),     INTENT ( OUT ) ::  ireq_minimal
    3224     INTEGER(iwp), INTENT ( OUT ) ::  nerr
    32253188!
    32263189!-- Type declaration for internal varables
     
    32543217
    32553218
    3256 !------------------------------------------------------------------------------!
     3219!--------------------------------------------------------------------------------------------------!
    32573220! Description:
    32583221! ------------
    3259 !> The SUBROUTINE surface area calculates the surface area of the individual
    3260 !> according to its height (m), weight (kg), and age (y)
    3261 !------------------------------------------------------------------------------!
     3222!> The SUBROUTINE surface area calculates the surface area of the individual according to its height
     3223!> (m), weight (kg), and age (y)
     3224!--------------------------------------------------------------------------------------------------!
    32623225 SUBROUTINE surface_area( height_cm, weight, age, surf )
    32633226
    32643227    IMPLICIT NONE
    32653228
     3229    INTEGER(iwp), INTENT(in)  ::  age
     3230
     3231    REAL(wp)    , INTENT(in)  ::  height_cm
    32663232    REAL(wp)    , INTENT(in)  ::  weight
    3267     REAL(wp)    , INTENT(in)  ::  height_cm
    3268     INTEGER(iwp), INTENT(in)  ::  age
     3233
    32693234    REAL(wp)    , INTENT(out) ::  surf
     3235
    32703236    REAL(wp)                  ::  height
    32713237
    3272     height = height_cm * 100._wp
     3238    height = height_cm * 100.0_wp
    32733239!
    32743240!-- According to Gehan-George, for children
     
    32823248    ENDIF
    32833249!
    3284 !-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if
    3285 !-- height and weight be known. In: Arch. Int. Med.. 17, 1916, pp. 863:871.
     3250!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if height and weight be
     3251!> known. In: Arch. Int. Med.. 17, 1916, pp. 863:871.
    32863252    surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
    32873253    RETURN
     
    32893255 END SUBROUTINE surface_area
    32903256
    3291 !------------------------------------------------------------------------------!
     3257!--------------------------------------------------------------------------------------------------!
    32923258! Description:
    32933259! ------------
     
    33033269!>  - work load (W)
    33043270!> for a sample human.
    3305 !------------------------------------------------------------------------------!
     3271!--------------------------------------------------------------------------------------------------!
    33063272 SUBROUTINE persdat( age, weight, height, sex, work, a_surf, actlev )
    33073273
    33083274    IMPLICIT NONE
    33093275
     3276    INTEGER(iwp), INTENT(in) ::  sex
     3277
     3278
    33103279    REAL(wp), INTENT(in) ::  age
     3280    REAL(wp), INTENT(in) ::  height
    33113281    REAL(wp), INTENT(in) ::  weight
    3312     REAL(wp), INTENT(in) ::  height
    33133282    REAL(wp), INTENT(in) ::  work
    3314     INTEGER(iwp), INTENT(in) ::  sex
     3283
    33153284    REAL(wp), INTENT(out) ::  actlev
     3285
    33163286    REAL(wp) ::  a_surf
     3287    REAL(wp) ::  basic_heat_prod
    33173288    REAL(wp) ::  energy_prod
     3289    REAL(wp) ::  factor
    33183290    REAL(wp) ::  s
    3319     REAL(wp) ::  factor
    3320     REAL(wp) ::  basic_heat_prod
     3291
    33213292
    33223293    CALL surface_area( height, weight, INT( age ), a_surf )
    3323     s = height * 100._wp / ( weight**( 1._wp / 3._wp ) )
    3324     factor = 1._wp + .004_wp  * ( 30._wp - age )
    3325     basic_heat_prod = 0.
     3294    s = height * 100.0_wp / ( weight**( 1.0_wp / 3.0_wp ) )
     3295    factor = 1.0_wp + .004_wp  * ( 30.0_wp - age )
     3296    basic_heat_prod = 0.0_wp
    33263297    IF ( sex == 1_iwp )  THEN
    3327        basic_heat_prod = 3.45_wp * weight**( 3._wp / 4._wp ) * ( factor +      &
    3328                      .01_wp  * ( s - 43.4_wp ) )
     3298       basic_heat_prod = 3.45_wp * weight**( 3.0_wp / 4.0_wp ) * ( factor + 0.01_wp                &
     3299                         * ( s - 43.4_wp ) )
    33293300    ELSE IF ( sex == 2_iwp )  THEN
    3330        basic_heat_prod = 3.19_wp * weight**( 3._wp / 4._wp ) * ( factor +      &
    3331                     .018_wp * ( s - 42.1_wp ) )
     3301       basic_heat_prod = 3.19_wp * weight**( 3.0_wp / 4.0_wp ) * ( factor + 0.018_wp               &
     3302                        * ( s - 42.1_wp ) )
    33323303    ENDIF
    33333304
     
    33383309
    33393310
    3340 !------------------------------------------------------------------------------!
     3311!--------------------------------------------------------------------------------------------------!
    33413312! Description:
    33423313! ------------
    33433314!> SUBROUTINE ipt_init
    33443315!> initializes the instationary perceived temperature
    3345 !------------------------------------------------------------------------------!
    3346 
    3347  SUBROUTINE ipt_init( age, weight, height, sex, work, actlev, clo,             &
    3348      ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                          &
    3349      ipt )
     3316!--------------------------------------------------------------------------------------------------!
     3317
     3318 SUBROUTINE ipt_init( age, weight, height, sex, work, actlev, clo, ta, vp, ws, tmrt, pair, dt,     &
     3319                      storage, t_clothing, ipt )
    33503320
    33513321    IMPLICIT NONE
    33523322!
    33533323!-- Input parameters
     3324
     3325    INTEGER(iwp), INTENT(in)  :: sex    !< Persons sex (1 = male, 2 = female)
     3326
    33543327    REAL(wp), INTENT(in) ::  age        !< Persons age          (years)
     3328    REAL(wp), INTENT(in) ::  dt         !< Timestep             (s)
     3329    REAL(wp), INTENT(in) ::  height     !< Persons height       (m)7
     3330    REAL(wp), INTENT(in) ::  pair       !< Air pressure         (hPa)
     3331    REAL(wp), INTENT(in) ::  ta         !< Air Temperature      (degree_C)
     3332    REAL(wp), INTENT(in) ::  tmrt       !< Mean radiant temperature   (degree_C)
     3333    REAL(wp), INTENT(in) ::  vp         !< Vapor pressure       (hPa)
    33553334    REAL(wp), INTENT(in) ::  weight     !< Persons weight       (kg)
    3356     REAL(wp), INTENT(in) ::  height     !< Persons height       (m)
    33573335    REAL(wp), INTENT(in) ::  work       !< Current workload     (W)
    3358     REAL(wp), INTENT(in) ::  ta         !< Air Temperature      (degree_C)
    3359     REAL(wp), INTENT(in) ::  vp         !< Vapor pressure       (hPa)
    33603336    REAL(wp), INTENT(in) ::  ws         !< Wind speed in approx. 1.1m (m/s)
    3361     REAL(wp), INTENT(in) ::  tmrt       !< Mean radiant temperature   (degree_C)
    3362     REAL(wp), INTENT(in) ::  pair       !< Air pressure         (hPa)
    3363     REAL(wp), INTENT(in) ::  dt         !< Timestep             (s)
    3364     INTEGER(iwp), INTENT(in)  :: sex    !< Persons sex (1 = male, 2 = female)
    33653337!
    33663338!-- Output parameters
    33673339    REAL(wp), INTENT(out) ::  actlev
    33683340    REAL(wp), INTENT(out) ::  clo
     3341    REAL(wp), INTENT(out) ::  ipt
    33693342    REAL(wp), INTENT(out) ::  storage
    33703343    REAL(wp), INTENT(out) ::  t_clothing
    3371     REAL(wp), INTENT(out) ::  ipt
    33723344!
    33733345!-- Internal variables
    3374     REAL(wp), PARAMETER :: eps = 0.0005_wp
    3375     REAL(wp), PARAMETER :: eta = 0._wp
    3376     REAL(wp) ::  sclo
    3377     REAL(wp) ::  wclo
    3378     REAL(wp) ::  d_pmv
    3379     REAL(wp) ::  svp_ta
    3380     REAL(wp) ::  sult_lim
    3381     REAL(wp) ::  dgtcm
    3382     REAL(wp) ::  dgtcstd
    3383     REAL(wp) ::  clon
    3384     REAL(wp) ::  ireq_minimal
    3385 !     REAL(wp) ::  clo_fanger
    3386     REAL(wp) ::  pmv_w
    3387     REAL(wp) ::  pmv_s
    3388     REAL(wp) ::  pmva
    3389     REAL(wp) ::  ptc
    3390     REAL(wp) ::  d_std
    3391     REAL(wp) ::  pmvs
    3392     REAL(wp) ::  a_surf
    3393 !     REAL(wp) ::  acti
    33943346    INTEGER(iwp) ::  ncount
    33953347    INTEGER(iwp) ::  nerr_cold
     
    33983350    LOGICAL ::  sultrieness
    33993351
    3400     storage = 0._wp
     3352    REAL(wp), PARAMETER :: eps = 0.0005_wp
     3353    REAL(wp), PARAMETER :: eta = 0.0_wp
     3354
     3355!    REAL(wp) ::  acti
     3356    REAL(wp) ::  a_surf
     3357!    REAL(wp) ::  clo_fanger
     3358    REAL(wp) ::  clon
     3359    REAL(wp) ::  d_pmv
     3360    REAL(wp) ::  d_std
     3361    REAL(wp) ::  dgtcm
     3362    REAL(wp) ::  dgtcstd
     3363    REAL(wp) ::  ireq_minimal
     3364    REAL(wp) ::  pmv_s
     3365    REAL(wp) ::  pmv_w
     3366    REAL(wp) ::  pmva
     3367    REAL(wp) ::  pmvs
     3368    REAL(wp) ::  ptc
     3369    REAL(wp) ::  sclo
     3370    REAL(wp) ::  sult_lim
     3371    REAL(wp) ::  svp_ta
     3372    REAL(wp) ::  wclo
     3373
     3374
     3375    storage = 0.0_wp
    34013376    CALL persdat( age, weight, height, sex, work, a_surf, actlev )
    34023377!
     
    34153390!
    34163391!-- Decision: firstly calculate for winter or summer clothing
    3417     IF ( ta <= 10._wp )  THEN
     3392    IF ( ta <= 10.0_wp )  THEN
    34183393!
    34193394!--    First guess: winter clothing insulation: cold stress
    34203395       clo = wclo
    34213396       t_clothing = bio_fill_value  ! force initial run
    3422        CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
    3423           t_clothing, storage, dt, pmva )
     3397       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt,    &
     3398                            pmva )
    34243399       pmv_w = pmva
    34253400
    3426        IF ( pmva > 0._wp )  THEN
    3427 !
    3428 !--       Case summer clothing insulation: heat load ?           
     3401       IF ( pmva > 0.0_wp )  THEN
     3402!
     3403!--       Case summer clothing insulation: heat load ?
    34293404          clo = sclo
    34303405          t_clothing = bio_fill_value  ! force initial run
    3431           CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,     &
    3432              t_clothing, storage, dt, pmva )
     3406          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, &
     3407                               pmva )
    34333408          pmv_s = pmva
    3434           IF ( pmva <= 0._wp )  THEN
    3435 !
    3436 !--          Case: comfort achievable by varying clothing insulation
    3437 !--          between winter and summer set values
    3438              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo,     &
    3439                             pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
     3409          IF ( pmva <= 0.0_wp )  THEN
     3410!
     3411!--          Case: comfort achievable by varying clothing insulation between winter and summer set
     3412!--                values
     3413             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo, pmv_s, wclo, pmv_w, eps,&
     3414                              pmva, ncount, clo )
    34403415             IF ( ncount < 0_iwp )  THEN
    34413416                nerr = -1_iwp
     
    34453420             clo = 0.5_wp
    34463421             t_clothing = bio_fill_value
    3447              CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
    3448                 t_clothing, storage, dt, pmva )
     3422             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage,  &
     3423                                  dt, pmva )
    34493424          ENDIF
    3450        ELSE IF ( pmva < -0.11_wp )  THEN
     3425       ELSE IF ( pmva < - 0.11_wp )  THEN
    34513426          clo = 1.75_wp
    34523427          t_clothing = bio_fill_value
    3453           CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,     &
    3454              t_clothing, storage, dt, pmva )
     3428          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, &
     3429                               pmva )
    34553430       ENDIF
    34563431
     
    34603435       clo = sclo
    34613436       t_clothing = bio_fill_value
    3462        CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
    3463           t_clothing, storage, dt, pmva )
     3437       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt,    &
     3438                            pmva )
    34643439       pmv_s = pmva
    34653440
    3466        IF ( pmva < 0._wp )  THEN
     3441       IF ( pmva < 0.0_wp )  THEN
    34673442!
    34683443!--       Case winter clothing insulation: cold stress ?
    34693444          clo = wclo
    34703445          t_clothing = bio_fill_value
    3471           CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,     &
    3472              t_clothing, storage, dt, pmva )
     3446          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, &
     3447                               pmva )
    34733448          pmv_w = pmva
    34743449
    3475           IF ( pmva >= 0._wp )  THEN
    3476 !
    3477 !--          Case: comfort achievable by varying clothing insulation
    3478 !--          between winter and summer set values
    3479              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,     &
    3480                                pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
     3450          IF ( pmva >= 0.0_wp )  THEN
     3451!
     3452!--          Case: comfort achievable by varying clothing insulation between winter and summer set
     3453!--                values
     3454             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, pmv_s, wclo, pmv_w, eps, &
     3455                               pmva, ncount, clo )
    34813456             IF ( ncount < 0_wp )  THEN
    34823457                nerr = -1_iwp
    34833458                RETURN
    34843459             ENDIF
    3485           ELSE IF ( pmva < -0.11_wp )  THEN
     3460          ELSE IF ( pmva < - 0.11_wp )  THEN
    34863461             clo = 1.75_wp
    34873462             t_clothing = bio_fill_value
    3488              CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
    3489                 t_clothing, storage, dt, pmva )
     3463             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage,  &
     3464                                  dt, pmva )
    34903465          ENDIF
    34913466       ELSE IF ( pmva > 0.06_wp )  THEN
    34923467          clo = 0.5_wp
    34933468          t_clothing = bio_fill_value
    3494           CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,     &
    3495              t_clothing, storage, dt, pmva )
     3469          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, &
     3470                               pmva )
    34963471       ENDIF
    34973472
     
    35023477    CALL perct_regression( pmva, clo, ipt )
    35033478    ptc = ipt
    3504     IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
     3479    IF ( clo >= 1.75_wp  .AND.  pmva <= - 0.11_wp )  THEN
    35053480!
    35063481!--    Adjust for cold conditions according to Gagge 1986
     
    35083483       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
    35093484       pmvs = pmva - d_pmv
    3510        IF ( pmvs > -0.11_wp )  THEN
    3511           d_pmv  = 0._wp
    3512           pmvs   = -0.11_wp
     3485       IF ( pmvs > - 0.11_wp )  THEN
     3486          d_pmv  = 0.0_wp
     3487          pmvs   = - 0.11_wp
    35133488       ENDIF
    35143489       CALL perct_regression( pmvs, clo, ipt )
     
    35183493    IF ( clo > 0.5_wp  .AND.  ipt <= 8.73_wp )  THEN
    35193494!
    3520 !--    Required clothing insulation (ireq) is exclusively defined for
    3521 !--    perceived temperatures (ipt) less 10 (C) for a
    3522 !--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
     3495!--    Required clothing insulation (ireq) is exclusively defined for perceived temperatures (ipt)
     3496!--    less 10 (C) for a reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
    35233497       clon = ireq_neutral ( ipt, ireq_minimal, nerr )
    35243498       clo = clon
     
    35263500    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
    35273501    sultrieness    = .FALSE.
    3528     d_std      = -99._wp
     3502    d_std      = - 99.0_wp
    35293503    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
    35303504!
     
    35343508       pmvs   = pmva + d_pmv
    35353509       CALL perct_regression( pmvs, clo, ipt )
    3536        IF ( sult_lim < 99._wp )  THEN
     3510       IF ( sult_lim < 99.0_wp )  THEN
    35373511          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
    35383512       ENDIF
    35393513    ENDIF
    35403514
    3541  
     3515
    35423516 END SUBROUTINE ipt_init
    3543  
    3544 !------------------------------------------------------------------------------!
     3517
     3518!--------------------------------------------------------------------------------------------------!
    35453519! Description:
    35463520! ------------
    35473521!> SUBROUTINE ipt_cycle
    3548 !> Calculates one timestep for the instationary version of perceived
    3549 !> temperature (iPT, degree_C) for
    3550 !>  - standard measured/predicted meteorological values and TMRT
    3551 !>    as input;
     3522!> Calculates one timestep for the instationary version of perceived temperature (iPT, degree_C) for
     3523!>  - standard measured/predicted meteorological values and TMRT as input;
    35523524!>  - regressions for determination of PT;
    3553 !>  - adjustment to Gagge's PMV* (2-node-model, 1986) as base of PT
    3554 !>    under warm/humid conditions (Icl= 0.50 clo) and under cold
    3555 !>    conditions (Icl= 1.75 clo)
    3556 !>
    3557 !------------------------------------------------------------------------------!
    3558  SUBROUTINE ipt_cycle( ta, vp, ws, tmrt, pair, dt, storage, t_clothing, clo,   &
    3559      actlev, work, ipt )
     3525!>  - adjustment to Gagge's PMV* (2-node-model, 1986) as base of PT under warm/humid conditions
     3526!>    (Icl= 0.50 clo) and under cold conditions (Icl= 1.75 clo)
     3527!--------------------------------------------------------------------------------------------------!
     3528 SUBROUTINE ipt_cycle( ta, vp, ws, tmrt, pair, dt, storage, t_clothing, clo, actlev, work, ipt )
    35603529
    35613530    IMPLICIT NONE
    35623531!
    35633532!-- Type of input of the argument list
     3533    REAL(wp), INTENT ( IN )  ::  actlev  !< Internal heat production    (W)
     3534    REAL(wp), INTENT ( IN )  ::  clo     !< Clothing index              (no dim)
     3535    REAL(wp), INTENT ( IN )  ::  dt      !< Timestep                    (s)
     3536    REAL(wp), INTENT ( IN )  ::  pair    !< Air pressure                (hPa)
    35643537    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature             (degree_C)
     3538    REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature    (degree_C)
    35653539    REAL(wp), INTENT ( IN )  ::  vp      !< Vapor pressure              (hPa)
    3566     REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature    (degree_C)
     3540    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical work load        (W)
    35673541    REAL(wp), INTENT ( IN )  ::  ws      !< Wind speed                  (m/s)
    3568     REAL(wp), INTENT ( IN )  ::  pair    !< Air pressure                (hPa)
    3569     REAL(wp), INTENT ( IN )  ::  dt      !< Timestep                    (s)
    3570     REAL(wp), INTENT ( IN )  ::  clo     !< Clothing index              (no dim)
    3571     REAL(wp), INTENT ( IN )  ::  actlev  !< Internal heat production    (W)
    3572     REAL(wp), INTENT ( IN )  ::  work    !< Mechanical work load        (W)
    35733542!
    35743543!-- In and output parameters
     
    35803549!
    35813550!-- Type of internal variables
     3551    INTEGER(iwp) ::  nerr
     3552    INTEGER(iwp) ::  nerr_cold
     3553
     3554    LOGICAL ::  sultrieness
     3555
    35823556    REAL(wp) ::  d_pmv
    3583     REAL(wp) ::  svp_ta
    3584     REAL(wp) ::  sult_lim
     3557    REAL(wp) ::  d_std
    35853558    REAL(wp) ::  dgtcm
    35863559    REAL(wp) ::  dgtcstd
    35873560    REAL(wp) ::  pmva
     3561    REAL(wp) ::  pmvs
    35883562    REAL(wp) ::  ptc
    3589     REAL(wp) ::  d_std
    3590     REAL(wp) ::  pmvs
    3591     INTEGER(iwp) ::  nerr_cold
    3592     INTEGER(iwp) ::  nerr
    3593 
    3594     LOGICAL ::  sultrieness
     3563    REAL(wp) ::  sult_lim
     3564    REAL(wp) ::  svp_ta
    35953565!
    35963566!-- Initialise
     
    36013571!
    36023572!-- Determine pmv_adjusted for current conditions
    3603     CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,            &
    3604        t_clothing, storage, dt, pmva )
     3573    CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, pmva )
    36053574!
    36063575!-- Determine perceived temperature by regression equation + adjustments
     
    36143583       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
    36153584       pmvs = pmva - d_pmv
    3616        IF ( pmvs > -0.11_wp )  THEN
    3617           d_pmv  = 0._wp
    3618           pmvs   = -0.11_wp
     3585       IF ( pmvs > - 0.11_wp )  THEN
     3586          d_pmv  = 0.0_wp
     3587          pmvs   = - 0.11_wp
    36193588       ENDIF
    36203589       CALL perct_regression( pmvs, clo, ipt )
     
    36253594    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
    36263595    sultrieness = .FALSE.
    3627     d_std       = -99._wp
     3596    d_std       = - 99.0_wp
    36283597    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
    36293598!
     
    36333602       pmvs  = pmva + d_pmv
    36343603       CALL perct_regression( pmvs, clo, ipt )
    3635        IF ( sult_lim < 99._wp )  THEN
     3604       IF ( sult_lim < 99.0_wp )  THEN
    36363605          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
    36373606       ENDIF
     
    36403609 END SUBROUTINE ipt_cycle
    36413610
    3642 !------------------------------------------------------------------------------!
     3611!--------------------------------------------------------------------------------------------------!
    36433612! Description:
    36443613! ------------
    3645 !> SUBROUTINE fanger_s calculates the
    3646 !> actual Predicted Mean Vote (dimensionless) according
    3647 !> to Fanger corresponding to meteorological (ta,tmrt,pa,ws,pair)
    3648 !> and individual variables (clo, actlev, eta) considering a storage
    3649 !> and clothing temperature for a given timestep.
    3650 !------------------------------------------------------------------------------!
    3651  SUBROUTINE fanger_s_acti( ta, tmrt, pa, in_ws, pair, in_clo, actlev,          &
    3652     activity, t_cloth, s, dt, pmva )
     3614!> SUBROUTINE fanger_s calculates the actual Predicted Mean Vote (dimensionless) according to Fanger
     3615!> corresponding to meteorological (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
     3616!> considering a storage and clothing temperature for a given timestep.
     3617!--------------------------------------------------------------------------------------------------!
     3618 SUBROUTINE fanger_s_acti( ta, tmrt, pa, in_ws, pair, in_clo, actlev, activity, t_cloth, s, dt,    &
     3619                           pmva )
    36533620
    36543621    IMPLICIT NONE
    36553622!
    36563623!--  Input argument types
     3624    REAL(wp), INTENT ( IN )  ::  activity !< Work load                (W/m²)
     3625    REAL(wp), INTENT ( IN )  ::  actlev   !< Metabolic + work energy  (W/m²)
     3626    REAL(wp), INTENT ( IN )  ::  dt       !< Timestep                 (s)
     3627    REAL(wp), INTENT ( IN )  ::  in_clo   !< Clothing index (clo)     (no dim)
     3628    REAL(wp), INTENT ( IN )  ::  in_ws    !< Wind speed               (m/s)
     3629    REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
     3630    REAL(wp), INTENT ( IN )  ::  pair     !< Air pressure             (hPa)
    36573631    REAL(wp), INTENT ( IN )  ::  ta       !< Air temperature          (degree_C)
    36583632    REAL(wp), INTENT ( IN )  ::  tmrt     !< Mean radiant temperature (degree_C)
    3659     REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
    3660     REAL(wp), INTENT ( IN )  ::  pair     !< Air pressure             (hPa)
    3661     REAL(wp), INTENT ( IN )  ::  in_ws    !< Wind speed               (m/s)
    3662     REAL(wp), INTENT ( IN )  ::  actlev   !< Metabolic + work energy  (W/m²)
    3663     REAL(wp), INTENT ( IN )  ::  dt       !< Timestep                 (s)
    3664     REAL(wp), INTENT ( IN )  ::  activity !< Work load                (W/m²)
    3665     REAL(wp), INTENT ( IN )  ::  in_clo   !< Clothing index (clo)     (no dim)
    36663633!
    36673634!-- Output argument types
     
    36723639!
    36733640!-- Internal variables
    3674     REAL(wp), PARAMETER  ::  time_equil = 7200._wp
    3675 
    3676     REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
    3677     REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
    3678     REAL(wp) ::  t_skin_aver  !< average skin temperature               (degree_C)
    3679     REAL(wp) ::  bc           !< preliminary result storage
    3680     REAL(wp) ::  cc           !< preliminary result storage
    3681     REAL(wp) ::  dc           !< preliminary result storage
    3682     REAL(wp) ::  ec           !< preliminary result storage
    3683     REAL(wp) ::  gc           !< preliminary result storage
    3684     REAL(wp) ::  t_clothing   !< clothing temperature                   (degree_C)
    3685 !     REAL(wp) ::  hr           !< radiational heat resistence
    3686     REAL(wp) ::  clo          !< clothing insulation index              (clo)
    3687     REAL(wp) ::  ws           !< wind speed                             (m/s)
    3688     REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat
    3689                               !< ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
    3690     REAL(wp) ::  z2           !< Water vapour diffution through the skin
    3691     REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
    3692     REAL(wp) ::  z4           !< Loss of latent heat through respiration
    3693     REAL(wp) ::  z5           !< Loss of radiational heat
    3694     REAL(wp) ::  z6           !< Heat loss through forced convection
    3695     REAL(wp) ::  en           !< Energy ballance                        (W)
    3696     REAL(wp) ::  d_s          !< Storage delta                          (W)
    3697     REAL(wp) ::  adjustrate   !< Max storage adjustment rate
    3698     REAL(wp) ::  adjustrate_cloth  !< max clothing temp. adjustment rate
    3699 
    37003641    INTEGER(iwp) :: i         !< running index
    37013642    INTEGER(iwp) ::  niter    !< Running index
    37023643
     3644    REAL(wp), PARAMETER  ::  time_equil = 7200.0_wp
     3645
     3646    REAL(wp) ::  adjustrate        !< Max storage adjustment rate
     3647    REAL(wp) ::  adjustrate_cloth  !< max clothing temp. adjustment rate
     3648    REAL(wp) ::  bc                !< preliminary result storage
     3649    REAL(wp) ::  cc                !< preliminary result storage
     3650    REAL(wp) ::  clo               !< clothing insulation index              (clo)
     3651    REAL(wp) ::  d_s               !< Storage delta                          (W)
     3652    REAL(wp) ::  dc                !< preliminary result storage
     3653    REAL(wp) ::  en                !< Energy ballance                        (W)
     3654    REAL(wp) ::  ec                !< preliminary result storage
     3655    REAL(wp) ::  f_cl              !< Increase in surface due to clothing    (factor)
     3656    REAL(wp) ::  gc                !< preliminary result storage
     3657    REAL(wp) ::  heat_convection   !< energy loss by autocnvection       (W)
     3658!    REAL(wp) ::  hr                !< radiational heat resistence
     3659    REAL(wp) ::  t_clothing        !< clothing temperature                   (degree_C)
     3660    REAL(wp) ::  t_skin_aver       !< average skin temperature               (degree_C)
     3661    REAL(wp) ::  ws                !< wind speed                             (m/s)
     3662    REAL(wp) ::  z1                !< Empiric factor for the adaption of the heat
     3663                                   !< ballance equation to the psycho-physical scale
     3664                                   !< (Equ. 40 in FANGER)
     3665    REAL(wp) ::  z2                !< Water vapour diffution through the skin
     3666    REAL(wp) ::  z3                !< Sweat evaporation from the skin surface
     3667    REAL(wp) ::  z4                !< Loss of latent heat through respiration
     3668    REAL(wp) ::  z5                !< Loss of radiational heat
     3669    REAL(wp) ::  z6                !< Heat loss through forced convection
     3670
     3671
     3672
     3673
    37033674!
    37043675!-- Clo must be > 0. to avoid div. by 0!
    37053676    clo = in_clo
    3706     IF ( clo < 001._wp )  clo = .001_wp
     3677    IF ( clo < 001.0_wp )  clo = 0.001_wp
    37073678!
    37083679!-- Increase in surface due to clothing
    3709     f_cl = 1._wp + .15_wp * clo
     3680    f_cl = 1.0_wp + 0.15_wp * clo
    37103681!
    37113682!-- Case of free convection (ws < 0.1 m/s ) not considered
    37123683    ws = in_ws
    3713     IF ( ws < .1_wp )  THEN
    3714        ws = .1_wp
     3684    IF ( ws < 0.1_wp )  THEN
     3685       ws = 0.1_wp
    37153686    ENDIF
    37163687!
     
    37193690!
    37203691!-- Average skin temperature
    3721     t_skin_aver = 35.7_wp - .0275_wp * activity
     3692    t_skin_aver = 35.7_wp - 0.0275_wp * activity
    37223693!
    37233694!-- Calculation of constants for evaluation below
    3724     bc = .155_wp * clo * 3.96_wp * 10._wp**( -8._wp ) * f_cl
     3695    bc = 0.155_wp * clo * 3.96_wp * 10.0_wp**( -8.0_wp ) * f_cl
    37253696    cc = f_cl * heat_convection
    3726     ec = .155_wp * clo
    3727     dc = ( 1._wp + ec * cc ) / bc
    3728     gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4._wp + ec * cc * ta ) / bc
    3729 !
    3730 !-- Calculation of clothing surface temperature (t_clothing) based on
    3731 !-- newton-approximation with air temperature as initial guess
    3732     niter = INT( dt * 10._wp, KIND=iwp )
     3697    ec = 0.155_wp * clo
     3698    dc = ( 1.0_wp + ec * cc ) / bc
     3699    gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4.0_wp + ec * cc * ta ) / bc
     3700!
     3701!-- Calculation of clothing surface temperature (t_clothing) based on Newton-approximation with air
     3702!-- temperature as initial guess
     3703    niter = INT( dt * 10.0_wp, KIND=iwp )
    37333704    IF ( niter < 1 )  niter = 1_iwp
    3734     adjustrate = 1._wp - EXP( -1._wp * ( 10._wp / time_equil ) * dt )
    3735     IF ( adjustrate >= 1._wp )  adjustrate = 1._wp
    3736     adjustrate_cloth = adjustrate * 30._wp
     3705    adjustrate = 1.0_wp - EXP( -1.0_wp * ( 10.0_wp / time_equil ) * dt )
     3706    IF ( adjustrate >= 1.0_wp )  adjustrate = 1.0_wp
     3707    adjustrate_cloth = adjustrate * 30.0_wp
    37373708    t_clothing = t_cloth
    37383709!
    3739 !-- Set initial values for niter, adjustrates and t_clothing if this is the
    3740 !-- first call
    3741     IF ( t_cloth <= -998._wp )  THEN  ! If initial run
     3710!-- Set initial values for niter, adjustrates and t_clothing if this is the first call
     3711    IF ( t_cloth <= -998.0_wp )  THEN  ! If initial run
    37423712       niter = 3_iwp
    3743        adjustrate = 1._wp
    3744        adjustrate_cloth = 1._wp
     3713       adjustrate = 1.0_wp
     3714       adjustrate_cloth = 1.0_wp
    37453715       t_clothing = ta
    37463716    ENDIF
     
    37483718!-- Update clothing temperature
    37493719    DO  i = 1, niter
    3750        t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
    3751           273.2_wp )**4._wp + t_clothing *                                     &
    3752           dc - gc ) / ( 4._wp * ( t_clothing + 273.2_wp )**3._wp + dc )
     3720       t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing + 273.2_wp )**4.0_wp  +         &
     3721                    t_clothing * dc - gc ) / ( 4.0_wp * ( t_clothing + 273.2_wp )**3.0_wp + dc )
    37533722    ENDDO
    37543723!
    3755 !-- Empiric factor for the adaption of the heat ballance equation
    3756 !-- to the psycho-physical scale (Equ. 40 in FANGER)
    3757     z1 = ( .303_wp * EXP( -.036_wp * actlev ) + .0275_wp )
     3724!-- Empiric factor for the adaption of the heat ballance equation to the psycho-physical scale
     3725!-- (Equ. 40 in FANGER)
     3726    z1 = ( 0.303_wp * EXP( - 0.036_wp * actlev ) + 0.0275_wp )
    37583727!
    37593728!-- Water vapour diffution through the skin
    3760     z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
     3729    z2 = 0.31_wp * ( 57.3_wp - 0.07_wp * activity-pa )
    37613730!
    37623731!-- Sweat evaporation from the skin surface
    3763     z3 = .42_wp * ( activity - 58._wp )
     3732    z3 = 0.42_wp * ( activity - 58.0_wp )
    37643733!
    37653734!-- Loss of latent heat through respiration
    3766     z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
    3767       ( 34._wp - ta )
     3735    z4 = 0.0017_wp * actlev * ( 58.7_wp - pa ) + 0.0014_wp * actlev * ( 34.0_wp - ta )
    37683736!
    37693737!-- Loss of radiational heat
    3770     z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + 273.2_wp )**4 - ( tmrt +         &
    3771        273.2_wp )**4 )
     3738    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + 273.2_wp )**4 - ( tmrt + 273.2_wp )**4 )
    37723739!
    37733740!-- Heat loss through forced convection
     
    37783745!
    37793746!-- Manage storage
    3780     d_s = adjustrate * en + ( 1._wp - adjustrate ) * s
     3747    d_s = adjustrate * en + ( 1.0_wp - adjustrate ) * s
    37813748!
    37823749!-- Predicted Mean Vote
     
    37913758
    37923759
    3793 !------------------------------------------------------------------------------!
     3760!--------------------------------------------------------------------------------------------------!
    37943761!
    37953762! Description:
     
    37983765!> stationary (calculated based on MEMI),
    37993766!> Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
    3800 !------------------------------------------------------------------------------!
     3767!--------------------------------------------------------------------------------------------------!
    38013768
    38023769 SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet_ij )
     
    38053772!
    38063773!-- Input arguments:
     3774    REAL(wp), INTENT( IN ) ::  pair  !< Air pressure                (hPa)
    38073775    REAL(wp), INTENT( IN ) ::  ta    !< Air temperature             (degree_C)
    38083776    REAL(wp), INTENT( IN ) ::  tmrt  !< Mean radiant temperature    (degree_C)
    38093777    REAL(wp), INTENT( IN ) ::  v     !< Wind speed                  (m/s)
    38103778    REAL(wp), INTENT( IN ) ::  vpa   !< Vapor pressure              (hPa)
    3811     REAL(wp), INTENT( IN ) ::  pair  !< Air pressure                (hPa)
    38123779!
    38133780!-- Output arguments:
     
    38273794    REAL(wp) ::  rtv
    38283795    REAL(wp) ::  vpts       !< Sat. vapor pressure over skin        (hPa)
     3796    REAL(wp) ::  tcl        !< Clothing temperature                 (degree_C)
    38293797    REAL(wp) ::  tsk        !< Skin temperature                     (degree_C)
    3830     REAL(wp) ::  tcl        !< Clothing temperature                 (degree_C)
    38313798    REAL(wp) ::  wetsk      !< Fraction of wet skin                 (factor)
    38323799!
     
    38363803!-- MEMI configuration
    38373804    REAL(wp) :: age         !< Persons age          (a)
     3805    REAL(wp) :: clo         !< Clothing insulation index (clo)
     3806    REAL(wp) :: eta         !< Work efficiency      (dimensionless)
     3807    REAL(wp) :: fcl         !< Surface area modification by clothing (factor)
     3808    REAL(wp) :: ht          !< Persons height       (m)
    38383809    REAL(wp) :: mbody       !< Persons body mass    (kg)
    3839     REAL(wp) :: ht          !< Persons height       (m)
    38403810    REAL(wp) :: work        !< Work load            (W)
    3841     REAL(wp) :: eta         !< Work efficiency      (dimensionless)
    3842     REAL(wp) :: clo         !< Clothing insulation index (clo)
    3843     REAL(wp) :: fcl         !< Surface area modification by clothing (factor)
    3844 !     INTEGER(iwp) :: pos     !< Posture: 1 = standing, 2 = sitting
    3845 !     INTEGER(iwp) :: sex     !< Sex: 1 = male, 2 = female
     3811!    INTEGER(iwp) :: pos     !< Posture: 1 = standing, 2 = sitting
     3812!    INTEGER(iwp) :: sex     !< Sex: 1 = male, 2 = female
    38463813!
    38473814!-- Configuration, keep standard parameters!
    3848     age   = 35._wp
    3849     mbody = 75._wp
     3815    age   = 35.0_wp
     3816    mbody = 75.0_wp
    38503817    ht    =  1.75_wp
    3851     work  = 80._wp
    3852     eta   =  0._wp
     3818    work  = 80.0_wp
     3819    eta   =  0.0_wp
    38533820    clo   =  0.9_wp
    38543821    fcl   =  1.15_wp
    38553822!
    38563823!-- Call subfunctions
    3857     CALL in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,     &
    3858             vpa, work )
    3859 
    3860     CALL heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht,  &
    3861             int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,     &
    3862             vpts, wetsk )
    3863 
    3864     CALL pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair,       &
    3865             rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
     3824    CALL in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, vpa, work )
     3825
     3826    CALL heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, int_heat, mbody,     &
     3827                    pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
     3828
     3829    CALL pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair, rdcl, rdsk, rtv, ta, tcl, &
     3830                        tsk, pet_ij, vpts, wetsk )
    38663831
    38673832
     
    38693834
    38703835
    3871 !------------------------------------------------------------------------------!
     3836!--------------------------------------------------------------------------------------------------!
    38723837! Description:
    38733838! ------------
    38743839!> Calculate internal energy ballance
    3875 !------------------------------------------------------------------------------!
    3876  SUBROUTINE in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
    3877     vpa, work )
     3840!--------------------------------------------------------------------------------------------------!
     3841 SUBROUTINE in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, vpa, work )
    38783842!
    38793843!-- Input arguments:
     3844    REAL(wp), INTENT( IN )  ::  age       !< Persons age              (a)
     3845    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency     (dimensionless)
     3846    REAL(wp), INTENT( IN )  ::  ht        !< Persons height           (m)
     3847    REAL(wp), INTENT( IN )  ::  mbody     !< Persons body mass        (kg)
    38803848    REAL(wp), INTENT( IN )  ::  pair      !< air pressure             (hPa)
    38813849    REAL(wp), INTENT( IN )  ::  ta        !< air temperature          (degree_C)
    38823850    REAL(wp), INTENT( IN )  ::  vpa       !< vapor pressure           (hPa)
    3883     REAL(wp), INTENT( IN )  ::  age       !< Persons age              (a)
    3884     REAL(wp), INTENT( IN )  ::  mbody     !< Persons body mass        (kg)
    3885     REAL(wp), INTENT( IN )  ::  ht        !< Persons height           (m)
    38863851    REAL(wp), INTENT( IN )  ::  work      !< Work load                (W)
    3887     REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency     (dimensionless)
    38883852!
    38893853!-- Output arguments:
     
    39013865!
    39023866!-- Metabolic heat production
    3903     met = 3.45_wp * mbody**( 3._wp / 4._wp ) * (1._wp + 0.004_wp *             &
    3904           ( 30._wp - age) + 0.010_wp * ( ( ht * 100._wp /                      &
    3905           ( mbody**( 1._wp / 3._wp ) ) ) - 43.4_wp ) )
     3867    met = 3.45_wp * mbody**( 3.0_wp / 4.0_wp ) * (1.0_wp + 0.004_wp *                              &
     3868          ( 30.0_wp - age) + 0.010_wp * ( ( ht * 100.0_wp /                                        &
     3869          ( mbody**( 1.0_wp / 3.0_wp ) ) ) - 43.4_wp ) )
    39063870    met = work + met
    3907     int_heat = met * (1._wp - eta)
     3871    int_heat = met * (1.0_wp - eta)
    39083872!
    39093873!-- Sensible respiration energy
    39103874    tex  = 0.47_wp * ta + 21.0_wp
    3911     rtv  = 1.44_wp * 10._wp**(-6._wp) * met
     3875    rtv  = 1.44_wp * 10.0_wp**(-6.0_wp) * met
    39123876    eres = c_p * (ta - tex) * rtv
    39133877!
    39143878!-- Latent respiration energy
    3915     vpex = 6.11_wp * 10._wp**( 7.45_wp * tex / ( 235._wp + tex ) )
     3879    vpex = 6.11_wp * 10.0_wp**( 7.45_wp * tex / ( 235.0_wp + tex ) )
    39163880    erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv
    39173881!
     
    39223886
    39233887
    3924 !------------------------------------------------------------------------------!
     3888!--------------------------------------------------------------------------------------------------!
    39253889! Description:
    39263890! ------------
    39273891!> Calculate heat gain or loss
    3928 !------------------------------------------------------------------------------!
    3929  SUBROUTINE heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff,   &
    3930         ht, int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,     &
    3931         vpts, wetsk )
     3892!--------------------------------------------------------------------------------------------------!
     3893 SUBROUTINE heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, int_heat, mbody,  &
     3894                       pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
    39323895
    39333896!
    39343897!-- Input arguments:
     3898    REAL(wp), INTENT( IN )  ::  clo    !< clothing insulation      (clo)
     3899    REAL(wp), INTENT( IN )  ::  fcl    !< factor for surface area increase by clothing
    39353900    REAL(wp), INTENT( IN )  ::  ere    !< Energy ballance          (W)
    39363901    REAL(wp), INTENT( IN )  ::  erel   !< Latent energy ballance   (W)
     3902    REAL(wp), INTENT( IN )  ::  ht     !< height                   (m)
    39373903    REAL(wp), INTENT( IN )  ::  int_heat  !< internal heat production (W)
     3904    REAL(wp), INTENT( IN )  ::  mbody  !< body mass                (kg)
    39383905    REAL(wp), INTENT( IN )  ::  pair   !< Air pressure             (hPa)
    39393906    REAL(wp), INTENT( IN )  ::  ta     !< Air temperature          (degree_C)
     
    39413908    REAL(wp), INTENT( IN )  ::  v      !< Wind speed               (m/s)
    39423909    REAL(wp), INTENT( IN )  ::  vpa    !< Vapor pressure           (hPa)
    3943     REAL(wp), INTENT( IN )  ::  mbody  !< body mass                (kg)
    3944     REAL(wp), INTENT( IN )  ::  ht     !< height                   (m)
    3945     REAL(wp), INTENT( IN )  ::  clo    !< clothing insulation      (clo)
    3946     REAL(wp), INTENT( IN )  ::  fcl    !< factor for surface area increase by clothing
    39473910!
    39483911!-- Output arguments:
     
    39613924!
    39623925!-- Cconstants:
    3963 !     REAL(wp), PARAMETER :: cair = 1010._wp      !< replaced by c_p
    3964     REAL(wp), PARAMETER :: cb   = 3640._wp        !<
     3926!     REAL(wp), PARAMETER :: cair = 1010.0_wp      !< replaced by c_p
     3927    REAL(wp), PARAMETER :: cb   = 3640.0_wp        !<
    39653928    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
    39663929    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
    3967 !     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
    3968     REAL(wp), PARAMETER :: food =    0._wp        !< Heat gain by food        (W)
     3930!    REAL(wp), PARAMETER :: evap = 2.42_wp * 10.0_wp **6.0_wp  !< replaced by l_v
     3931    REAL(wp), PARAMETER :: food =    0.0_wp        !< Heat gain by food        (W)
    39693932    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
    3970     REAL(wp), PARAMETER :: rob  =    1.06_wp      !< 
     3933    REAL(wp), PARAMETER :: rob  =    1.06_wp      !<
    39713934!
    39723935!-- Internal variables
    3973     REAL(wp) ::  c(0:10)        !< Core temperature array           (degree_C)
     3936    INTEGER(iwp) ::  count1     !< running index
     3937    INTEGER(iwp) ::  count3     !< running index
     3938    INTEGER(iwp) ::  j          !< running index
     3939    INTEGER(iwp) ::  i          !< running index
     3940
     3941    LOGICAL ::  skipincreasecount   !< iteration control flag
     3942
    39743943    REAL(wp) ::  cbare          !< Convection through bare skin
    39753944    REAL(wp) ::  cclo           !< Convection through clothing
     
    39823951    REAL(wp) ::  eswphy         !< sweat created by physiology
    39833952    REAL(wp) ::  eswpot         !< potential sweat evaporation
    3984     REAL(wp) ::  fec            !< 
    3985     REAL(wp) ::  hc             !< 
    3986     REAL(wp) ::  he             !< 
    3987     REAL(wp) ::  htcl           !< 
    3988     REAL(wp) ::  r1             !< 
    3989     REAL(wp) ::  r2             !< 
     3953    REAL(wp) ::  fec            !<
     3954    REAL(wp) ::  hc             !<
     3955    REAL(wp) ::  he             !<
     3956    REAL(wp) ::  htcl           !<
     3957    REAL(wp) ::  r1             !<
     3958    REAL(wp) ::  r2             !<
    39903959    REAL(wp) ::  rbare          !< Radiational loss of bare skin    (W/m²)
    3991     REAL(wp) ::  rcl            !< 
     3960    REAL(wp) ::  rcl            !<
    39923961    REAL(wp) ::  rclo           !< Radiational loss of clothing     (W/m²)
    39933962    REAL(wp) ::  rclo2          !< Longwave radiation gain or loss  (W/m²)
    39943963    REAL(wp) ::  rsum           !< Radiational loss or gain         (W/m²)
    3995     REAL(wp) ::  sw             !<
    3996 !     REAL(wp) ::  swf            !< female factor, currently unused
    3997     REAL(wp) ::  swm            !<
    3998     REAL(wp) ::  tbody          !<
    3999     REAL(wp) ::  tcore(1:7)     !<
    4000     REAL(wp) ::  vb             !<
    4001     REAL(wp) ::  vb1            !<
    4002     REAL(wp) ::  vb2            !<
    4003     REAL(wp) ::  wd             !<
    4004     REAL(wp) ::  wr             !<
    4005     REAL(wp) ::  ws             !<
    4006     REAL(wp) ::  wsum           !<
     3964    REAL(wp) ::  sw             !<
     3965!    REAL(wp) ::  swf            !< female factor, currently unused
     3966    REAL(wp) ::  swm            !<
     3967    REAL(wp) ::  tbody          !<
     3968    REAL(wp) ::  vb             !<
     3969    REAL(wp) ::  vb1            !<
     3970    REAL(wp) ::  vb2            !<
     3971    REAL(wp) ::  wd             !<
     3972    REAL(wp) ::  wr             !<
     3973    REAL(wp) ::  ws             !<
     3974    REAL(wp) ::  wsum           !<
    40073975    REAL(wp) ::  xx             !< modification step                (K)
    40083976    REAL(wp) ::  y              !< fraction of bare skin
    4009     INTEGER(iwp) ::  count1     !< running index
    4010     INTEGER(iwp) ::  count3     !< running index
    4011     INTEGER(iwp) ::  j          !< running index
    4012     INTEGER(iwp) ::  i          !< running index
    4013     LOGICAL ::  skipIncreaseCount   !< iteration control flag
     3977
     3978    REAL(wp) ::  c(0:10)        !< Core temperature array           (degree_C)
     3979    REAL(wp) ::  tcore(1:7)     !<
    40143980
    40153981!
    40163982!-- Initialize
    4017     wetsk = 0._wp  !< skin is dry everywhere on init (no non-evaporated sweat)
     3983    wetsk = 0.0_wp  !< skin is dry everywhere on init (no non-evaporated sweat)
    40183984!
    40193985!-- Set Du Bois Area for the sample person
     
    40283994!
    40293995!-- Set surface modification by clothing
    4030     facl = ( - 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp    &
    4031           * ( clo**3._wp ) ) / 100._wp
    4032     IF ( facl > 1._wp )  facl = 1._wp
     3996    facl = ( - 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp * ( clo**3.0_wp ) )    &
     3997           / 100.0_wp
     3998    IF ( facl > 1.0_wp )  facl = 1.0_wp
    40333999!
    40344000!-- Initialize heat resistences
    40354001    rcl = ( clo / 6.45_wp ) / facl
    4036     IF ( clo >= 2._wp )  y  = 1._wp
    4037     IF ( ( clo > 0.6_wp )   .AND.  ( clo < 2._wp ) )   y = ( ht - 0.2_wp ) / ht
     4002    IF ( clo >= 2.0_wp )  y  = 1.0_wp
     4003    IF ( ( clo > 0.6_wp )   .AND.  ( clo < 2.0_wp ) )   y = ( ht - 0.2_wp ) / ht
    40384004    IF ( ( clo <= 0.6_wp )  .AND.  ( clo > 0.3_wp ) )  y = 0.5_wp
    4039     IF ( ( clo <= 0.3_wp )  .AND.  ( clo > 0._wp ) )   y = 0.1_wp
    4040     r2   = adu * ( fcl - 1._wp + facl ) / ( 2._wp * 3.14_wp * ht * y )
    4041     r1   = facl * adu / ( 2._wp * 3.14_wp * ht * y )
     4005    IF ( ( clo <= 0.3_wp )  .AND.  ( clo > 0.0_wp ) )   y = 0.1_wp
     4006    r2   = adu * ( fcl - 1.0_wp + facl ) / ( 2.0_wp * 3.14_wp * ht * y )
     4007    r1   = facl * adu / ( 2.0_wp * 3.14_wp * ht * y )
    40424008    di   = r2 - r1
    40434009
     
    40464012    DO  j = 1, 7
    40474013
    4048        tsk    = 34._wp
     4014       tsk    = 34.0_wp
    40494015       count1 = 0_iwp
    4050        tcl    = ( ta + tmrt + tsk ) / 3._wp
     4016       tcl    = ( ta + tmrt + tsk ) / 3.0_wp
    40514017       count3 = 1_iwp
    4052        enbal2 = 0._wp
     4018       enbal2 = 0.0_wp
    40534019
    40544020       DO  i = 1, 100  ! allow for 100 iterations max
    4055           acl   = adu * facl + adu * ( fcl - 1._wp )
    4056           rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4._wp -             &
    4057             ( tmrt + degc_to_k )**4._wp ) * feff
     4021          acl   = adu * facl + adu * ( fcl - 1.0_wp )
     4022          rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4.0_wp -                                &
     4023                  ( tmrt + degc_to_k )**4.0_wp ) * feff
    40584024          htcl  = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
    4059           tsk   = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
     4025          tsk   = 1.0_wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
    40604026!
    40614027!--       Radiation saldo
    40624028          aeff  = adu * feff
    4063           rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
    4064             ( ( tmrt + degc_to_k )**4._wp - ( tsk + degc_to_k )**4._wp )
    4065           rclo  = feff * acl * emcl * sigma_sb *                               &
    4066             ( ( tmrt + degc_to_k )**4._wp - ( tcl + degc_to_k )**4._wp )
     4029          rbare = aeff * ( 1.0_wp - facl ) * emsk * sigma_sb *                                     &
     4030                  ( ( tmrt + degc_to_k )**4.0_wp - ( tsk + degc_to_k )**4.0_wp )
     4031          rclo  = feff * acl * emcl * sigma_sb *                                                   &
     4032                  ( ( tmrt + degc_to_k )**4.0_wp - ( tcl + degc_to_k )**4.0_wp )
    40674033          rsum  = rbare + rclo
    40684034!
    40694035!--       Convection
    4070           cbare = hc * ( ta - tsk ) * adu * ( 1._wp - facl )
     4036          cbare = hc * ( ta - tsk ) * adu * ( 1.0_wp - facl )
    40714037          cclo  = hc * ( ta - tcl ) * acl
    40724038          csum  = cbare + cclo
     
    40754041          c(0)  = int_heat + ere
    40764042          c(1)  = adu * rob * cb
    4077           c(2)  = 18._wp - 0.5_wp * tsk
     4043          c(2)  = 18.0_wp - 0.5_wp * tsk
    40784044          c(3)  = 5.28_wp * adu * c(2)
    40794045          c(4)  = 0.0208_wp * c(1)
     
    40814047          c(6)  = c(3) - c(5) - tsk * c(4)
    40824048          c(7)  = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
    4083           c(8)  = c(6) * c(6) - 4._wp * c(4) * c(7)
     4049          c(8)  = c(6) * c(6) - 4.0_wp * c(4) * c(7)
    40844050          c(9)  = 5.28_wp * adu - c(5) - c(4) * tsk
    4085           c(10) = c(9) * c(9) - 4._wp * c(4) *                                 &
    4086                   ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
    4087 
    4088           IF ( ABS( tsk - 36._wp ) < 0.00001_wp )  tsk = 36.01_wp
    4089           tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk
    4090           tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) /   &
    4091             ( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk
    4092           IF ( c(10) >= 0._wp )  THEN
    4093              tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2._wp * c(4) )
    4094              tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2._wp * c(4) )
     4051          c(10) = c(9) * c(9) - 4.0_wp * c(4) * ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
     4052
     4053          IF ( ABS( tsk - 36.0_wp ) < 0.00001_wp )  tsk = 36.01_wp
     4054          tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600.0_wp ) + tsk
     4055          tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600.0_wp ) /   &
     4056                     ( 1.0_wp + 0.5_wp * ( 34.0_wp - tsk ) ) ) + tsk
     4057          IF ( c(10) >= 0.0_wp )  THEN
     4058             tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2.0_wp * c(4) )
     4059             tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2.0_wp * c(4) )
    40954060          ENDIF
    40964061
    4097           IF ( c(8) >= 0._wp )  THEN
    4098              tcore(2) = ( - c(6) + ABS( c(8) )**0.5_wp ) / ( 2._wp * c(4) )
    4099              tcore(5) = ( - c(6) - ABS( c(8) )**0.5_wp ) / ( 2._wp * c(4) )
    4100              tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk
     4062          IF ( c(8) >= 0.0_wp )  THEN
     4063             tcore(2) = ( - c(6) + ABS( c(8) )**0.5_wp ) / ( 2.0_wp * c(4) )
     4064             tcore(5) = ( - c(6) - ABS( c(8) )**0.5_wp ) / ( 2.0_wp * c(4) )
     4065             tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1.0_wp / 40.0_wp ) + tsk
    41014066          ENDIF
    41024067!
    41034068!--       Transpiration
    41044069          tbody = 0.1_wp * tsk + 0.9_wp * tcore(j)
    4105           swm   = 304.94_wp * ( tbody - 36.6_wp ) * adu / 3600000._wp
    4106           vpts  = 6.11_wp * 10._wp**( 7.45_wp * tsk / ( 235._wp + tsk ) )
    4107 
    4108           IF ( tbody <= 36.6_wp )  swm = 0._wp  !< no need for sweating
     4070          swm   = 304.94_wp * ( tbody - 36.6_wp ) * adu / 3600000.0_wp
     4071          vpts  = 6.11_wp * 10.0_wp**( 7.45_wp * tsk / ( 235.0_wp + tsk ) )
     4072
     4073          IF ( tbody <= 36.6_wp )  swm = 0.0_wp  !< no need for sweating
    41094074
    41104075          sw = swm
    41114076          eswphy = - sw * l_v
    41124077          he     = 0.633_wp * hc / ( pair * c_p )
    4113           fec    = 1._wp / ( 1._wp + 0.92_wp * hc * rcl )
     4078          fec    = 1.0_wp / ( 1.0_wp + 0.92_wp * hc * rcl )
    41144079          eswpot = he * ( vpa - vpts ) * adu * l_v * fec
    41154080          wetsk  = eswphy / eswpot
    41164081
    4117           IF ( wetsk > 1._wp )  wetsk = 1._wp
     4082          IF ( wetsk > 1.0_wp )  wetsk = 1.0_wp
    41184083!
    41194084!--       Sweat production > evaporation?
    41204085          eswdif = eswphy - eswpot
    41214086
    4122           IF ( eswdif <= 0._wp )  esw = eswpot  !< Limit is evaporation
    4123           IF ( eswdif > 0._wp )   esw = eswphy  !< Limit is sweat production
    4124           IF ( esw  > 0._wp )     esw = 0._wp   !< Sweat can't be evaporated, no more cooling effect
     4087          IF ( eswdif <= 0.0_wp )  esw = eswpot     !< Limit is evaporation
     4088          IF ( eswdif > 0.0_wp )   esw = eswphy     !< Limit is sweat production
     4089          IF ( esw  > 0.0_wp )     esw = 0.0_wp     !< Sweat can't be evaporated, no more cooling
     4090                                                    !< effect
    41254091!
    41264092!--       Diffusion
    4127           rdsk = 0.79_wp * 10._wp**7._wp
    4128           rdcl = 0._wp
    4129           ed   = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa -     &
    4130              vpts )
     4093          rdsk = 0.79_wp * 10.0_wp**7.0_wp
     4094          rdcl = 0.0_wp
     4095          ed   = l_v / ( rdsk + rdcl ) * adu * ( 1.0_wp - wetsk ) * ( vpa - vpts )
    41314096!
    41324097!--       Max vb
    4133           vb1 = 34._wp - tsk
     4098          vb1 = 34.0_wp - tsk
    41344099          vb2 = tcore(j) - 36.6_wp
    41354100
    4136           IF ( vb2 < 0._wp )  vb2 = 0._wp
    4137           IF ( vb1 < 0._wp )  vb1 = 0._wp
    4138           vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 )
     4101          IF ( vb2 < 0.0_wp )  vb2 = 0.0_wp
     4102          IF ( vb1 < 0.0_wp )  vb1 = 0.0_wp
     4103          vb = ( 6.3_wp + 75.0_wp * vb2 ) / ( 1.0_wp + 0.5_wp * vb1 )
    41394104!
    41404105!--       Energy ballence
     
    41434108!--       Clothing temperature
    41444109          xx = 0.001_wp
    4145           IF ( count1 == 0_iwp )  xx = 1._wp
     4110          IF ( count1 == 0_iwp )  xx = 1.0_wp
    41464111          IF ( count1 == 1_iwp )  xx = 0.1_wp
    41474112          IF ( count1 == 2_iwp )  xx = 0.01_wp
    41484113          IF ( count1 == 3_iwp )  xx = 0.001_wp
    41494114
    4150           IF ( enbal > 0._wp )  tcl = tcl + xx
    4151           IF ( enbal < 0._wp )  tcl = tcl - xx
    4152 
    4153           skipIncreaseCount = .FALSE.
    4154           IF ( ( (enbal <= 0._wp )  .AND.  (enbal2 > 0._wp ) )  .OR.           &
    4155              ( ( enbal >= 0._wp )  .AND.  ( enbal2 < 0._wp ) ) )  THEN
    4156              skipIncreaseCount = .TRUE.
     4115          IF ( enbal > 0.0_wp )  tcl = tcl + xx
     4116          IF ( enbal < 0.0_wp )  tcl = tcl - xx
     4117
     4118          skipincreasecount = .FALSE.
     4119          IF ( ( (enbal <= 0.0_wp )  .AND.  (enbal2 > 0.0_wp ) )  .OR.                             &
     4120             ( ( enbal >= 0.0_wp )   .AND.  ( enbal2 < 0.0_wp ) ) )  THEN
     4121             skipincreasecount = .TRUE.
    41574122          ELSE
    41584123             enbal2 = enbal
     
    41604125          ENDIF
    41614126
    4162           IF ( ( count3 > 200_iwp )  .OR.  skipIncreaseCount )  THEN
     4127          IF ( ( count3 > 200_iwp )  .OR.  skipincreasecount )  THEN
    41634128             IF ( count1 < 3_iwp )  THEN
    41644129                count1 = count1 + 1_iwp
    4165                 enbal2 = 0._wp
     4130                enbal2 = 0.0_wp
    41664131             ELSE
    41674132                EXIT
     
    41724137       IF ( count1 == 3_iwp )  THEN
    41734138          SELECT CASE ( j )
    4174              CASE ( 2, 5)
    4175                 IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.                    &
    4176                    ( tsk <= 34.050_wp ) ) )  CYCLE
     4139             CASE ( 2, 5)
     4140                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.  ( tsk <= 34.050_wp ) ) )  CYCLE
    41774141             CASE ( 6, 1 )
    4178                 IF ( c(10) < 0._wp ) CYCLE
    4179                 IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.                    &
    4180                    ( tsk > 33.850_wp ) ) )  CYCLE
     4142                IF ( c(10) < 0.0_wp ) CYCLE
     4143                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.  ( tsk > 33.850_wp ) ) )  CYCLE
    41814144             CASE ( 3 )
    4182                 IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.                     &
    4183                    ( tsk <= 34.000_wp ) ) )  CYCLE
     4145                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.  ( tsk <= 34.000_wp ) ) )  CYCLE
    41844146             CASE ( 7 )
    4185                 IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.                     &
    4186                    ( tsk > 34.000_wp ) ) )  CYCLE
     4147                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.  ( tsk > 34.000_wp ) ) )  CYCLE
    41874148             CASE default
    41884149          END SELECT
    41894150       ENDIF
    41904151
    4191        IF ( ( j /= 4_iwp )  .AND.  ( vb >= 91._wp ) )  CYCLE
    4192        IF ( ( j == 4_iwp )  .AND.  ( vb < 89._wp ) )  CYCLE
    4193        IF ( vb > 90._wp ) vb = 90._wp
     4152       IF ( ( j /= 4_iwp )  .AND.  ( vb >= 91.0_wp ) )  CYCLE
     4153       IF ( ( j == 4_iwp )  .AND.  ( vb < 89.0_wp ) )  CYCLE
     4154       IF ( vb > 90.0_wp ) vb = 90.0_wp
    41944155!
    41954156!--    Loses by water
    4196        ws = sw * 3600._wp * 1000._wp
    4197        IF ( ws > 2000._wp )  ws = 2000._wp
    4198        wd = ed / l_v * 3600._wp * ( -1000._wp )
    4199        wr = erel / l_v * 3600._wp * ( -1000._wp )
     4157       ws = sw * 3600.0_wp * 1000.0_wp
     4158       IF ( ws > 2000.0_wp )  ws = 2000.0_wp
     4159       wd = ed / l_v * 3600.0_wp * ( -1000.0_wp )
     4160       wr = erel / l_v * 3600.0_wp * ( -1000.0_wp )
    42004161
    42014162       wsum = ws + wr + wd
     
    42054166 END SUBROUTINE heat_exch
    42064167
    4207 !------------------------------------------------------------------------------!
     4168!--------------------------------------------------------------------------------------------------!
    42084169! Description:
    42094170! ------------
    42104171!> Calculate PET
    4211 !------------------------------------------------------------------------------!
    4212  SUBROUTINE pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair,    &
    4213         rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
     4172!--------------------------------------------------------------------------------------------------!
     4173 SUBROUTINE pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair, rdcl, rdsk, rtv, ta,   &
     4174                          tcl, tsk, pet_ij, vpts, wetsk )
    42144175!
    42154176!-- Input arguments:
    4216     REAL(wp), INTENT( IN ) ::  acl   !< clothing surface area        (m²)
    4217     REAL(wp), INTENT( IN ) ::  adu   !< Du-Bois area                 (m²)
    4218     REAL(wp), INTENT( IN ) ::  esw   !< energy-loss through sweat evap. (W)
    4219     REAL(wp), INTENT( IN ) ::  facl  !< surface area extension through clothing (factor)
    4220     REAL(wp), INTENT( IN ) ::  feff  !< surface modification by posture (factor)
     4177    REAL(wp), INTENT( IN ) ::  acl       !< clothing surface area        (m²)
     4178    REAL(wp), INTENT( IN ) ::  adu       !< Du-Bois area                 (m²)
     4179    REAL(wp), INTENT( IN ) ::  esw       !< energy-loss through sweat evap. (W)
     4180    REAL(wp), INTENT( IN ) ::  facl      !< surface area extension through clothing (factor)
     4181    REAL(wp), INTENT( IN ) ::  feff      !< surface modification by posture (factor)
    42214182    REAL(wp), INTENT( IN ) ::  int_heat  !< internal heat production (W)
    4222     REAL(wp), INTENT( IN ) ::  pair  !< air pressure                 (hPa)
    4223     REAL(wp), INTENT( IN ) ::  rdcl  !< diffusion resistence of clothing (factor)
    4224     REAL(wp), INTENT( IN ) ::  rdsk  !< diffusion resistence of skin (factor)
    4225     REAL(wp), INTENT( IN ) ::  rtv   !< respiratory volume
    4226     REAL(wp), INTENT( IN ) ::  ta    !< air temperature              (degree_C)
    4227     REAL(wp), INTENT( IN ) ::  tcl   !< clothing temperature         (degree_C)
    4228     REAL(wp), INTENT( IN ) ::  tsk   !< skin temperature             (degree_C)
    4229     REAL(wp), INTENT( IN ) ::  vpts  !< sat. vapor pressure over skin (hPa)
    4230     REAL(wp), INTENT( IN ) ::  wetsk !< fraction of wet skin (dimensionless)
     4183    REAL(wp), INTENT( IN ) ::  pair      !< air pressure                 (hPa)
     4184    REAL(wp), INTENT( IN ) ::  rdcl      !< diffusion resistence of clothing (factor)
     4185    REAL(wp), INTENT( IN ) ::  rdsk      !< diffusion resistence of skin (factor)
     4186    REAL(wp), INTENT( IN ) ::  rtv       !< respiratory volume
     4187    REAL(wp), INTENT( IN ) ::  ta        !< air temperature              (degree_C)
     4188    REAL(wp), INTENT( IN ) ::  tcl       !< clothing temperature         (degree_C)
     4189    REAL(wp), INTENT( IN ) ::  tsk       !< skin temperature             (degree_C)
     4190    REAL(wp), INTENT( IN ) ::  vpts      !< sat. vapor pressure over skin (hPa)
     4191    REAL(wp), INTENT( IN ) ::  wetsk     !< fraction of wet skin (dimensionless)
    42314192!
    42324193!-- Output arguments: